text:
# terminals
isidentifier(x::EXPR) = headof(x) === :IDENTIFIER || headof(x) === :NONSTDIDENTIFIER
isoperator(x::EXPR) = headof(x) === :OPERATOR
isnonstdid(x::EXPR) = headof(x) === :NONSTDIDENTIFIER
iskeyword(x::EXPR) = headof(x) in (:ABSTRACT, :BAREMODULE, :BEGIN, :BREAK, :CATCH, :CONST, :CONTINUE, :DO, :ELSE, :ELSEIF, :END, :EXPORT, :FINALLY, :FOR, :FUNCTION, :GLOBAL, :IF, :IMPORT, :importall, :LET, :LOCAL, :MACRO, :MODULE, :MUTABLE, :NEW, :OUTER, :PRIMITIVE, :QUOTE, :RETURN, :STRUCT, :TRY, :TYPE, :USING, :WHILE)
isliteral(x::EXPR) = isstringliteral(x) || iscmd(x) || ischar(x) || headof(x) in (:INTEGER, :BININT, :HEXINT, :OCTINT, :FLOAT,  :NOTHING, :TRUE, :FALSE)
ispunctuation(x::EXPR) = is_comma(x) || is_lparen(x) || is_rparen(x) || is_lsquare(x) || is_rsquare(x) || is_lbrace(x) || is_rbrace(x) || headof(x) === :ATSIGN  || headof(x) === :DOT
isstringliteral(x) = headof(x) === :STRING || headof(x) === :TRIPLESTRING
isstring(x) = headof(x) === :string || isstringliteral(x)
iscmd(x) = headof(x) === :CMD || headof(x) === :TRIPLECMD
ischar(x) = headof(x) === :CHAR
isinteger(x) = headof(x) === :INTEGER
isfloat(x) = headof(x) === :FLOAT
isnumber(x) = isinteger(x) || isfloat(x)
is_nothing(x) = headof(x) === :NOTHING
is_either_id_op_interp(x::EXPR) = isidentifier(x) || isoperator(x) || isinterpolant(x)
is_id_or_macroname(x::EXPR) = isidentifier(x) || ismacroname(x)

# expressions
iscall(x::EXPR) = headof(x) === :call
isunarycall(x::EXPR) = (headof(x) === :call && length(x) == 2 && (isoperator(x.args[1]) || isoperator(x.args[2])))
isunarysyntax(x::EXPR) =  (isoperator(x.head) && length(x.args) == 1)
isbinarycall(x::EXPR) = headof(x) === :call && length(x) == 3 && isoperator(x.args[1])
isbinarycall(x::EXPR, op) = headof(x) === :call && length(x) == 3 && isoperator(x.args[1]) && valof(x.args[1]) == op
isbinarysyntax(x::EXPR) =  (isoperator(x.head) && length(x.args) == 2)
ischainedcall(x::EXPR) = headof(x) === :call && isoperator(x.args[1]) && (valof(x.args[1]) == "+" || valof(x.args[1]) == "*") && hastrivia(x) && isoperator(first(x.trivia))
iswhere(x::EXPR) = headof(x) === :where
istuple(x::EXPR) = headof(x) === :tuple
ismacrocall(x::EXPR) = headof(x) === :macrocall
ismacroname(x::EXPR) = (headof(x) === :NONSTDIDENTIFIER && startswith(valof(first(x.args)), "@") || (isidentifier(x) && valof(x) !== nothing && !isempty(valof(x)) && first(valof(x)) == '@')) || (is_getfield_w_quotenode(x) && (ismacroname(unquotenode(rhs_getfield(x))) || ismacroname(x.args[1])))

iskwarg(x::EXPR) = headof(x) === :kw
isparameters(x::EXPR) = headof(x) === :parameters
iscurly(x::EXPR) = headof(x) === :curly
isbracketed(x::EXPR) = headof(x) === :brackets

isassignment(x::EXPR) = isbinarysyntax(x) && valof(x.head) == "="
isdeclaration(x::EXPR) = isbinarysyntax(x) && valof(x.head) == "::"
isinterpolant(x::EXPR) = isunarysyntax(x) && valof(x.head) == "\$"
issplat(x::EXPR) = isunarysyntax(x) && valof(x.head) == "..."

isbeginorblock(x::EXPR) = headof(x) === :begin || headof(unwrapbracket(x)) == :block

"""
    is_getfield(x::EXPR)

Is this an expression of the form `a.b`.
"""
is_getfield(x::EXPR) = x.head isa EXPR && isoperator(x.head) && valof(x.head) == "." && length(x.args) == 2
is_getfield_w_quotenode(x) = is_getfield(x) && headof(x.args[2]) === :quotenode && length(x.args[2].args) > 0
rhs_getfield(x) = x.args[2]
unquotenode(x) = x.args[1]

function get_rhs_of_getfield(ret::EXPR)
    if headof(ret.args[2]) === :quotenode || headof(ret.args[2]) === :quote
        ret.args[2].args[1]
    else
        ret.args[2]
    end
end

"""
    is_wrapped_assignment(x::EXPR)

Is `x` an assignment expression, ignoring any surrounding parentheses.
"""
is_wrapped_assignment(x::EXPR) = isassignment(x) || (isbracketed(x) && is_wrapped_assignment(x.args[1]))

function is_func_call(x::EXPR)
    if isoperator(x.head) && !issplat(x)
        if length(x.args) == 2
            return is_decl(x.head) && is_func_call(x.args[1])
        elseif length(x.args) == 1
            return !(is_exor(x.head) || is_decl(x.head))
        end
    elseif x.head === :call
        return true
    elseif x.head === :where || isbracketed(x)
        return is_func_call(x.args[1])
    end
    return false
end

fcall_name(x::EXPR) = iscall(x) && length(x.args) > 0 && valof(x.args[1])
hasparent(x::EXPR) = parentof(x) isa EXPR

# OPERATOR
is_approx(x::EXPR) = isoperator(x) && valof(x) == "~"
is_exor(x) = isoperator(x) && valof(x) == "\$"
is_decl(x) = isoperator(x) && valof(x) == "::"
is_issubt(x) = isoperator(x) && valof(x) == "<:"
is_issupt(x) = isoperator(x) && valof(x) == ">:"
is_and(x) = isoperator(x) && valof(x) == "&"
is_not(x) = isoperator(x) && valof(x) == "!"
is_plus(x) = isoperator(x) && valof(x) == "+"
is_minus(x) = isoperator(x) && valof(x) == "-"
is_star(x) = isoperator(x) && valof(x) == "*"
is_eq(x) = isoperator(x) && valof(x) == "="
is_dot(x) = isoperator(x) && valof(x) == "."
is_ddot(x) = isoperator(x) && valof(x) == ".."
is_dddot(x) = isoperator(x) && valof(x) == "..."
is_pairarrow(x) = isoperator(x) && valof(x) == "=>"
is_in(x) = isoperator(x) && valof(x) == "in"
is_elof(x) = isoperator(x) && valof(x) == "∈"
is_colon(x) = isoperator(x) && valof(x) == ":"
is_prime(x) = isoperator(x) && maybe_strip_suffix(valof(x)) == "'"
is_cond(x) = isoperator(x) && valof(x) == "?"
is_where(x) = isoperator(x) && valof(x) == "where"
is_anon_func(x) = isoperator(x) && valof(x) == "->"

is_comma(x) = headof(x) === :COMMA
is_lparen(x) = headof(x) === :LPAREN
is_rparen(x) = headof(x) === :RPAREN
is_lbrace(x) = headof(x) === :LBRACE
is_rbrace(x) = headof(x) === :RBRACE
is_lsquare(x) = headof(x) === :LSQUARE
is_rsquare(x) = headof(x) === :RSQUARE

# KEYWORD
is_if(x) = iskeyword(x) && headof(x) === :IF
is_import(x) = iskeyword(x) && headof(x) === :IMPORT

# Literals

issubtypedecl(x::EXPR) = isoperator(x.head) && valof(x.head) == "<:"

rem_subtype(x::EXPR) = issubtypedecl(x) ? x.args[1] : x
rem_decl(x::EXPR) = isdeclaration(x) ? x.args[1] : x
rem_curly(x::EXPR) = headof(x) === :curly ? x.args[1] : x
rem_call(x::EXPR) = headof(x) === :call ? x.args[1] : x
rem_where(x::EXPR) = iswhere(x) ? x.args[1] : x
rem_wheres(x::EXPR) = iswhere(x) ? rem_wheres(x.args[1]) : x
rem_where_subtype(x::EXPR) = (iswhere(x) || issubtypedecl(x)) ? x.args[1] : x
rem_wheres_subtypes(x::EXPR) = (iswhere(x) || issubtypedecl(x)) ? rem_wheres_subtypes(x.args[1]) : x
rem_where_decl(x::EXPR) = (iswhere(x) || isdeclaration(x)) ? x.args[1] : x
rem_wheres_decls(x::EXPR) = (iswhere(x) || isdeclaration(x)) ? rem_wheres_decls(x.args[1]) : x
rem_invis(x::EXPR) = isbracketed(x) ? rem_invis(x.args[1]) : x
rem_dddot(x::EXPR) = issplat(x) ? x.args[1] : x
const rem_splat = rem_dddot
rem_kw(x::EXPR) = headof(x) === :kw ? x.args[1] : x
unwrapbracket(x::EXPR) = isbracketed(x) ? unwrapbracket(x.args[1]) : x

is_some_call(x) = headof(x) === :call || isunarycall(x)
is_eventually_some_call(x) = is_some_call(x) || ((isdeclaration(x) || iswhere(x)) && is_eventually_some_call(x.args[1]))

defines_function(x::EXPR) = headof(x) === :function || (isassignment(x) && is_eventually_some_call(x.args[1]))
defines_macro(x) = headof(x) == :macro
defines_datatype(x) = defines_struct(x) || defines_abstract(x) || defines_primitive(x)
defines_struct(x) = headof(x) === :struct
defines_mutable(x) = defines_struct(x) && x.args[1].head == :TRUE
defines_abstract(x) = headof(x) === :abstract
defines_primitive(x) = headof(x) === :primitive
defines_module(x) = headof(x) === :module
defines_anon_function(x) = isoperator(x.head) && valof(x.head) == "->"

has_sig(x::EXPR) = defines_datatype(x) || defines_function(x) || defines_macro(x) || defines_anon_function(x)

"""
    get_sig(x)

Returns the full signature of function, macro and datatype definitions.
Should only be called when has_sig(x) == true.
"""
function get_sig(x::EXPR)
    if headof(x) isa EXPR # headof(headof(x)) === :OPERATOR valof(headof(x)) == "="
        return x.args[1]
    elseif headof(x) === :struct || headof(x) === :mutable
        return x.args[2]
    elseif headof(x) === :abstract || headof(x) === :primitive || headof(x) === :function || headof(x) === :macro
        return x.args[1]
    end
end

function get_name(x::EXPR)
    if defines_datatype(x)
        sig = get_sig(x)
        sig = rem_subtype(sig)
        sig = rem_wheres(sig)
        sig = rem_subtype(sig)
        sig = rem_curly(sig)
    elseif defines_module(x)
        sig = x.args[2]
    elseif defines_function(x) || defines_macro(x)
        sig = get_sig(x)
        sig = rem_wheres(sig)
        sig = rem_decl(sig)
        sig = rem_call(sig)
        sig = rem_curly(sig)
        sig = rem_invis(sig)
        # if isbinarysyntax(sig) && is_dot(sig.head)
        if is_getfield_w_quotenode(sig)
            sig = sig.args[2].args[1]
        end
        return sig
    elseif is_getfield_w_quotenode(x)
        return x.args[2].args[1]
    elseif isbinarycall(x)
        sig = x.args[1]
        if isunarycall(sig)
            return get_name(sig.args[1])
        end
        sig = rem_wheres(sig)
        sig = rem_decl(sig)
        sig = rem_call(sig)
        sig = rem_curly(sig)
        sig = rem_invis(sig)
        return get_name(sig)
    elseif isbinarysyntax(x) && valof(x.head) == "<:"
        return get_name(x.args[1])
    elseif isunarysyntax(x) && valof(x.head) == "..."
        return get_name(x.args[1])
    else
        sig = x
        if isunarycall(sig)
            sig = sig.args[1]
        end
        sig = rem_wheres(sig)
        sig = rem_decl(sig)
        sig = rem_call(sig)
        sig = rem_curly(sig)
        rem_invis(sig)
    end
end

function get_arg_name(arg::EXPR)
    arg = rem_kw(arg)
    arg = rem_dddot(arg)
    arg = rem_where(arg)
    arg = rem_decl(arg)
    arg = rem_subtype(arg)
    arg = rem_curly(arg)
    rem_invis(arg)
end


text:
#=
"TwiceOneish" - a play on words for "bimonoidal". The package focuses on structures
with both a direct sum and tensor product, such as that seen in a symmetric bimonoidal category.
We focus on simple Cartesian structures.
=#

module TwiceOneish

import Base: size, getindex, setindex!, IndexStyle, IndexCartesian, tail,
    @propagate_inbounds, @_inline_meta

export ⊕, DirectSumArray

# Multidimensional direct sums are modelled as nested structures which
# concatenate their outermost dimensions

struct DirectSumArray{T, N} <: AbstractArray{T, N}
    a::AbstractArray{T, N}
    b::AbstractArray{T, N}

    function DirectSumArray{T, N}(a::AbstractArray{<:T, N}, b::AbstractArray{<:T, N}) where {T, N}
        @_inline_meta
        @boundscheck check_leading_sizes(size(a), size(b))
        return new{T,N}(a, b)
    end
end

const DirectSumVector{T} = DirectSumArray{T, 1}
const DirectSumMatrix{T} = DirectSumArray{T, 2}

@inline check_leading_sizes(::Tuple{}, ::Tuple{}) = nothing
@inline function check_leading_sizes(sa::Tuple, sb::Tuple)
    if sa[1] == sb[1]
        check_leading_sizes(tail(sa), tail(sb))
    else
        throw(DimensionMismatch("DirectSumArray expected leading dimensions to match"))
    end
end

DirectSumArray(a::AbstractArray{T1, N}, b::AbstractArray{T2, N}) where {T1, T2, N} = DirectSumArray{typejoin(T1, T2), N}(a, b)

# AbstractArray interface
IndexStyle(::Union{DSA, Type{DSA}} where DSA <: DirectSumArray) = IndexCartesian()

@inline size(dsa::DirectSumArray) = _size_concat((), size(dsa.a), size(dsa.b))
@inline _size_concat(out::Tuple, ::Tuple{}, ::Tuple{}) = out
@inline function _size_concat(out::Tuple, sa::Tuple, sb::Tuple)
    _size_concat((out..., sa[1] + sb[1]), tail(sa), tail(sb))
end

@propagate_inbounds function getindex(dsa::DirectSumArray{<:Any, N}, inds::Vararg{Int, N}) where {N}
    n = size(dsa.a, N)
    if inds[end] <= n
        return dsa.a[inds...]
    else
        return dsa.b[subtract_last(inds, n)...]
    end
end

@propagate_inbounds function setindex!(dsa::DirectSumArray{<:Any, N}, v, inds::Vararg{Int, N}) where {N}
    n = size(dsa.a, N)
    if inds[end] <= n
        return dsa.a[inds...] = v
    else
        return dsa.b[subtract_last(inds, n)] = v
    end
end

@inline subtract_last(inds::Tuple{Vararg{Int}}, n::Int) = _subtract_last((), n, inds...)
@inline _subtract_last(out, n, i) = (out..., i-n)
@inline _subtract_last(out, n, i, inds...) = _subtract_last((out..., i), n, inds...)

# Operator ⊕ notation
directsumcat(a::AbstractArray{<:Any, N}, b::AbstractArray{<:Any, N}) where {N} = DirectSumArray(a, b)
const ⊕ = directsumcat

end # module

ext: jl

max_stars_repo_name: andyferris/TwiceOneish


text:
using Delaunay
using Test

@testset "Delaunay.jl" begin
    # Write your tests here.
end


text:
using ForTheBadge
using Test

testdir = "$(@__DIR__)/../assets/tests"

ispath(testdir) && rm(testdir, recursive = true)
ispath(testdir) || mkpath(testdir)

cd(testdir)

@testset "create badges" begin
    badge("FOR THE", "BADGE")

    @test true
end

print("Created 1 badge.  Tests over for now!")


text:
##################################################################################################
## Read length
##################################################################################################
"""

    `fastq_len_store(FASTQ)`

    Store length of reads in FASTQ file.

    # Examples
    ```julia-repl
    julia> CpelNano.fastq_len_store(fastq)
    ```
"""
function fastq_len_store(fastq::String, outdir::String)::Nothing

    # Get prefix
    pref = split(basename(fastq), ".fastq")[1]

    # Loop over records in FASTQ file
    lens = Vector{Int64}()
    open(FASTQ.Reader, fastq) do reader
        for record in reader
            # Get lengths
            push!(lens, length(FASTQ.sequence(record)))
        end
    end

    # Write lengths
    print_log("Storing lengths...")
    open("$(outdir)/$(pref)_read_lens.txt", "w") do io
        writedlm(io, lens)
    end

    # Return nothing
    return nothing

end
"""
    `gen_read_len(NREADS)`

    Returns a vector of read lengths following an exponential distribution obtained performing
    ML estimation to read length of dataset from:

        https://github.com/nanopore-wgs-consortium/NA12878/blob/master/Genome.md

    # Examples
    ```julia-repl
    julia> CpelNano.gen_read_len(10)
    ```
"""
function gen_read_len(n_reads::Int64)::Vector{Int64}

    # Return read lengths
    return Int.(round.(rand(Exponential(8484.84), n_reads)))

end
##################################################################################################
## Generate Nucleotide Sequences
##################################################################################################
"""

    `gen_seq(seq_len,n)`

    Generates a random sequence of length `seq_len` with `n` CpG sites.

    # Examples
    ```julia-repl
    julia> CpelNano.gen_seq(10,2)
    "CCGGGGCGGA"
    ```
"""
function gen_seq(seq_len::Int64, n::Int64)::String

    # Create random sequence
    seq = randstring("ACGT", seq_len)

    # Get rid of coincidental CpGs
    seq = replace(seq, r"CG" => (c) -> "C" * randstring("AT", 1))

    # Determine spacing required
    spacing = floor(Int64, seq_len / n)
    spacing >= 2 || return ""

    # Add equally spaced CG's
    seq = split(seq, "")
    ind = floor(Int64, spacing / 2):spacing:seq_len
    for i in ind
        seq[i:(i + 1)] = ["C","G"]
    end

    # Return sequence with n CpG sites
    return join(seq)

end
"""

    `gen_read(seq,read_len)`

    Returns a random read of length `read_len` from sequence `seq`.

    # Examples
    ```julia-repl
    julia> CpelNano.gen_read("CCGGGGCGGA",5)
    (5,"GGCGG")
    ```
"""
function gen_read(chr::String, fasta::String, read_len::Int64)::Tuple{Int64,String}

    # Get position
    fa_record = CpelNano.get_chr_fa_rec(chr, fasta)
    chrstart = fa_record.sequence[1]
    chrend = fa_record.sequence[length(fa_record.sequence)];
    pos = rand(chrstart:(chrend - read_len + 1))

    # Extract read
    seq = convert(String, FASTA.sequence(fa_record, pos:(pos + read_len - 1)));

    # Return tuple
    return pos, seq

end
"""

    `count_cpg(seq)`

    Count the number of CGs in sequence `seq`.

    # Examples
    ```julia-repl
    julia> CpelNano.count_cpg("CCGGGGCGGA")
    2
    ```
"""
function count_cpg(seq::String)::Int64
    ind = minimum.(findall(r"CG", seq))
    return length(ind)
end
"""

    `count_mpg(seq)`

    Count the number of MGs in sequence `seq`.

    # Examples
    ```julia-repl
    julia> CpelNano.count_mpg("CMGGGGCGGA")
    1
    ```
"""
function count_mpg(seq::String)::Int64
    ind = minimum.(findall(r"MG", seq))
    return length(ind)
end
"""

    `meth_seq(seq,methyl_vec)`

    Converts a CG in `seq` to MG based on its methylation status from the corresponding index in `methyl_vec`.

    # Examples
    ```julia-repl
    julia> CpelNano.meth_seq("CCGGGGCGGA",CpelNano.metro_hast_inst(2,10000,0.5,0.5))
    "CMGGGGMGGA"
    ```
"""
function meth_seq(seq::String, methyl_vec::Vector{Int8})::String

    # Get number of CpG sites and spacing
    n = length(methyl_vec)
    spacing = floor(Int64, length(seq) / n)

    # Find indices of C's in CG context
    ind = minimum.(findall(r"CG", seq))
    length(ind) == n || return ""

    # Modify sequence
    seq = split(seq, "")
    seq[ind[methyl_vec .== Int8(1)]] .= "M"

    # Return modified sequence
    return join(seq)

end
"""
    `aug_mod_file(reg_file,chr_name,fasta)`

    Augments model file. If analysis region is in `reg_file`,uses parameter vector from `reg_file`. Otherwise,
    randomly samples parameter vector from analysis regions in `reg_file`. This function produces a text file
    with suffix `.aug`.

    # Examples
    ```julia-repl
    julia> CpelNano.aug_mod_file(reg_file,chr_name,fasta)
    ```
"""
function aug_mod_file(reg_file::String, chr_name::String, fasta::String)::Nothing

    # Read in model file
    reg_lines = readdlm(reg_file)

    # Get FASTA record for chr
    fa_record = get_chr_fa_rec(chr_name, fasta)

    # Initialize vectors to store parameter vector and partitions from model file
    model_part = Vector{UnitRange{Int64}}()
    ϕhat_1 = Vector{Float64}()
    ϕhat_2 = Vector{Float64}()
    ϕhat_3 = Vector{Float64}()

    # Read in parameter vectors and partitions from model files
    for i in 1:size(reg_lines)[1]
        chrst = reg_lines[i,2]
        chrend = reg_lines[i,3]
        ϕhat = parse.(Float64, split(reg_lines[i,4], ','));
        push!(model_part, chrst:chrend)
        push!(ϕhat_1, ϕhat[1])
        push!(ϕhat_2, ϕhat[2])
        push!(ϕhat_3, ϕhat[3])
    end

    # Obtain full set of partitions from whole chromosome
    config = CpelNano.CpelNanoConfig()
    chr_part = CpelNano.get_chr_part(chr_name, config, fasta)

    # Initialize vector to store analysis regions
    model_regs = Vector{CpelNano.RegStruct}()

    # Iterate through every partition in the chromosome
    for i in 1:length(chr_part)

        # Check if partition is in model file
        modelpart_index = findfirst(x -> x == chr_part[i], model_part)

        # If partition is not in model file,then randomly sample parameter vector from model file
        if isnothing(modelpart_index)
            # Randomly sample parameters
            new_ϕhat_1 = ϕhat_1[rand(1:length(ϕhat_1))]
            new_ϕhat_2 = ϕhat_2[rand(1:length(ϕhat_2))]
            new_ϕhat_3 = ϕhat_2[rand(1:length(ϕhat_3))]
        else
            # Keep the estimated vectors
            new_ϕhat_1 = ϕhat_1[modelpart_index]
            new_ϕhat_2 = ϕhat_2[modelpart_index]
            new_ϕhat_3 = ϕhat_3[modelpart_index]
        end

        # Construct RegStruct object for this analysis region
        reg_data = CpelNano.RegStruct()
        reg_data.chr = chr_name
        reg_data.chrst = minimum(chr_part[i])
        reg_data.chrend = maximum(chr_part[i])
        reg_data.ϕhat = [new_ϕhat_1,new_ϕhat_2,new_ϕhat_3]

        # Get group info assuming singletons
        CpelNano.get_grp_info!(reg_data, fa_record, 1)

        # Set processed flag as true to store
        reg_data.proc = true

        # Add RegStruct to list of RegStructs for the chromosome
        push!(model_regs, reg_data)

    end

    # Store
    write_output_ϕ("$(reg_file).aug", model_regs)

    # Return nothing
    return nothing

end
"""
    `reads_reg_chr(reg_file,chr_name,fasta,fa_record)`

    Returns a list of RegStructs for every partition in chromosome given `reg_file` containing all
    estimation regions in the chromosome. WARNING: this function assumes there is an existing model
    for each single region.

    # Examples
    ```julia-repl
    julia> model_regs = CpelNano.read_reg_chr(reg_file,chr_name,fasta,fa_record)
    ```
"""
function read_reg_chr(reg_file::String, chr_name::String, fasta::String, fa_record::FASTA.Record)::Vector{CpelNano.RegStruct}

    # Read in model file
    reg_lines = readdlm(reg_file)

    # Initialize vectors to store parameter vector and partitions from model file
    model_part = Vector{UnitRange{Int64}}()
    ϕhat_1 = Vector{Float64}()
    ϕhat_2 = Vector{Float64}()
    ϕhat_3 = Vector{Float64}()

    # Read in parameter vectors and partitions from model files
    for i in 1:size(reg_lines)[1]
        chrst = reg_lines[i,2]
        chrend = reg_lines[i,3]
        ϕhat = parse.(Float64, split(reg_lines[i,4], ','));
        push!(model_part, chrst:chrend)
        push!(ϕhat_1, ϕhat[1])
        push!(ϕhat_2, ϕhat[2])
        push!(ϕhat_3, ϕhat[3])
    end

    # Initialize vector to store analysis regions
    model_regs = Vector{CpelNano.RegStruct}()

    # Obtain full set of partitions from whole chromosome
    config = CpelNano.CpelNanoConfig();
    chr_part = CpelNano.get_chr_part(chr_name, config, fasta);

    # Iterate through every partition in the chromosome
    for i in 1:length(chr_part)

        # Get model index
        modelpart_index = findfirst(x -> x == chr_part[i], model_part)

        # Get corresponding ϕ
        new_ϕhat_1 = ϕhat_1[modelpart_index]
        new_ϕhat_2 = ϕhat_2[modelpart_index]
        new_ϕhat_3 = ϕhat_3[modelpart_index]

        # Construct RegStruct object for this analysis region
        reg_data = CpelNano.RegStruct()
        reg_data.chr = chr_name
        reg_data.chrst = minimum(chr_part[i])
        reg_data.chrend = maximum(chr_part[i])
        reg_data.ϕhat = [new_ϕhat_1,new_ϕhat_2,new_ϕhat_3]

        # Get group info assuming singletons
        CpelNano.get_grp_info!(reg_data, fa_record, 1)

        # Add RegStruct to list of RegStructs for the chromosome
        push!(model_regs, reg_data)

    end

     # Return list of RegStructs in chromosome in same order as chromosome partitions
    return model_regs

end
"""
    `linearsearch_readstart(startpos,parts)`

    Find index of analysis region corresponding to start of read `startpos` given list of chromosome partitions `parts`

    # Examples
    ```julia-repl
    julia> CpelNano.linearsearch_readstart(14996959,CpelNano.get_chr_part("chr22",CpelNano.CpelNanoConfig(),"/Users/sandeepk/Downloads/chr22.fa"))
    5000
    ```
"""
function linearsearch_readstart(startpos::Int64, parts::Vector{UnitRange{Int64}})::Int64
    for i in 1:length(parts)
        chrst = parts[i][1]
        chrend = parts[i][length(parts[i])]
        if startpos >= chrst && startpos <= chrend
            return i
        end
    end
    return -1
end
"""
    `linearsearch_readend(startpos,idx,parts)`

    Find index of analysis region corresponding to end of read `startpos` given list of chromosome partitions `parts`
    and index of analysis region corresponding to start of read `idx`

    # Examples
    ```julia-repl
    julia> CpelNano.linearsearch_readend(20999958,5000,CpelNano.get_chr_part("chr22",CpelNano.CpelNanoConfig(),"/Users/sandeepk/Downloads/chr22.fa"))
    7001
    ```
"""
function linearsearch_readend(endpos::Int64, idx::Int64, parts::Vector{UnitRange{Int64}})::Int64
    for i in idx:length(parts)
        chrst = parts[i][1]
        chrend = parts[i][length(parts[i])]
        if endpos >= chrst && endpos <= chrend
            return i
        end
    end
    return -1
end
"""
    `get_methyl_vector(rd_st,rd_end,startidx,endidx,model_regs)`

    Obtain methylation vector for a given read given start of read `rd_st`,end of read `rd_end`,
    index of analysis region corresponding to start of read `startidx`,index of analysis region
    corresponding to end of read `endidx`,and list of RegStructs in chromosome `model_regs`

    # Examples
    ```julia-repl
    julia> CpelNano.get_methyl_vector(rd_st,rd_end,startidx,endidx,model_regs)
    ```
"""
function get_methyl_vector(rd_st::Int64, rd_end::Int64, startidx::Int64, endidx::Int64, model_regs::Vector{RegStruct})::Tuple{Vector{Int8},Vector{Int64}}

    # Initialize
    x_out = Vector{Int8}()
    cpg_pos = Vector{Int64}()

    # Concatenate all αs and βs
    @inbounds for i in startidx:endidx

        # Get region struct
        rs = model_regs[i]

        # If no CpG sites, skip
        length(rs.cpg_pos) > 0 || continue

        # Get parameter vectors
        α, β = CpelNano.get_αβ_from_ϕ(rs.ϕhat, rs)

        # Generate realization x
        x = length(rs.cpg_pos) == 1 ? [Int8(1)] : CpelNano.gen_x_mc(α, β)
        x_out = vcat(x_out, x)

        # Add CpG sites
        cpg_pos = vcat(cpg_pos, rs.cpg_pos)

    end

    # Return nothing
    length(cpg_pos) > 0 || return [], []

    # Keep only overlapping part
    fst_ind = findfirst(x -> x >= rd_st, cpg_pos)
    isnothing(fst_ind) && return [], []
    lst_ind = findlast(x -> x <= rd_end, cpg_pos)
    isnothing(lst_ind) && return [], []
    cpg_pos = cpg_pos[fst_ind:lst_ind]
    x_out = x_out[fst_ind:lst_ind]

    # Return x & cpg positions
    return x_out, cpg_pos

end
##################################################################################################
## Markov chain sampler of fine grain model
##################################################################################################
"""
    `get_log_g_x(n,α,β)`

    Compute function log(g_n(x_n)).

    # Examples
    ```julia-repl using CpelNano
    julia> CpelNano.get_log_g_x(1,zeros(10),zeros(9))
    (256.0,256.0)
    ```
"""
function get_log_g_x(n::Int64, α::Vector{Float64}, β::Vector{Float64})::NTuple{2,Float64}

    # Get size
    N = length(α)

    # If n=N return 0 (log(1))
    n == N && return 0.0, 0.0

    # Define u_n for xn = ±1
    unm1 = [-(α[n + 1] - β[n]) / 2.0; +(α[n + 1] - β[n]) / 2.0]
    unp1 = [-(α[n + 1] + β[n]) / 2.0; +(α[n + 1] + β[n]) / 2.0]

    # Define u_N
    uN = [-α[end] / 2.0; α[end] / 2.0]

    # If longer gap...
    n == (N - 1) && return log_vec_vec_mult(unm1, uN), log_vec_vec_mult(unp1, uN)

    # Define Wn
    Wnm1 = [-(α[n + 1] - β[n]) / 2.0 + β[n + 1] - α[n + 2] / 2.0 -(α[n + 1] - β[n]) / 2.0 - β[n + 1] + α[n + 2] / 2.0;
            +(α[n + 1] - β[n]) / 2.0 - β[n + 1] - α[n + 2] / 2.0 +(α[n + 1] - β[n]) / 2.0 + β[n + 1] + α[n + 2] / 2.0]
    Wnp1 = [-(α[n + 1] + β[n]) / 2.0 + β[n + 1] - α[n + 2] / 2.0 -(α[n + 1] + β[n]) / 2.0 - β[n + 1] + α[n + 2] / 2.0;
            +(α[n + 1] + β[n]) / 2.0 - β[n + 1] - α[n + 2] / 2.0 +(α[n + 1] + β[n]) / 2.0 + β[n + 1] + α[n + 2] / 2.0]

    # If longer gap...
    n == (N - 2) && return log_vec_vec_mult(log_vec_mat_mult(unm1, Wnm1), uN), log_vec_vec_mult(log_vec_mat_mult(unp1, Wnp1), uN)

    # Define Wn for n'>n+1
    Wns = Vector{Array{Float64,2}}()
    @inbounds for l = (n + 2):(N - 1)
        W = [-α[l] / 2.0 + β[l] - α[l + 1] / 2.0 -α[l] / 2.0 - β[l] + α[l + 1] / 2.0;
             +α[l] / 2.0 - β[l] - α[l + 1] / 2.0 +α[l] / 2.0 + β[l] + α[l + 1] / 2.0]
        push!(Wns, W)
    end

    # Return log g_l(x_l) for x_l∈{-1,+1}
    Wns_prod = mult_log_mats(Wns)
    gl_m1 = log_vec_vec_mult(log_vec_mat_mult(log_vec_mat_mult(unm1, Wnm1), Wns_prod), uN)
    gl_p1 = log_vec_vec_mult(log_vec_mat_mult(log_vec_mat_mult(unp1, Wnp1), Wns_prod), uN)

    # return log unm1'*Wnm1*prod(Wns)*uN,log unp1'*Wnp1*prod(Wns)*uN
    return gl_m1, gl_p1

end
"""
    `gen_x_mc(α,β)`

    Generate an realization of X=[X1,…,XN] using α_n and β_n,or X̄=[X̄1,…,X̄L] if using α_l and β_l.

    # Examples
    ```julia-repl
    julia> CpelNano.gen_x_mc(ones(2),ones(1))
    2-element Array{Int8,1}:
     1
     1
    ```
"""
function gen_x_mc(α::Vector{Float64}, β::Vector{Float64})::Vector{Int8}

    # Size of x
    N = length(α)

    # Initialize x
    x = zeros(Int8, N)

    # Sample first CpG site (y1/(y1+y2))
    g = get_log_g_x(1, α, β)
    log_y1 = g[1] - α[1]
    log_y2 = g[2] + α[1]
    c = max(log_y1, log_y2)
    log_den = c + log(exp(log_y1 - c) + exp(log_y2 - c))
    p = exp(log_y2 - log_den)
    x[1] = rand() < p ? Int8(1) : Int8(-1)
        # print_log("g: $(g)")
        # print_log("p: $(p)")
        # readline()

    # Loop over all CpG sites
    @inbounds for n in 2:(N - 1)

        # Get g
        g = get_log_g_x(n, α, β)
            # g = get_g_x(n,α,β)

        # Get transition probability
        log_y1 = g[1] - α[n] - β[n - 1] * x[n - 1]
        log_y2 = g[2] + α[n] + β[n - 1] * x[n - 1]
        c = max(log_y1, log_y2)
        log_den = c + log(exp(log_y1 - c) + exp(log_y2 - c))
        p = exp(log_y2 - log_den)
            # p = g[2]*exp(α[n]+β[n-1]*x[n-1])
            # p /= (p + g[1]*exp(-α[n]-β[n-1]*x[n-1]))
            # print_log("p: $(p)")

        # Sample n-th CpG site
        x[n] = rand() < p ? Int8(1) : Int8(-1)

        # Print for debug
        # print_log("g: $(g)")
        # print_log("p: $(p)")
        # readline()

    end

    # Sample last CpG site
    log_y1 = -α[N] - β[N - 1] * x[N - 1]
    log_y2 = α[N] + β[N - 1] * x[N - 1]
    c = max(log_y1, log_y2)
    log_den = c + log(exp(log_y1 - c) + exp(log_y2 - c))
    p = exp(log_y2 - log_den)
    x[N] = rand() < p ? Int8(1) : Int8(-1)
        # p = exp(α[N]+β[N-1]*x[N-1])
        # p /= (p + exp(-α[N]-β[N-1]*x[N-1]))
        # x[N] = rand()<p ? Int8(1) : Int8(-1)

    # Return x in final state
    return x

end
##################################################################################################
## Permutation testing simulations
##################################################################################################
"""

    `gen_grp_comp_data(S1,S2,ϕ1,ϕ2,SEED)`

    Simulates data for group comparison.

    # Examples
    ```julia-repl
    julia> CpelNano.gen_grp_comp_data(s1,s2,ϕ1,ϕ2,seed)
    ```
"""
function gen_grp_comp_data(s1::Int64, s2::Int64, ϕ1::Vector{Float64}, ϕ2::Vector{Float64}, nse_sc::Float64, seed::Int64)::Tuple{Vector{RegStruct},Vector{RegStruct}}

    # Configuration
    config = CpelNano.CpelNanoConfig()

    # Random genomic properties
    Random.seed!(seed)
    chrst = 1
    chrend = 3000
    cpg_pos = collect(25:25:2975)
    dn = fill(25.0, length(cpg_pos) - 1)
    ρn = rand(0.001:0.001:0.1, length(cpg_pos))

    ## Region genomic properties

    # Group 2
    mods_g1 = Vector{RegStruct}()
    @inbounds for s in 1:s1

        # Parameter vectors
        ϕ1_s = ϕ1 .+ nse_sc .* [rand(-0.1:0.01:0.1),rand(-0.5:0.01:0.5),rand(0.0:0.01:0.2)]

        # Init structure
        rs1 = RegStruct();
        rs1.cpg_pos = cpg_pos; rs1.N = length(cpg_pos); rs1.chrst = chrst;
        rs1.chrend = chrend; rs1.ϕhat = ϕ1_s; rs1.ρn = ρn; rs1.dn = dn;

        # Re-scale models
        rscle_grp_mod!(rs1)

        # Get statistical summaries
        get_stat_sums!(rs1, config)

        # Push
        push!(mods_g1, rs1)

    end

    # Group 2
    mods_g2 = Vector{RegStruct}()
    @inbounds for s in 1:s2

        # Parameter vectors
        ϕ2_s = ϕ2 .+ nse_sc .* [rand(-0.1:0.01:0.1),rand(-0.5:0.01:0.5),rand(0.0:0.01:0.2)]

        # Init structure
        rs2 = RegStruct();
        rs2.cpg_pos = cpg_pos; rs2.N = length(cpg_pos); rs2.chrst = chrst;
        rs2.chrend = chrend; rs2.ϕhat = ϕ2_s; rs2.ρn = ρn; rs2.dn = dn;

        # Re-scale model
        rscle_grp_mod!(rs2)

        # Get statistical summaries
        get_stat_sums!(rs2, config)

        # Push
        push!(mods_g2, rs2)

    end

    # Return tuple of arrays
    return mods_g1, mods_g2

end
###################################################################################################################
# Simulate 1. X̄~q(x̄;ϕ) 2. Y|X̄ ~ 𝒩(μ(X̄),σ(X̄))
###################################################################################################################
"""
    `cpel_samp_ont(m,Nl,αs,βs,pobs,μ_m,σ_m,μ_u,σ_u,𝒜s,ℬs)`

    Simulate entire model where observable is distributed according to a Gaussian distribution. The same
    Gaussian distributions are assumed regardless of the k-mer. The mean is given in terms of mA. Here we
    assume that we are passed

        α_l = (a + b⋅ρ̄_l)⋅N_l ,β_l = c/d_l

    where ρ̄_l and N_l is the average density and the number of CpG sites in the l-th group,respectively.

    # Examples
    ```julia-repl
    julia> L = 10; M = 25; αs = fill(0.0,L); βs = fill(0.0,L-1);
    julia> pobs = 1.0; μ_m = 60.0; σ_m = 1.0; μ_u = 70.0; σ_u = 1.0;
    julia> data = CpelNano.cpel_samp_ont(M,αs,βs,pobs,μ_m,σ_m,μ_u,σ_u);
    ```
"""
function cpel_samp_ont(m::Int64,α::Vector{Float64},β::Vector{Float64},pobs::Float64,μ_m::Float64,σ_m::Float64,
    μ_u::Float64,σ_u::Float64)::RegStruct

    # Distributions p(y|x̄)
    p_y_x_m = Normal(μ_m, σ_m)
    p_y_x_u = Normal(μ_u, σ_u)

    # Generate calls
    L = length(α)
    calls = Vector{Vector{MethCallCpgGrp}}()
    while length(calls) < m
        # Init call vector
        call = [MethCallCpgGrp() for l = 1:L]
        # Get x̄ realization
        xcpel = gen_x_mc(α, β)
        # Fill call vector
        for l in 1:L
            # Check if group is observed
            rand() < pobs || continue
            # Methylation call (convert to bool)
            xbar = xcpel[l] == 1
            # Sample normal observable
            y = xbar ? rand(p_y_x_m) : rand(p_y_x_u)
            # Compute log p(y_l|x̄_l)
            log_pyx_u = log(pdf(p_y_x_u, y))
            log_pyx_m = log(pdf(p_y_x_m, y))
            # Set call for x̄_l
            call[l].obs = true
            call[l].log_pyx_u = log_pyx_u
            call[l].log_pyx_m = log_pyx_m
        end
        push!(calls, call)
    end

    # Create simulation structure
    sim_struct = CpelNano.RegStruct(calls)

    # Return data object
    return sim_struct

end
"""
    `cpel_samp_wgbs(m,Nl,αs,βs,pobs,μ_m,σ_m,μ_u,σ_u,𝒜s,ℬs)`

    Simulate entire model where observable is distributed according to a Gaussian distribution. The same
    Gaussian distributions are assumed regardless of the k-mer. The mean is given in terms of mA. Here we
    assume that we are passed

        α_l = (a + b⋅ρ̄_l)⋅N_l ,β_l = c/d_l

    where ρ̄_l and N_l is the average density and the number of CpG sites in the l-th group,respectively.

    # Examples
    ```julia-repl
    julia> L = 20; M = 20; αs = fill(0.0,L); βs = fill(0.0,L-1); pobs = 0.5;
    julia> data = CpelNano.cpel_samp_wgbs(M,αs,βs,pobs);
    ```
"""
function cpel_samp_wgbs(m::Int64, α::Vector{Float64}, β::Vector{Float64}, pobs::Float64)::RegStruct

    # Generate calls
    L = length(α)
    calls = Vector{Vector{MethCallCpgGrp}}()
    while length(calls) < m
        # Init call vector
        call = [MethCallCpgGrp() for l = 1:L]
        # Init obs counter
        obs_count = 15
        # Get x̄ realization
        xcpel = gen_x_mc(α, β)
        # Fill call vector
        for l in 1:L
            # Check if group is observed
            if rand() < pobs && obs_count >= 15
                obs_count = 0
            end
            # Increase obs counter
            obs_count += 1
            obs_count < 10 || continue
            # Methylation call (convert to bool)
            xbar = xcpel[l] == 1
            # Compute log p(y_l|x̄_l)
            log_pyx_u = xbar ? log_pyx_wrong_x : log_pyx_right_x
            log_pyx_m = xbar ? log_pyx_right_x : log_pyx_wrong_x
            # Set call for x̄_l
            call[l].obs = true
            call[l].log_pyx_u = log_pyx_u
            call[l].log_pyx_m = log_pyx_m
        end
        push!(calls, call)
    end

    # Create simulation structure
    sim_struct = CpelNano.RegStruct(calls)

    # Return data object
    return sim_struct

end
"""
    `bin_calls!(RS)`

    Binarizes calls.

    # Examples
    ```julia-repl
    julia> CpelNano.bin_calls!(rs)
    ```

"""
function bin_calls!(rs::CpelNano.RegStruct)::Nothing

    # Loop over calls
    for call in rs.calls
        # Loop over sites
        for meth_call in call
            # If observed
            if meth_call.obs
                if meth_call.log_pyx_u >= meth_call.log_pyx_m
                    meth_call.log_pyx_u = log_pyx_right_x
                    meth_call.log_pyx_m = log_pyx_wrong_x
                else
                    meth_call.log_pyx_u = log_pyx_wrong_x
                    meth_call.log_pyx_m = log_pyx_right_x
                end
            end
        end
    end

    # Return nothing
    return nothing

end
##################################################################################################
## Metropolis Hastings (NOT USED)
##################################################################################################
"""
    `init_lattice(lat_length)`

     Randomly initializes vector of length `lat_length` with 1 and -1 representing methylated
     and unmethylated CGs respectively.

    # Examples
    ```julia-repl
    julia> CpelNano.init_lattice(2)
    2-element Array{Int8,1}:
     -1
     -1
    ```
"""
function init_lattice(lat_length::Int64)::Vector{Int8}

    # Return initialized lattice
    return rand([Int8(-1),Int8(1)], lat_length)

end
"""
    `get_ΔU(lattice,i,α,β)`

     Computes the change in potential energy U(X) by flipping the methylation state of the `i`th CG
     in a lattice with vectors `[α1,α2,…,αN]` and `[β1,β2,…,β_{N-1}]`.

    # Examples
    ```julia-repl
    julia> CpelNano.get_ΔU([Int8(1),Int8(1)],2,[0.0,0.0],[0.0])
    0.0
    ```
"""
function get_ΔU(lattice::Vector{Int8}, i::Int64, α::Vector{Float64}, β::Vector{Float64})::Float64

    # Modified lattice
    mod_lattice = copy(lattice)
    mod_lattice[i] = -mod_lattice[i]

    # Original contribution of i-th CpG site
    U = -α[i] * lattice[i]
    if i == 1
        U -= β[1] * lattice[i] * lattice[i + 1]
    elseif i == length(lattice)
        U -= β[end] * lattice[i] * lattice[i - 1]
    else
        U -= β[i - 1] * lattice[i - 1] * lattice[i] + β[i] * lattice[i] * lattice[i + 1]
    end

    # Original contribution of i-th CpG site
    Uflip = -α[i] * mod_lattice[i]
    if i == 1
        Uflip -= β[1] * mod_lattice[i] * mod_lattice[i + 1]
    elseif i == length(mod_lattice)
        Uflip -= β[end] * mod_lattice[i] * mod_lattice[i - 1]
    else
        Uflip -= β[i - 1] * mod_lattice[i - 1] * mod_lattice[i] + β[i] * mod_lattice[i] * mod_lattice[i + 1]
    end

    # Note: if ΔU>0 → new config is less likely under θ
# Return ΔU
    return Uflip - U

end
"""
    `metro_hast_inst(lat_length,n_iters,α,β)`

    Run MCMC for `n_iters` iterations on a lattice of length `lat_length` with vectors `[α1,α2,…,αN]`
    and `[β1,β2,…,β_{N-1}]`.

    # Examples
    ```julia-repl
    julia> CpelNano.metro_hast_inst(2,1000,[1.0,1.0],[1.0])
    2-element Array{Int8,1}:
     1
     1
    ```
"""
function metro_hast_inst(lat_length::Int64, n_iters::Int64, α::Vector{Float64}, β::Vector{Float64})::Vector{Int8}

    # Initialize lattice
    lattice = init_lattice(lat_length)

    # Iterate til convergence of MC chain
    @inbounds for i in 1:n_iters
        @inbounds for j in 1:lat_length
            # Get energy Hamiltonian differential
            ΔU = get_ΔU(lattice, j, α, β)
            # If flipping leads to lower energy,flip
            if ΔU < 0
                lattice[j] = -lattice[j]
            else
                # If flipping leads to ΔU>0,flip based on weighted likelihood from Boltzmann distribution
                if rand() < exp(-ΔU)
                    lattice[j] = -lattice[j]
                end
            end
        end
    end

    # Return lattice in final state
    return lattice

end


text:
using Test
using LinearAlgebra
using Eigen

include("powermethod.jl")
include("inverseiteration.jl")
include("rayleighiteration.jl")
include("arnoldi.jl")


text:
using NamedTuples
using Base.Test

@test @NT( a = 1 ).a == 1
@test @NT( a = 1 )[1] == 1
@test @NT( a = 1 )[:a] == 1
@test @NT( a = 12)[:b,2] == 2

@test @NT( :a = 1 ).a == 1
@test @NT( :a = 1 )[1] == 1
@test @NT( :a = 1 )[:a] == 1

@test length( @NT( a = 1)) == 1
@test length( @NT( a = 1, b = 2.0)) == 2

@test values(@NT(a=1,b=2)) == [1,2]
@test get(@NT(a=1,b=2),:a,3) == 1
@test get(@NT(a=1,b=2),:c,3) == 3

@test first( @NT( a = 1, b = 2.0 )) == 1
@test last( @NT( a = 1, b = "hello", c = 2.0 )) == 2.0
@test [ v for v in @NT( a = 1.0, b = 2.0 ) ] == [ 1.0, 2.0 ]

@test ( x = @NT( a::Int64, b::Float64 )( 1, 2.0 ) ; typeof(x.a) == Int64 && typeof(x.b) == Float64 )
@test @NT( a = 1, b = "hello")  ==  @NT( a, b )( 1, "hello")
@test @NT( a = 1) != @NT( b = 1 )
@test @NT(a=1) == @NT(a=1.)

@test hash( @NT( a = 1, b = "hello"))  ==  hash( @NT( a, b )( 1, "hello") )
@test hash( @NT( a = 1, b = "hello")) != hash( @NT( a = 1, b = 2.0 ))

@test @NT( a ) ==  @NT( a )
@test @NT( a ) !=  @NT( b )

@test @NT(a::Int, b) == @NT(a::Int, b::Any)

@test typeof( @NT( a::Int, b::Float64 )(1, 3.0) ) == typeof( @NT( a = 1, b = 2.0 ))

# Syntax tests, including anon named tuples
@test @NT( a, b ) <: NamedTuple
@test @NT( ::Int64, ::Float64 ) <: NamedTuple
@test typeof( @NT( 1, 2, "3" )) <: NamedTuple
@test typeof( @NT( 1 + 2, "hello")) <: NamedTuple
@test length( @NT( 1, 2, 3, "hello")[ 1:2 ]) == 2

@test isbits( @NT( ::Int64, ::Float64)) == true

nt = @NT( a=1, b=2, c=3 )
@test nt.a == 1
@test nt.b == 2
@test nt.c == 3
@test haskey( nt, :x ) == false
x = setindex( nt, :x, 123 )
@test x.x == 123
@test x.a == 1
@test x.b == 2
@test x.c == 3
@test nt[[:c,:b]] == @NT( c=3, b=2 )

y = delete( x, :a)
@test x != y
@test haskey( nt, :a ) == true
@test haskey( y, :a ) == false
@test x.x == 123
@test x.a == 1
@test x.b == 2
@test x.c == 3

@test map(-, @NT(x=1, y=2)) == @NT(x=-1, y=-2)
@test map(+, @NT(x=1, y=2), @NT(x=1, y=2)) == @NT(x=2, y=4)
@test_throws ArgumentError map(+, @NT(x=1, y=2), @NT(y=1, x=2))
@test map(string, @NT(x=1, y=2)) == @NT(x="1", y="2")
@test map(round, @NT(x=1//3, y=Int), @NT(x=3, y=2//3)) == @NT(x=0.333, y=1)

@test @NT(x=1, y=2) <  @NT(x=1, y=2.5)
@test @NT(x=1, y=2) >= @NT(x=1, y=2.0)

@test merge( nt, @NT( d = "hello", e = "world"))  == @NT( a=1,b=2,c=3,d="hello",e="world")

@test get.(@NT( a = Nullable(3), b = Nullable("world") )) == @NT( a = 3, b = "world")
@test_throws MethodError @NT( a = 3) .+ [4]
@test_throws MethodError [4] .+ @NT( a = 3)

# serialize and deserialize
addprocs(1)
@everywhere using NamedTuples
x = @NT(a=1, b=2)
y = @fetchfrom 2 identity(x)
@test isa(y,NamedTuple)
@test y.a == 1
@test y.b == 2

io = IOBuffer()
serialize(io, Union{})
@test deserialize(seekstart(io)) === Union{}

# allow custom types
struct Empty end
nt = @NT(a::Empty, b::Int)
@test nt.parameters[1] == Empty

# type reuse
module A
    using NamedTuples
    const NT = @NT(x)
end

module B
    using NamedTuples
    const NT = @NT(x)
end

@test A.NT === B.NT


text:
@doc """
	had!(A, B) -> A

In-place Hadamard product: Replace `A` with `A.*B`.
"""->
function had!{T<:Number}(A::StridedArray{T}, B::AbstractArray{T})
	size(A) == size(B) || throw(DimensionMismatch())
	for idx in eachindex(A,B)
		@inbounds A[idx] *= B[idx]
	end
end

@doc """
	had!(A, B, C) -> A

In-place Hadamard product: Replace `A` with `B.*C`.
"""->
function had!{T<:Number}(A::StridedArray{T}, B::AbstractArray{T}, C::AbstractArray{T})
	size(A) == size(B) == size(C) || throw(DimensionMismatch())
	for idx in eachindex(A,B)
		@inbounds A[idx] = B[idx] * C[idx]
	end
end

@doc """
	hadc!(A, B) -> A

In-place Hadamard product with complex conjugation: Replace `A` with
`A.*conj(B)`.
"""->
function hadc!{T<:Number}(A::StridedArray{T}, B::AbstractArray{T})
	size(A) == size(B) || throw(DimensionMismatch())
	for idx in eachindex(A,B)
		@inbounds A[idx] *= conj(B[idx])
	end
end

@doc """
	hadc!(A, B, C) -> A

In-place Hadamard product with complex conjugation: Replace `A` with
`B.*conj(C)`.
"""->
function hadc!{T<:Number}(A::StridedArray{T}, B::AbstractArray{T}, C::AbstractArray{T})
	size(A) == size(B) == size(C) || throw(DimensionMismatch())
	for idx in eachindex(A,B)
		@inbounds A[idx] = B[idx] * conj(C[idx])
	end
end

@doc """
	yphad!(y, a, b) -> y

Replace `y` with `y + a.*b`.
"""->
function yphad!{T<:Number}(y::StridedVector{T}, a::AbstractVector{T}, b::AbstractVector{T})
	(Ny = length(y)) == length(a) == length(b) || throw(DimensionMismatch())
	for ny in 1:Ny
		@inbounds y[ny] += a[ny] * b[ny]
	end
end

@doc """
	conj!(A, B) -> A

Replace `A` with `conj(B)`.
"""->
function Base.conj!(A::StridedArray{Complex{Float64}}, B::AbstractArray{Complex{Float64}})
	size(A) == size(B) || throw(DimensionMismatch())
	for idx in eachindex(A,B)
		@inbounds A[idx] = conj(B[idx])
	end
	return A
end

@doc """
	isuniform( points::Matrix ) -> Bool

Test if the rows in the `M`-by-2 matrix `points` are on a uniform 2D grid.
"""->
function isuniform{T<:Real}( points::AbstractMatrix{T} )
	M, D = size(points)
	D == 2 || throw(DimensionMismatch())

	M <= 2 && return true

	x = sort( points[:,1] )
	uniquex = unique(x)
	Mx = length(uniquex)

	y = sort( points[:,2] )
	uniquey = unique(y)
	My = length(uniquey)

	if isuniform(uniquex) && isuniform(uniquey) && Mx*My == M
		return true
	else
		return false
	end
end

@doc """
	grid(M, scale) -> Vector

Equidistant 1D points centered around the origin.
With even `M` the grid has one extra point on the negative values.

The points are scaled by `scale` which by default is 1.
"""->
function grid(M::Int, grid_dist::Real=1.0)
	M >= 2 || throw(AssertionError())
	grid_dist > 0 || throw(DomainError())

	startx = -div(M,2)
	endx = (isodd(M) ? -startx : -startx-1)

	x = collect( float(startx):float(endx) )

	if grid_dist != one(grid_dist)
		scale!(grid_dist, x)
	end

	return x
end

@doc """
	grid( (Mx,My), scale) -> Matrix

2D points on an `Mx`-by-`My` grid centered around the origin.

The points are scaled by `scale` which by default is 1.
"""->
function grid( M::Tuple{Integer,Integer}, grid_dist::Real=1.0)
	minimum(M) >= 2 || throw(AssertionError())
	grid_dist > 0 || throw(DomainError())

	# The points are sorted by the x coordinate
	gridx = grid(M[1], grid_dist)
	x = kron(gridx, ones(M[2]))

	gridy = grid(M[2], grid_dist)
	y = repmat(gridy, M[1])

	return hcat(x, y)
end

@doc """
	weights(xi, bandwidth) -> Vector

Compute weights for sampling points `xi`.

- For 1D points `xi` must be a vector. The function is faster if `xi` is sorted.
- For 2D points `xi` must be a matrix with 2 columns.
"""->
function weights(xi::AbstractVector{Float64}, bandwidth::Real)
	bandwidth > 0 || throw(DomainError())
	maxabs(xi) <= bandwidth || throw(AssertionError())
	(Nx = length(xi)) >= 2 || throw(AssertionError())

	is_xi_sorted = issorted(xi)
	if !is_xi_sorted
		sort_idx = sortperm(xi)
		inv_sort_idx = invperm(sort_idx)
		xi = xi[sort_idx]
	end

	mu = Array{Float64}(Nx)

	# Boundary cases
	L = xi[Nx] - 2*bandwidth
	mu[1] = 0.5*(xi[2] - L)

	U = xi[1] + 2*bandwidth
	mu[Nx] = 0.5*(U - xi[Nx])

	# Non-boundary cases
	for n in 2:Nx-1
		@inbounds mu[n] = 0.5*(xi[n+1] - xi[n-1])
	end

	!is_xi_sorted && permute!(mu, inv_sort_idx)

	return mu
end

function weights(xi::DenseMatrix{Float64}, bandwidth::Real)
	bandwidth > 0 || throw(DomainError())
	size(xi,2) == 2 || throw(DimensionMismatch())
	size(xi,1) >= 2 || throw(DimensionMismatch())
	maxabs(xi) <= bandwidth || throw(AssertionError())

	x = slice(xi, :, 1)
	y = slice(xi, :, 2)
	rw = [-bandwidth; bandwidth; -bandwidth; bandwidth]
	voronoiarea(x, y, rw)
end

@doc """
	density(xi, K) -> δ

Compute the density of the sampling points `xi` with bandwidth `K`:

- In 1D, `xi` must be a vector and δ is the largest difference between two numerically consecutive points in `xi`.
- In 2D, `xi` must be an `M-by-2` matrix and the density is the smallest number δ that allows a covering of the "bandwidth area" with circles centered in the points of `xi` and radii δ.
Currently, the "bandwidth area" is a square centered at the origin and with sidelength `2*K`.

To ensure numerical stability in computations with the associated change of basis matrix, `xi` must contain both negative and positive elements.
"""->
function density(xi::AbstractVector{Float64}, bandwidth::Real)
	(minx = minimum(xi)) < 0 || throw(AssertionError())
	0 < (maxx = maximum(xi)) <= bandwidth || throw(AssertionError())
	(N = length(xi)) >= 2 || throw(DimensionMismatch())

	# Boundary cases
	lower_boundary = maxx - 2*bandwidth
	upper_boundary = minx + 2*bandwidth
	density = max(minx-lower_boundary, upper_boundary-maxx)

	# Non-boundary cases
	xi_sorted = sort(xi)
	for n in 2:N
		@inbounds diff = xi_sorted[n] - xi_sorted[n-1]
		density = max(density, diff)
	end

	return density
end

function density(xi::AbstractMatrix{Float64}, bandwidth::Real)
	bandwidth > 0 || throw(DomainError())
	maxabs(xi) <= bandwidth || throw(AssertionError())
	M, dim = size(xi)
	M >= 2 || throw(DimensionMismatch())
	dim == 2 || throw(DimensionMismatch())

	x = slice(xi, :, 1)
	y = slice(xi, :, 2)
	rw = [-bandwidth; bandwidth; -bandwidth; bandwidth]
	density(x, y, rw)
end

@doc """
	frac(x)
	frac!(x)

The fractional part of `x` as a number in [-0.5, 0.5).
"""->
function frac(x::DenseArray{Float64})
	y = copy(x)
	frac!(y)

	return y
end

function frac!(x::DenseArray{Float64})
	for idx in eachindex(x)
		@inbounds x[idx] -= round( x[idx] )
	end
end

@doc """
	isdaubechies(wavename::AbstractString) -> Bool

Return `true` if `wavename` is of the form `dbN`, where `N` is an integer.
"""->
function isdaubechies(wavename::AbstractString)
	ishaar(wavename) && return true

	prefix = lowercase(wavename[1:2])
	return prefix == "db" && isnumber(wavename[3:end])
end

@doc """
	ishaar(wavename::AbstractString) -> Bool

Return `true` if `wavename` is "haar" or "db1".
"""->
function ishaar(wavename::AbstractString)
	lowername = lowercase(wavename)

	if lowername == "haar" || lowername == "db1"
		return true
	else
		return false
	end
end

@doc """
	hasboundary(wavename::String) -> Bool

Does the `wavename` scaling function have boundary correction.
"""->
function hasboundary(wavename::AbstractString)
	lowername = lowercase(wavename)

	ishaar(lowername) && return false

	if isdaubechies(lowername)
		return true
	else
		error("Only Daubechies scaling functions are valid")
	end
end

@doc """
	split(x::Vector, border) -> L, I, R

Split `x` into 3 slices:
Left, internal and right, where left and right are the `border` outmost entries.
"""->
function Base.split(x::DenseVector, border::Integer)
	border >= 1 || throw(DomainError())
	(N = length(x)) > 2*border || throw(AssertionError())

	L = slice(x, 1:border)
	I = slice(x, border+1:N-border)
	R = slice(x, N-border+1:N)

	return L, I, R
end

type SplitMatrix{T, A<:AbstractMatrix}
	LL::A
	IL::A
	RL::A
	LI::A
	II::A
	RI::A
	LR::A
	IR::A
	RR::A

	SplitMatrix(LL::AbstractMatrix{T}, IL::AbstractMatrix{T},
	RL::AbstractMatrix{T}, LI::AbstractMatrix{T}, II::AbstractMatrix{T},
	RI::AbstractMatrix{T}, LR::AbstractMatrix{T}, IR::AbstractMatrix{T},
	RR::AbstractMatrix{T}) = new(LL, IL, RL, LI, II, RI, LR, IR, RR)
end
SplitMatrix(LL, IL, RL, LI, II, RI, LR, IR, RR) =
SplitMatrix{eltype(LL), typeof(LL)}(LL, IL, RL, LI, II, RI, LR, IR, RR)

@doc """
	split(A::Matrix, border) -> SplitMatrix

Split `A` into 9 parts, where the outer parts are `border` wide/high.
`LL` is e.g. `border`-by-`border`.

	 ________________
	|    |      |    |
	| LL |  LI  | LR |
	|____|______|____|
	|    |      |    |
	| IL |  II  | IR |
	|____|______|____|
	|    |      |    |
	| RL |  RI  | RR |
	|____|______|____|
"""->
function Base.split{T}(A::DenseMatrix{T}, border::Integer)
	border >= 2 || throw(DomainError())
	N = size(A)
	minimum(N) > 2*border || throw(AssertionError())

	Lidx = 1:border
	I1idx = border+1:N[1]-border
	R1idx = N[1]-border+1:N[1]
	I2idx = border+1:N[2]-border
	R2idx = N[2]-border+1:N[2]

	LL = slice(A, Lidx, Lidx)
	IL = slice(A, I1idx, Lidx)
	RL = slice(A, R1idx, Lidx)

	LI = slice(A, Lidx, I2idx)
	II = slice(A, I1idx, I2idx)
	RI = slice(A, R1idx, I2idx)

	LR = slice(A, Lidx, R2idx)
	IR = slice(A, I1idx, R2idx)
	RR = slice(A, R1idx, R2idx)

	SplitMatrix{T, typeof(LL)}( LL, IL, RL, LI, II, RI, LR, IR, RR )
end


text:
struct Grads
    g::Dict
    ps::Dict
end

function Base.getindex(G::Grads, t::T) where T <:AbstractArray
    return G.g[t]
end

c_t(x::Tuple) = x
c_t(x) = Tuple(x)

function grad(m, d, loss, ps)
    layers = m.net
    record = m(d[1], rf = true)
    grads = Dict()
    l = length(layers)
    r = record[l+1]
    Δ = loss(r[2], d[2], back = true)
    σ = layers[l].activation
    Δ = Δ .* σ.(r[1])
    for i in l : -1 : 1
        r = record[i]
        if typeof(layers[i]) <: Conv
            C = layers[i]
            Δ = c_dims(layers[i], Δ)
            _r = Im2Col(r[2], size(C.W)[1:end-2], C.stride, C.padding).x
            grads[ps.ps[i*2-1]] = Δ * _r'
            grads[ps.ps[i*2]] = dropdims(sum(Δ, dims = 2), dims = 2)
        elseif typeof(layers[i]) <: Param
            grads[ps.ps[i*2-1]] = Δ * r[2]'
            grads[ps.ps[i*2]] = Δ
        end
        i == 1 && continue
        Δ = layers[i](Δ, r[1], true)
    end
    Grads(grads, ps.ps)
end

function c_dims(C::Conv, x)
    oc = size(C.W)[end]
    other = Int(length(x) / oc)
    return reshape(x, oc, other)
end


text:
using EzXML
import EzXML: Document, Node

export findalllocal, findfirstlocal, findfirstcontent, findallcontent, addelementOne!, addelementVect!, updateallcontent!, updatefirstcontent!, docOrElmInit, UN
################################################################
UN{T}= Union{T, Nothing}
################################################################
# Extractors

"""
    findfirstlocal(string, node)

findfirst with ignoring namespaces. It considers element.name for returning the elements
"""
function findfirstlocal(s::String, node::Union{Node, Document})
    out = nothing # return nothing if nothing is found
    for child in eachelement(node)
        if child.name == s
            out = child
            break
        end
    end
    return out
end

"""
    findfirstlocal(string, node)

findalllocal with ignoring namespaces. It considers element.name for returning the elements
"""
function findalllocal(s::String, node::Union{Node, Document})
    out = Node[]
    for child in eachelement(node)
        if child.name == s
            push!(out, child)
        end
    end
    if !isempty(out)
        return out
    else # return nothing if nothing is found
        return nothing
    end
end
################################################################
# Single extraction
"""
    findfirstcontent(element,node, amlType)
    findfirstcontent(type,element,node, amlType)

Returns first element content. It also convert to the desired format by passing type. element is given as string.
```julia
findfirstcontent("/instrument-name",node, 0)
findfirstcontent(UInt8,"/midi-channel",node, 0)
```
"""
function findfirstcontent(::Type{T}, s::String, node::Union{Node, Document}, amlType::Int64) where{T<:String} # for strings

    if amlType == 0 # normal elements

        if typeof(node) == Document || hasdocument(node)
            elm = findfirst(s,node)
        else
            elm = findfirstlocal(s,node)
        end

        if isnothing(elm) # return nothing if nothing is found
            return nothing
        else
            return elm.content
        end

    elseif amlType == 2 # Attributes

        if haskey(node, s)
            elm = node[s]
            return elm

        else # return nothing if nothing is found
            elm = nothing
            return elm
        end

    end

end

# if no type is provided consider it to be string
findfirstcontent(s::String,node::Union{Node, Document}, amlType::Int64) = findfirstcontent(Union{String, Nothing}, s, node, amlType)

# for numbers
function findfirstcontent(::Type{T},s::String,node::Union{Node, Document}, amlType::Int64) where {T<:Union{Number,Bool}}

    if amlType == 0 # normal elements

        if typeof(node) == Document || hasdocument(node)
            elm = findfirst(s,node)
        else
            elm = findfirstlocal(s,node)
        end

        if isnothing(elm) # return nothing if nothing is found
            return nothing
        else
            return parse(T, elm.content)
        end

    elseif amlType == 2 # Attributes

        if haskey(node, s)
            elm = parse(T, node[s])
            return elm

        else # return nothing if nothing is found
            elm = nothing
            return elm
        end

    end

end

# for defined types
function findfirstcontent(::Type{T},s::String,node::Union{Node, Document}, amlType::Int64) where {T}

    if amlType == 0 # normal elements

        if typeof(node) == Document || hasdocument(node)
            elm = findfirst(s,node)
        else
            elm = findfirstlocal(s,node)
        end

        if isnothing(elm) # return nothing if nothing is found
            return nothing
        else
            return T(elm)
        end

    elseif amlType == 2 # Attributes

        if haskey(node, s)
            elm = node[s]
            return elm

        else # return nothing if nothing is found
            elm = nothing
            return elm
        end

    end

end

# Union with Nothing
findfirstcontent(::Type{UN{T}},s::String,node::Union{Node, Document}, amlType::Int64) where {T} = findfirstcontent(T,s,node, amlType)

# Nothing Alone
findfirstcontent(::Type{Nothing},s::String,node::Union{Node, Document}, amlType::Int64) = nothing
################################################################
# Vector extraction
"""
    findallcontent(type, string, node, amlType)

Finds all the elements with the address of string in the node, and converts the elements to Type object.
```julia
findallcontent(UInt8,"/midi-channel",node, 0)
```
"""
function findallcontent(::Type{Vector{T}}, s::String, node::Union{Node, Document}, amlType::Int64) where{T<:String} # for strings

    if amlType == 0 # normal elements

        if typeof(node) == Document || hasdocument(node)
            elmsNode = findall(s, node) # a vector of Node elements
        else
            elmsNode = findalllocal(s, node) # a vector of Node elements
        end

        if isnothing(elmsNode)  # return nothing if nothing is found
            return nothing
        else
            elmsType = Vector{T}(undef, length(elmsNode)) # a vector of Type elements
            i=1
            for elm in elmsNode
                elmsType[i]=elm.content
                i+=1
            end
            return elmsType
        end

    elseif amlType == 2 # Attributes

        if haskey(node, s)
            elmsNode = node[s]
            elmsType = Vector{T}(undef, length(elmsNode)) # a vector of Type elements
            i=1
            for elm in elmsNode
                elmsType[i]=elm
                i+=1
            end
            return elmsType
        else  # return nothing if nothing is found
            elmsNode = nothing
        end
    end
end
# if no type is provided consider it to be string
findallcontent(s::String, node::Union{Node, Document}, amlType::Int64) = findallcontent(Vector{Union{String, Nothing}},s, node, amlType)

# for numbers
function findallcontent(::Type{Vector{T}}, s::String, node::Union{Node, Document}, amlType::Int64) where{T<:Union{Number,Bool}}

    if amlType == 0 # normal elements

        if typeof(node) == Document || hasdocument(node)
            elmsNode = findall(s, node) # a vector of Node elements
        else
            elmsNode = findalllocal(s, node) # a vector of Node elements
        end

        if isnothing(elmsNode) # return nothing if nothing is found
            return nothing
        else
            elmsType = Vector{T}(undef, length(elmsNode)) # a vector of Type elements
            i=1
            for elm in elmsNode
                elmsType[i]=parse(T, elm.content)
                i+=1
            end
            return elmsType
        end

    elseif amlType == 2 # Attributes

        if haskey(node, s)
            elmsNode = parse(T, node[s])
            elmsType = Vector{T}(undef, length(elmsNode)) # a vector of Type elements
            i=1
            for elm in elmsNode
                elmsType[i]=parse(T, elm)
                i+=1
            end
            return elmsType
        else  # return nothing if nothing is found
            elmsNode = nothing
        end

    end

end

# for defined types
function findallcontent(::Type{Vector{T}}, s::String, node::Union{Node, Document}, amlType::Int64) where{T}

    if amlType == 0 # normal elements

        if typeof(node) == Document || hasdocument(node)
            elmsNode = findall(s, node) # a vector of Node elements
        else
            elmsNode = findalllocal(s, node) # a vector of Node elements
        end

    elseif amlType == 2 # Attributes

        if haskey(node, s)
            elmsNode = node[s]
        else  # return nothing if nothing is found
            elmsNode = nothing
        end
    end

    if isnothing(elmsNode) # return nothing if nothing is found
        return nothing
    else
        elmsType = Vector{T}(undef, length(elmsNode)) # a vector of Type elements
        i=1
        for elm in elmsNode
            elmsType[i]=T(elm)
            i+=1
        end
        return elmsType
    end

end

# Union with Nothing
findallcontent(::Type{Vector{UN{T}}},s::String,node::Union{Node, Document}, amlType::Int64) where {T} = findallcontent(Vecotr{T},s,node, amlType)

# Nothing Alone
findallcontent(::Type{Vector{Nothing}},s::String,node::Union{Node, Document}, amlType::Int64) = nothing
################################################################
# Constructors

################################################################
# Document
#  defined or nothing for Documents # add strings and others for documents
"""
    addelementOne!(node, name, value, amlType)

Add one element to a node/document
```
"""
function addelementOne!(aml::Document, name::String, value, amlType::Int64)

    if !isnothing(value) # do nothing if value is nothing

        if hasroot(aml)
            amlNode = root(aml)
            if hasmethod(string, Tuple{T})
                if amlType == 0 # normal elements

                    addelement!(aml, name, string(value))
                elseif amlType == 2 # Attributes

                    link!(aml, AttributeNode(name, string(value)))

                end
            else
                link!(amlNode,value.aml)
            end
        else
            setroot!(aml, value.aml)
        end

    end

end

# strings
function addelementOne!(aml::Document, name::String, value::String, amlType::Int64)

    if !isnothing(value) # do nothing if value is nothing

        if hasroot(aml)
            amlNode = root(aml)

            if amlType == 0 # normal elements

                addelement!(amlNode, name, value)

            elseif amlType == 2 # Attributes

                link!(amlNode, AttributeNode(name, value))

            end
        else
            error("You cannot insert a string in the document directly. Define a @aml defined field for xd/hd struct")
        end

    end
end

# number
function addelementOne!(aml::Document, name::String, value::T, amlType::Int64) where {T<:Union{Number, Bool}}

    if !isnothing(value) # do nothing if value is nothing

        if hasroot(aml)
            amlNode = root(aml)

            if amlType == 0 # normal elements

                addelement!(amlNode, name, string(value))
            elseif amlType == 2 # Attributes

                link!(amlNode, AttributeNode(name, string(value)))

            end
        else
            error("You cannot insert a number in the document directly. Define a @aml defined field for xd/hd struct")
        end

    end
end
################################################################
# vector of strings
"""
    addelementVect!(node, name, value, amlType)

Add a vector to a node/document
```
"""
function addelementVect!(aml::Document, name::String, value::Vector{String}, amlType::Int64)

    if hasroot(aml)
        amlNode = root(aml)

        if amlType == 0 # normal elements

            for ii = 1:length(value)
                if !isnothing(value[ii]) # do nothing if value is nothing
                    addelement!(amlNode, name, value[ii])
                end
            end

        elseif amlType == 2 # Attributes

            for ii = 1:length(value)
                if !isnothing(value[ii]) # do nothing if value is nothing
                    link!(amlNode, AttributeNode(name, value[ii]))
                end
            end
        end

    else
        error("You cannot insert a vector in the document directly. Define a @aml defined field for xd/hd struct")
    end

end

# vector of numbers
function addelementVect!(aml::Document, name::String, value::Vector{T}, amlType::Int64) where {T<:Union{Number, Bool}}

    if hasroot(aml)
        amlNode = root(aml)

        if amlType == 0 # normal elements

            for ii = 1:length(value)
                if !isnothing(value[ii]) # do nothing if value is nothing
                    addelement!(amlNode, name, string(value[ii]))
                end
            end

        elseif amlType == 2 # Attributes

            for ii = 1:length(value)
                if !isnothing(value[ii]) # do nothing if value is nothing
                    link!(amlNode, AttributeNode(name, string(value[ii])))
                end
            end
        end

    else
        error("You cannot put string in the document directly.Define a @aml defined field for xd/hd struct")
    end

end

#  vector of defined or nothing
function addelementVect!(aml::Document, name::String, value::Vector{T}, amlType::Int64) where {T}
    if hasroot(aml)
        amlNode = root(aml)

        for ii = 1:length(value)
            if !isnothing(value[ii]) # do nothing if value is nothing
                if hasmethod(string, Tuple{T})
                    if amlType == 0 # normal elements

                        addelement!(amlNode, name, string(value[ii]))
                    elseif amlType == 2 # Attributes

                        link!(amlNode, AttributeNode(name, string(value[ii])))

                    end
                else
                    link!(amlNode,value[ii].aml)
                end
            end
        end

    else
        error("You cannot put string in the document directly. Define a @aml defined field for xd/hd struct")
    end

end

################################################################
# Nodes
# strings
function addelementOne!(aml::Node, name::String, value::String, amlType::Int64)

    if !isnothing(value) # do nothing if value is nothing

        if amlType == 0 # normal elements

            addelement!(aml, name, value)

        elseif amlType == 2 # Attributes

            link!(aml, AttributeNode(name, value))

        end
    end
end

# number
function addelementOne!(aml::Node, name::String, value::T, amlType::Int64) where {T<:Union{Number, Bool}}

    if !isnothing(value) # do nothing if value is nothing

        if amlType == 0 # normal elements

            addelement!(aml, name, string(value))
        elseif amlType == 2 # Attributes

            link!(aml, AttributeNode(name, string(value)))

        end
    end
end

#  defined or nothing
function addelementOne!(aml::Node, name::String, value::T, amlType::Int64) where {T}
    if !isnothing(value)
        if hasmethod(string, Tuple{T})
            if amlType == 0 # normal elements

                addelement!(aml, name, string(value))
            elseif amlType == 2 # Attributes

                link!(aml, AttributeNode(name, string(value)))

            end
        else
            link!(aml,value.aml)
        end
    end
end

# vector of strings
function addelementVect!(aml::Node, name::String, value::Vector{String}, amlType::Int64)

    if amlType == 0 # normal elements

        for ii = 1:length(value)
            if !isnothing(value[ii]) # do nothing if value is nothing
                addelement!(aml, name, value[ii])
            end
        end

    elseif amlType == 2 # Attributes

        for ii = 1:length(value)
            if !isnothing(value[ii]) # do nothing if value is nothing
                link!(aml, AttributeNode(name, value[ii]))
            end
        end
    end
end

# vector of numbers
function addelementVect!(aml::Node, name::String, value::Vector{T}, amlType::Int64) where {T<:Union{Number, Bool}}

    if amlType == 0 # normal elements

        for ii = 1:length(value)
            if !isnothing(value[ii]) # do nothing if value is nothing
                addelement!(aml, name, string(value[ii]))
            end
        end

    elseif amlType == 2 # Attributes

        for ii = 1:length(value)
            if !isnothing(value[ii]) # do nothing if value is nothing
                link!(aml, AttributeNode(name, string(value[ii])))
            end
        end
    end
end

#  vector of defined or nothing
function addelementVect!(aml::Node, name::String, value::Vector{T}, amlType::Int64) where {T}
    for ii = 1:length(value)
        if !isnothing(value[ii]) # do nothing if value is nothing
            if hasmethod(string, Tuple{T})
                if amlType == 0 # normal elements

                    addelement!(aml, name, string(value[ii]))
                elseif amlType == 2 # Attributes

                    link!(aml, AttributeNode(name, string(value[ii])))

                end
            else
                link!(aml,value[ii].aml)
            end
        end
    end
end

################################################################
# Updaters
################################################################
# Single Updater
"""
    updatefirstcontent(value, string, node, amlType)

Updates first element content. It also converts any type to string. element is given as string.
"""
function updatefirstcontent!(value::T, s::String, node::Union{Node, Document}, amlType::Int64) where{T<:Union{String, Number, Bool}} # for strings, number and bool

    if amlType == 0 # normal elements

        if typeof(node) == Document || hasdocument(node)
            elm = findfirst(s,node)
        else
            elm = findfirstlocal(s,node)
        end

        if isnothing(elm) # error if nothing is found
            return error("field not found in aml")
        else
            elm.content = value
        end

    elseif amlType == 2 # Attributes

        if haskey(node, s)
            node[s] = value

        else # error if nothing is found
            return error("field not found in aml")
        end

    end

end

# for defined types
function updatefirstcontent!(value::T, s::String,node::Union{Node, Document}, amlType::Int64) where {T}

    if amlType == 0 # normal elements

        if typeof(node) == Document || hasdocument(node)
            elm = findfirst(s,node)
        else
            elm = findfirstlocal(s,node)
        end

        if isnothing(elm) # error if nothing is found
            return error("field not found in aml")
        else
            if hasmethod(string, Tuple{T})
                elm.content = string(value)
            else
                unlink!(elm)
                link!(node, value.aml)
            end
        end

    elseif amlType == 2 # Attributes

        if haskey(node, s)
            elm = node[s]
            unlink!(elm)
            link!(node, value.aml)

        else # error if nothing is found
            return error("field not found in aml")
        end

    end

end

# Nothing Alone
function updatefirstcontent!(value::Nothing, s::String,node::Union{Node, Document}, amlType::Int64)

    if amlType == 0 # normal elements

        if typeof(node) == Document || hasdocument(node)
            elm = findfirst(s,node)
        else
            elm = findfirstlocal(s,node)
        end

        if isnothing(elm) # error if nothing is found
            return error("field not found in aml")
        else
            unlink!(elm)
        end

    elseif amlType == 2 # Attributes

        if haskey(node, s)
            elm = node[s]
            unlink!(elm)

        else # error if nothing is found
            return error("field not found in aml")
        end

    end

end
################################################################
# Vector update
"""
    updateallcontent!(value, string, node, amlType)

Finds all the elements with the address of string in the node, and updates the content
"""
function updateallcontent!(value::Vector{T}, s::String, node::Union{Node, Document}, amlType::Int64) where{T<:Union{String, Number, Bool}} # for stringsm numbers, and bool

    if amlType == 0 # normal elements

        if typeof(node) == Document || hasdocument(node)
            elmsNode = findall(s, node) # a vector of Node elements
        else
            elmsNode = findalllocal(s, node) # a vector of Node elements
        end

        if isnothing(elmsNode) # error if nothing is found
            return error("field not found in aml")
        else
            i = 1
            for elm in elmsNode
                elm = value[i]
                i+=1
            end
        end

    elseif amlType == 2 # Attributes

        if haskey(node, s)
            elmsNode = node[s]
            i = 1
            for elm in elmsNode
                elm = value[i]
                i+=1
            end
        else # error if nothing is found
            return error("field not found in aml")
        end
    end
end

# for defined types and nothing
function updateallcontent!(value::Vector{T}, s::String, node::Union{Node, Document}, amlType::Int64) where{T}

    if amlType == 0 # normal elements

        if typeof(node) == Document || hasdocument(node)
            elmsNode = findall(s, node) # a vector of Node elements
        else
            elmsNode = findalllocal(s, node) # a vector of Node elements
        end

    elseif amlType == 2 # Attributes

        if haskey(node, s)
            elmsNode = node[s]
        else # error if nothing is found
            return error("field not found in aml")
        end
    end

    if isnothing(elmsNode) # error if nothing is found
        return error("field not found in aml")
    else
        i = 1
        for elm in elmsNode
            if isnothing(value[i])
                unlink!(elm)
            else
                if hasmethod(string, Tuple{T})
                    elm.content = string(value[i])
                else
                    unlink!(elm)
                    link!(node, value[i].aml)
                end
            end
            i+=1
        end
        return elmsType
    end

end

################################################################
# doc or element initialize
"""
    docOrElmInit(name)
    docOrElmInit(type, name)

Function to initialize the aml

type:
0 : element node # default
10: empty element node
-1: xml
-2: html
"""
function docOrElmInit(type::Int64 = 0, name::String = nothing)

    if type == 0 # element node

        out = ElementNode(name)

    elseif type == 10 # empty element node

        out = ElementNode(name)

    elseif type == -1 # xml

        out = XMLDocument() # version 1

    elseif type == -2 # html

        out = HTMLDocument() # no URI and external id
    end

    return out
end
################################################################
# moved to my fork of EzXML
# function Base.print(x::Node)
#     println("")
#     prettyprint(x)
# end
# function Base.print(x::Document)
#     println("")
#     prettyprint(x)
# end


text:
#
# This is application file of Satellite. It tracks some process
# and rerun it in case of crash. This is needed, because Julia has
# an issue with long term eval() function running. Github issue:
# https://github.com/JuliaLang/julia/issues/15017
#
# @author DeadbraiN
#
include("global/ImportFolders.jl")
import Config.@if_debug
import Helper
import Backup
import Config
#
# This function starts the satellite and begin watching on
# Manager application (process).
#
function main()
  if length(ARGS) < 1
    Helper.warn("Application file required. e.g.: julia AppSatellite.jl App.jl")
    return false
  end
  local stamp::Float64 = time()
  local cfg::Config.ConfigData = Config.create()

  if Config.isEmpty(cfg) quit() end
  #
  # In case of error or if non zero exit code will be returned
  # an exeption will be thrown.
  #
  while true
    try
      while true
        run(`julia --color=yes $(ARGS[1]) recover $(ARGS[2:end])`)
        break
      end
      break
    catch e
      stamp = _removeBrokenBackup(stamp, cfg)
      @if_debug showerror(STDOUT, e, catch_backtrace())
    end
  end

  true
end
#
# This is a fix for broken backup files. It's related to issue:
# https://github.com/JuliaLang/julia/issues/15017
# If app crashes less then backupPeriod, then we have to remove last
# backup file, because it contains broken organisms code or somthing
# similar. The fact that we are inside this function means that our
# application has crashed a moment ago.
# @param stamp Timestamp of previous crash
# @param cfg Global configuration object
# @return Updated last crash timestamp
#
function _removeBrokenBackup(stamp::Float64, cfg::Config.ConfigData)
  local last::String

  if (time() - stamp) < cfg.backupPeriod
    try
      if (last = Backup.lastFile()) != ""
        Helper.warn("Removing broken backup file: $last")
        rm(Backup.FOLDER_NAME * "/" * last)
      end
    catch e
      Helper.error("Something wrong with backup file removing: $last")
      @if_debug showerror(STDOUT, e, catch_backtrace())
    end
  end

  (stamp = time())
end
#
# Application entry point
#
main()

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo


text:
VERSION >= v"0.4.0-dev+6521" && __precompile__(true)

module Packing
using GeometryTypes

include("rectangle.jl")
export RectanglePacker

end # module


text:
function main()
end

if nameof(@__MODULE__) == :Main
end


text:
using Revise
using Random
using POMDPs
using Spot
using POMDPModelChecking
using POMDPPolicies
using BeliefUpdaters
using SARSOP
using DroneSurveillance
using POMDPSimulators
import Cairo
using POMDPGifs
using ProgressMeter

# pomdp = DroneSurveillancePOMDP(size=(5,5), camera=PerfectCam(), agent_policy=:restricted)

pomdp = DroneSurveillancePOMDP(size=(8,8), camera=PerfectCam())

# pomdp = DroneSurveillancePOMDP(size=(10,10), camera=PerfectCam())

function POMDPModelChecking.labels(pomdp::DroneSurveillancePOMDP, s::DSState, a::Int64)
    return labels(pomdp, s)
end
function POMDPModelChecking.labels(pomdp::DroneSurveillancePOMDP, s::DSState)
    if s.quad == pomdp.region_A
        # return (:a,)
        return ()
    elseif s.quad == pomdp.region_B
        return (:b,)
        # return ()
    elseif s.quad == s.agent && !(isterminal(pomdp, s))
        return (:det,)
        # return ()
    elseif isterminal(pomdp, s)
        # return (:exit,)
        return ()
    else
        return ()
    end
end

@show n_states(pomdp)
@show n_actions(pomdp)
@show n_observations(pomdp)

prop = ltl"!det U b"
# prop = ltl"G!det & G!exit"
# prop = ltl"F b"
# prop = ltl"G !hover & G !det & G!exit"
# prop = ltl"!det U b & !det U a & G!exit"
# prop = ltl"!det U b & !det U a & b => !det U a"

run(`rm model.pomdpx`)
solver = ModelCheckingSolver(property = prop,
                             solver = SARSOPSolver(precision=1e-2), verbose=true)

policy = solve(solver, pomdp);

rng = MersenneTwister(2)
up = DiscreteUpdater(policy.problem)
b0 = initialize_belief(up, initialstate_distribution(policy.problem))
s0 = initialstate(pomdp, rng)
s0 = DSState([1,1], [1,2])
hr = HistoryRecorder(max_steps=20)
hist = simulate(hr, policy.problem, policy, up, b0);
success = undiscounted_reward(hist) > 0.0
state_hist = [s.s for s in hist.state_hist]
shist = POMDPHistory(state_hist, [getfield(hist, f) for f in fieldnames(POMDPHistory) if f != :state_hist]...);
makegif(pomdp, shist, filename="test.gif", spec="(s,a)")

@showprogress for n=1:100
    hr = HistoryRecorder(max_steps=10)
    hist = simulate(hr, policy.problem, policy, up, b0);
    success = undiscounted_reward(hist) > 0.0
    println("discounted reward: ", discounted_reward(hist))
    if !success
        println("Crash")
        state_hist = [s.s for s in hist.state_hist]
        shist = POMDPHistory(state_hist, [getfield(hist, f) for f in fieldnames(POMDPHistory) if f != :state_hist]...);
        makegif(pomdp, shist, filename="test.gif", spec="(s,a)")
    end
end

actionvalues(policy, b0)

undiscounted_reward(hist)

b0 = b

for (i,s) in enumerate(b0.state_list)
    if b0.b[i] != 0
        println(s)
    end
end

# predict
using Printf
a = 1
for (i,s) in enumerate(b0.state_list)
    if b0.b[i] != 0
        d = transition(policy.problem, s, a)
        for ss in d.vals
            sss = (ss.s.quad, ss.s.agent, s.q)
            @printf("%s", ss)
        end
        @printf("\n")
    end
end

s0 = rand(rng, b0)
sp, o, r = generate_sor(policy.problem, s0, 1, rng)

b = update(up, b0, 5, o);

a = actionvalues(policy, b)

empty!(policy.problem.accepting_states)

s0 = initialstate(pomdp, rng)
hr = HistoryRecorder(max_steps=10)
hist = simulate(hr, pomdp, policy, up, b0, s0);
state_hist = [s for s in hist.state_hist]
hist = POMDPHistory(state_hist, [getfield(hist, f) for f in fieldnames(POMDPHistory) if f != :state_hist]...);
makegif(pomdp, hist, filename="test.gif", spec="(s,a)")

function deterministic_belief(pomdp, s)
    b = zeros(length(states(pomdp)))
    si = stateindex(pomdp, s)
    b[si] = 1.0
    return DiscreteBelief(pomdp, b)
end

dra = DeterministicRabinAutomata(prop)
sink_state = ProductState(first(states(pomdp)), -1)
pmdp = ProductPOMDP(pomdp, dra, sink_state)

mecs = maximal_end_components(pmdp, verbose=true)

stateindex(pmdp, ProductState(DSState([5,5], [1,1]), 1))

for si in mecs[2]
    s = states(pmdp)[si]
end


text:
using Reexport

@reexport module DelaunayTriangulations
    include("AbstractDelaunayTriangulation.jl")
    include("composite_types/DelaunayTriangulation.jl")

    include("plot_recipes/plot_recipes.jl")
end

"""
    DelaunayTriangulations

A module handling the creation of Delaunay triangulations from data.
"""
DelaunayTriangulations


text:
using Pkg
Pkg.activate(pwd())

# # Monte Carlo sampling
# ## Gamma function

using SpecialFunctions
using Plots

plot(0:0.1:10, gamma;
    xlabel="x",
    ylabel="gamma(x): log scale",
    label="",
    yscale=:log10,
)

# ## Volume of m-dimensional ball
# ### Exercise:
# Use the [formula](https://en.wikipedia.org/wiki/Volume_of_an_n-ball) to compute the volume
# of a $m$-dimensional ball. Plot the dependence of the volume on the dimension
# $m=1, \dots, 100$.
#
# ---
# ### Solution:

# ---
#
# ### Exercise:
# Write the function `volume_monte_carlo`, which estimates the volume of the
# $m$-dimensional ball based on $n$ randomly sampled points.
#
# **Hint**: function `rand(m,n)` creates a $m\times n$ matrix, which can be understood as
# $n$ randomly sampled points in $[0,1]^m$. Transform them to $[-1,1]^m$.
#
# ---
# ### Solution:

# ---

ms = 1:15
ns = Int64.([1e1; 1e3; 1e5])

Random.seed!(666)

plt = plot(ms, m -> volume_true(m, 1);
    xlabel="dimension",
    ylabel="unit ball volume",
    legend=:topleft,
    label="True",
    line=(4,:black),
)

#+

for n in ns
    plot!(plt, ms, m -> volume_monte_carlo.(m; n=n); label="n = $n")
end

display(plt)

#+

using Distributions

rand(Uniform(-1, 1), 10, 5)

# ## Sampling from distributions

using Distributions

d1 = Normal()
d2 = Normal(1, 1)
d3 = Normal(0, 0.01)

#+

f1(x) = pdf(d1, x)
f2(x) = pdf(d2, x)
f3(x) = pdf(d3, x)

#+

function plot_histogram(xs, f; kwargs...)
    plt = histogram(xs;
        label="Sampled density",
        xlabel = "x",
        ylabel = "pdf(x)",
        nbins = 85,
        normalize = :pdf,
        opacity = 0.5,
        kwargs...
    )

    plot!(plt, range(minimum(xs), maximum(xs); length=100), f;
        label="True density",
        line=(4,:black),
    )

    return plt
end

#+

plot_histogram(rand(d1, 1000000), f1)

# ### Exercise:
# Implement the `rejection_sampling` function. It should generate $n$ trial points and
# return all accepted points.
#
# ---
# ### Solution:

# ---

xlims = (-10, 10)

for (f, d) in zip((f1, f2, f3), (d1, d2, d3))
    Random.seed!(666)
    xs = rejection_sampling(f, f(d.μ), xlims...)

    plt = plot_histogram(xs, f)
    display(plt)
end

#+

h(x) = cos(100*x)

Δx = 0.001
xs = range(xlims...; step=Δx)
e0 = Δx * sum(f3.(xs) .* h.(xs))

#+

expectation1(h, d; n = 1000000) = mean(h, rand(d, n))

function expectation2(h, f, f_max, xlims; n=1000000)
    return mean(h, rejection_sampling(f, f_max, xlims...; n))
end

function expectation3(h, f, d_gen; n=1000000)
    g(x) = h(x)*f(x)/pdf(d_gen, x)
    return mean(g, rand(d_gen, n))
end

#+

n = 100000
n_rep = 20

Random.seed!(666)
e1 = [expectation1(h, d3; n=n) for _ in 1:n_rep]
e2 = [expectation2(h, f3, f3(d3.μ), xlims; n=n) for _ in 1:n_rep]
e3 = [expectation3(h, f3, d1; n=n) for _ in 1:n_rep]

#+

scatter([1], [e0]; label="Integral discretization", legend=:topleft);
scatter!(2*ones(n_rep), e1; label="Generating from Distributions.jl");
scatter!(3*ones(n_rep), e2; label="Generating from rejection sampling");
scatter!(4*ones(n_rep), e3; label="Generating from other distribution")

using Pkg
Pkg.activate(pwd())

# # How many samples do we need?
# ### Exercise:
#
# Sample $n=1000$ in the $m=9$-dimensional space. What is the minimum distance of these
# points? Before implementing the exercise, try to guess the answer.
#
# ---
# ### Solution:

# ---

quantile_sampled1(d, n::Int, α) = partialsort(rand(d, n), floor(Int, α*n))
quantile_sampled2(d, n::Int, α) = quantile(rand(d, n), α)
quantile_sampled1(d, ns::AbstractVector, α) = quantile_sampled1.(d, ns, α)
quantile_sampled2(d, ns::AbstractVector, α) = quantile_sampled2.(d, ns, α)

#+

α = 0.99
n_rep = 20
ns = round.(Int, 10 .^ (1:0.05:5))

Random.seed!(666)
qs1 = hcat([quantile_sampled1(d1, ns, α) for _ in 1:n_rep]...)

#+

Random.seed!(666)
qs2 = hcat([quantile_sampled2(d1, ns, α) for _ in 1:n_rep]...)

#+

plt1 = plot([0.9*minimum(ns); 1.1*maximum(ns)], quantile(d1, α)*ones(2);
    xlabel="n: log scale",
    ylabel="sampled quantile",
    xscale=:log10,
    label="True quantile",
    line=(4,:black),
    ylims=(0,3.5),
);
plt2 = deepcopy(plt1);

for i in 1:size(qs1,1)
    scatter!(plt1, ns[i]*ones(size(qs1,2)), qs1[i,:];
        label="",
        color=:blue,
        markersize = 2,
    )
    scatter!(plt2, ns[i]*ones(size(qs2,2)), qs2[i,:];
        label="",
        color=:blue,
        markersize = 2,
    )
end

plot!(plt1, ns, mean(qs1; dims=2);
    label="Sampled mean",
    line=(4,:red),
);
plot!(plt2, ns, mean(qs2; dims=2);
    label="Sampled mean",
    line=(4,:red),
);

display(plt1)
display(plt2)


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

mutable struct IoK8sApiDiscoveryV1beta1EndpointPort <: SwaggerModel
    appProtocol::Any # spec type: Union{ Nothing, String } # spec name: appProtocol
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    port::Any # spec type: Union{ Nothing, Int32 } # spec name: port
    protocol::Any # spec type: Union{ Nothing, String } # spec name: protocol

    function IoK8sApiDiscoveryV1beta1EndpointPort(;appProtocol=nothing, name=nothing, port=nothing, protocol=nothing)
        o = new()
        validate_property(IoK8sApiDiscoveryV1beta1EndpointPort, Symbol("appProtocol"), appProtocol)
        setfield!(o, Symbol("appProtocol"), appProtocol)
        validate_property(IoK8sApiDiscoveryV1beta1EndpointPort, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(IoK8sApiDiscoveryV1beta1EndpointPort, Symbol("port"), port)
        setfield!(o, Symbol("port"), port)
        validate_property(IoK8sApiDiscoveryV1beta1EndpointPort, Symbol("protocol"), protocol)
        setfield!(o, Symbol("protocol"), protocol)
        o
    end
end # type IoK8sApiDiscoveryV1beta1EndpointPort

const _property_map_IoK8sApiDiscoveryV1beta1EndpointPort = Dict{Symbol,Symbol}(Symbol("appProtocol")=>Symbol("appProtocol"), Symbol("name")=>Symbol("name"), Symbol("port")=>Symbol("port"), Symbol("protocol")=>Symbol("protocol"))
const _property_types_IoK8sApiDiscoveryV1beta1EndpointPort = Dict{Symbol,String}(Symbol("appProtocol")=>"String", Symbol("name")=>"String", Symbol("port")=>"Int32", Symbol("protocol")=>"String")
Base.propertynames(::Type{ IoK8sApiDiscoveryV1beta1EndpointPort }) = collect(keys(_property_map_IoK8sApiDiscoveryV1beta1EndpointPort))
Swagger.property_type(::Type{ IoK8sApiDiscoveryV1beta1EndpointPort }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiDiscoveryV1beta1EndpointPort[name]))}
Swagger.field_name(::Type{ IoK8sApiDiscoveryV1beta1EndpointPort }, property_name::Symbol) =  _property_map_IoK8sApiDiscoveryV1beta1EndpointPort[property_name]

function check_required(o::IoK8sApiDiscoveryV1beta1EndpointPort)
    true
end

function validate_property(::Type{ IoK8sApiDiscoveryV1beta1EndpointPort }, name::Symbol, val)
end


text:
using Test: @test, @testset
using LinearAlgebra: norm
using Statistics: mean
using RestrictedBoltzmannMachines: zerosum, zerosum!, zerosum_weights, extensive_sample, free_energy
using RestrictedBoltzmannMachines: RBM, Potts, Binary, Spin, Gaussian, ReLU, dReLU, pReLU, xReLU

@testset "zerosum (visible Potts)" begin
    q = 3
    N = (2, 3)
    M = (2, 3)
    rbm = RBM(Potts(randn(q, N...)), Binary(randn(M...)), randn(q, N..., M...))
    rbm1 = zerosum(rbm)
    @test norm(mean(rbm1.w; dims=1)) < 1e-13
    @test norm(mean(rbm1.visible.θ; dims=1)) < 1e-13
    v = extensive_sample(rbm.visible)
    F0 = free_energy(rbm, v)
    F1 = free_energy(rbm1, v)
    @test all(F0 - F1 .≈ mean(F0 - F1))
    @test zerosum_weights(rbm.w, rbm) ≈ zerosum(rbm).w
    zerosum!(rbm)
    @test free_energy(rbm, v) ≈ F1
end

@testset "zerosum (hidden Potts)" begin
    p = 3
    N = (2, 3)
    M = (2, 3)
    rbm = RBM(Binary(randn(N...)), Potts(randn(p, M...)), randn(N..., p, M...))
    rbm1 = zerosum(rbm)
    @test norm(mean(rbm1.w; dims=3)) < 1e-13
    @test norm(mean(rbm1.hidden.θ; dims=1)) < 1e-13
    v = extensive_sample(rbm.visible)
    F0 = free_energy(rbm, v)
    F1 = free_energy(rbm1, v)
    @test all(F0 - F1 .≈ mean(F0 - F1))
    @test zerosum_weights(rbm.w, rbm) ≈ zerosum(rbm).w
    zerosum!(rbm)
    @test free_energy(rbm, v) ≈ F1
end

@testset "zerosum (visible and hidden Potts)" begin
    q = 3
    p = 3
    N = (2, 3)
    M = (2, 3)
    rbm = RBM(Potts(randn(q, N...)), Potts(randn(p, M...)), randn(q, N..., p, M...))
    rbm1 = zerosum(rbm)
    @test norm(mean(rbm1.w; dims=1)) < 1e-13
    @test norm(mean(rbm1.w; dims=4)) < 1e-13
    @test norm(mean(rbm1.visible.θ; dims=1)) < 1e-13
    @test norm(mean(rbm1.hidden.θ; dims=1)) < 1e-13
    v = extensive_sample(rbm.visible)
    F0 = free_energy(rbm, v)
    F1 = free_energy(rbm1, v)
    @test all(F0 - F1 .≈ mean(F0 - F1))
    @test zerosum_weights(rbm.w, rbm) ≈ zerosum(rbm).w
    zerosum!(rbm)
    @test free_energy(rbm, v) ≈ F1
end

@testset "zerosum (visible Potts, various hidden layers)" begin
    q = 3
    N = (2, 3)
    M = (2, 3)
    hidden_layers = (
        Binary(M...), Spin(M...), Gaussian(M...), ReLU(M...),
        dReLU(M...), pReLU(M...), xReLU(M...)
    )
    for hidden in hidden_layers
        rbm = RBM(Potts(randn(q, N...)), hidden, randn(q, N..., M...))
        rbm1 = zerosum(rbm)
        @test norm(mean(rbm1.w; dims=1)) < 1e-13
        @test norm(mean(rbm1.visible.θ; dims=1)) < 1e-13
        v = extensive_sample(rbm.visible)
        F0 = free_energy(rbm, v)
        F1 = free_energy(rbm1, v)
        @test all(F0 - F1 .≈ mean(F0 - F1))
        zerosum!(rbm)
        @test free_energy(rbm, v) ≈ F1
    end
end


text:
export
    namfrm,
    ncpos,
    ncposr,
    nearpt,
    npedln,
    npelpt,
    nplnpt,
    nvc2pl,
    nvp2pl

"""
    namfrm(frname)

Look up the frame ID code associated with a string.

### Arguments ###

- `frname`: The name of some reference frame

### Output ###

The SPICE ID code of the frame.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/namfrm_c.html)
"""
function namfrm(frname)
    frcode = Ref{SpiceInt}()
    ccall((:namfrm_c, libcspice), Cvoid, (Cstring, Ref{SpiceInt}), frname, frcode)
    handleerror()
    Int(frcode[])
end

"""
    ncpos(str, chars, start)

Find the first occurrence in a string of a character NOT belonging to a
collection of characters, starting at a specified location, searching forward.

### Arguments ###

- `str`: A string
- `chars`: A collection of characters
- `start`: Position to begin looking for a character not in `chars`

### Output ###

Returns the index of the first character of `str` at or following index `start`
that is not in the collection `chars`.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/ncpos_c.html)
"""
function ncpos(str, chars, start)
    res = ccall((:ncpos_c, libcspice), SpiceInt, (Cstring, Cstring, Cint),
          str, String(chars), start - 1)
    handleerror()
    Int(res) + 1
end

"""
    ncposr(str, chars, start)

Find the first occurrence in a string of a character NOT belonging to a
collection of characters, starting at a specified location, searching
in reverse.

### Arguments ###

- `str`: A string
- `chars`: A collection of characters
- `start`: Position to begin looking for a character not in `chars`

### Output ###

Returns the index of the last character of `str` at or before index `start`
that is not in the collection `chars`.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/ncposr_c.html)
"""
function ncposr(str, chars, start)
    res = ccall((:ncposr_c, libcspice), SpiceInt, (Cstring, Cstring, SpiceInt),
          str, String(chars), start - 1)
    handleerror()
    Int(res) + 1
end

"""
    nearpt(positn, a, b, c)

This routine locates the point on the surface of an ellipsoid that is nearest
to a specified position. It also returns the altitude of the position above
the ellipsoid.

### Arguments ###

- `positn`: Position of a point in the bodyfixed frame
- `a`: Length of semi-axis parallel to x-axis
- `b`: Length of semi-axis parallel to y-axis
- `c`: Length on semi-axis parallel to z-axis

### Output ###

Returns a tuple consisting of `npoint` and `alt`.

- `npoint`: Point on the ellipsoid closest to `positn`
- `alt`: Altitude of `positn` above the ellipsoid

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/nearpt_c.html)
"""
function nearpt(positn, a, b, c)
    @checkdims 3 positn
    npoint = Array{SpiceDouble}(undef, 3)
    alt = Ref{SpiceDouble}()
    ccall((:nearpt_c, libcspice), Cvoid,
          (Ref{SpiceDouble}, SpiceDouble, SpiceDouble, SpiceDouble,
           Ref{SpiceDouble}, Ref{SpiceDouble}),
          collect(positn), a, b, c, npoint, alt)
    handleerror()
    npoint, alt[]
end

"""
    npedln(a, b, c, linept, linedr)

Find nearest point on a triaxial ellipsoid to a specified line, and the
distance from the ellipsoid to the line.

### Arguments ###

- `a`: Length of semi-axis in the x direction
- `b`: Length of semi-axis in the y direction
- `c`: Length of semi-axis in the z direction
- `linept`: Point on line
- `linedr`: Direction vector of line

### Output ###

Returns a tuple consisting of `pnear` and `dist`.

- `pnear`: Nearest point on ellipsoid to line
- `dist`: Distance of ellipsoid from line

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/npedln_c.html)
"""
function npedln(a, b, c, linept, linedr)
    @checkdims 3 linept linedr
    pnear = Array{SpiceDouble}(undef, 3)
    dist = Ref{SpiceDouble}()
    ccall((:npedln_c, libcspice), Cvoid,
          (SpiceDouble, SpiceDouble, SpiceDouble, Ref{SpiceDouble},
           Ref{SpiceDouble}, Ref{SpiceDouble}, Ref{SpiceDouble}),
          a, b, c, collect(linept), collect(linedr), pnear, dist)
    handleerror()
    pnear, dist[]
end

"""
    npelpt(point, ellips)

Find the nearest point on an ellipse to a specified point, both in
three-dimensional space, and find the distance between the ellipse
and the point.

### Arguments ###

- `point`: Point whose distance to an ellipse is to be found
- `ellips`: A SPICE ellipse

### Output ###

Returns a tuple consisting of `pnear` and `dist`.

- `pnear`: Nearest point on ellipse to input point
- `dist`: Distance of input point to ellipse

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/npelpt_c.html)
"""
function npelpt(point, ellips)
    @checkdims 3 point
    pnear = Array{SpiceDouble}(undef, 3)
    dist = Ref{SpiceDouble}()
    ccall((:npelpt_c, libcspice), Cvoid,
          (Ref{SpiceDouble}, Ref{Ellipse}, Ref{SpiceDouble}, Ref{SpiceDouble}),
          collect(point), ellips, pnear, dist)
    handleerror()
    pnear, dist[]
end

"""
    nplnpt(linept, linedr, point)

Find the nearest point on a line to a specified point, and find the distance
between the two points.

### Arguments ###

- `linept`: Point on line
- `linedr`: Direction vector of line
- `point`: A second point

### Output ###

Returns a tuple consisting of `pnear` and `dist`.

- `pnear`: Nearest point on the line to `point`
- `dist`: Distance between `point` and `pnear`

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/nplnpt_c.html)
"""
function nplnpt(linept, linedr, point)
    @checkdims 3 linept linedr point
    pnear = Array{SpiceDouble}(undef, 3)
    dist = Ref{SpiceDouble}()
    ccall((:nplnpt_c, libcspice), Cvoid,
          (Ref{SpiceDouble}, Ref{SpiceDouble}, Ref{SpiceDouble},
           Ref{SpiceDouble}, Ref{SpiceDouble}),
          collect(linept), collect(linedr), collect(point), pnear, dist)
    handleerror()
    pnear, dist[]
end

"""
    nvc2pl(norm, point)

Make a SPICE plane from a normal vector and a point.

### Arguments ###

- `norm`: A normal vector...
- `constant`: ...and a constant defining a plane

### Output ###

Returns a struct representing the plane.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/nvc2pl_c.html)
"""
function nvc2pl(norm, constant)
    @checkdims 3 norm
    plane = Ref{Plane}()
    ccall((:nvc2pl_c, libcspice), Cvoid,
          (Ref{SpiceDouble}, SpiceDouble, Ref{Plane}), collect(norm), constant, plane)
    handleerror()
    plane[]
end

"""
    nvp2pl(norm, point)

Make a SPICE plane from a normal vector and a point.

### Arguments ###

- `norm`: A normal vector...
- `point`: ...and a point defining a plane

### Output ###

Returns a struct representing the plane.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/nvp2pl_c.html)
"""
function nvp2pl(norm, orig)
    @checkdims 3 norm orig
    plane = Ref{Plane}()
    ccall((:nvp2pl_c, libcspice), Cvoid,
          (Ref{SpiceDouble}, Ref{SpiceDouble}, Ref{Plane}),
          collect(norm), collect(orig), plane)
    handleerror()
    plane[]
end