text:
using DocStringExtensions
import Base
using DataStructures

"""
A node represents either a number or operator.
  If number, then the node is a leaf.
  If operator, then it has left and right nodes.
A mathematical expression can be represented by a binary tree of the nodes.
E.g.,
(1 + 2) * (7 - 2) is represented by

*
├─ +
│  ├─ 1
│  └─ 2
└─ -
   ├─ 7
   └─ 2
"""
struct Node
    value::Rational
    str::String
    isnumber::Bool
    left::Union{Node,Nothing}
    right::Union{Node,Nothing}

    Node(x::Rational) = new(x, replace(string(x), r"//" => "/"),
                            true, nothing, nothing)
    Node(x::Int) = new(Rational(x), string(x), true, nothing, nothing)
    Node(op, left::Node, right::Node) = new(op(left.value, right.value),
                                            op == (//) ? "/" : string(op),
                                            false, left, right)
end

Base.:+(a::Node, b::Node) = Node(+, a, b)
Base.:-(a::Node, b::Node) = Node(-, a, b)
Base.:*(a::Node, b::Node) = Node(*, a, b)
Base.:/(a::Node, b::Node) = Node(//, a, b)

"""
Convert a tree into a mathematical expression
"""
function toexpr(node::Node, opprev::String, isright::Bool)::String
    if node.isnumber
        if occursin("/", node.str)
            # rational number should be treated as a chunk
            return "($(node.str))"
        else
            return node.str
        end
    end

    paren = needparen(opprev, node.str, isright)
    ret = ""
    paren && (ret *= "(")
    ret *= toexpr(node.left, node.str, false)
    ret *= " $(node.str) "
    ret *= toexpr(node.right, node.str, true)
    paren && (ret *= ")")
    ret
end

toexpr(node::Node) = toexpr(node, "", false)

needparen(opprev::String, op::String, isright::Bool) =
    (op=="/" && opprev in ("*","/") && isright) ||
    (op in ("+", "-") && opprev in ("*","/")) ||
    (op in ("+", "-") && opprev=="-" && isright)

#===================================================================================#

"""
Reduction iterator, given a set of nodes, generates a set of nodes reduced by one-step computation
"""
struct Reductions
    ctr::Accumulator{Node,Int}
    nodes::Vector{Node}
    counts::Vector{Int}

    function Reductions(ctr::Accumulator{Node,Int})
        nodes = [k for k in keys(ctr)]
        counts = [v for v in values(ctr)]
        new(ctr, nodes, counts)
    end
    Reductions(xs::Node...) = Reductions(counter(xs))
end

function Base.iterate(r::Reductions, state::Int)
    n = length(r.nodes)  # number of unique nodes
    (state >= 4*n*n) && return nothing

    k = div(state, n*n)
    i = div(state-n*n*k, n)
    j = state-n*n*k-n*i
    i += 1; j += 1; k += 1

    n1, n2 = r.nodes[i], r.nodes[j]
    op = (+,-,*,/)[k]

    # invalid cases
    # same node, but count is smaller than 2
    (i==j) && (r.counts[i] < 2) && return iterate(r, state+1)
    # + and * are order independent
    (i>j) && (op in (+,*)) && return iterate(r, state+1)
    # zero division
    (n2.value==Rational(0)) && (op == (/)) && return iterate(r, state+1)

    remain = copy(r.ctr)
    n = op(n1, n2)
    inc!(remain, n)
    dec!(remain, n1)==0 && delete!(remain.map, n1)
    dec!(remain, n2)==0 && delete!(remain.map, n2)
    remain, state+1
end

Base.iterate(r::Reductions) = iterate(r, 0)

#===================================================================================#

"""
Extract only numbers from a set of nodes
"""
valueonly(xs::Accumulator{Node,Int}) = Set(Pair(n.value,c) for (n,c) = xs)

const NOTFOUND = "NOT FOUND"

"""
Solver of a ten puzzle
"""
function solve!(xs::Accumulator{Node,Int};
                tgt=Rational(10), failed=Set{Set{Pair{Rational,Int}}}())::String
    valueonly(xs) in failed && return NOTFOUND

    if (sum(xs)==1)
        node = collect(keys(xs))[1]
        return node.value==tgt ? toexpr(node) : NOTFOUND
    end
    for ys = Reductions(xs)
        ans = solve!(ys, tgt=tgt)
        ans == NOTFOUND ? push!(failed, valueonly(ys)) : return ans
    end
    NOTFOUND
end
#===================================================================================#

"""
$(SIGNATURES)

Make ten (or given target value) using given numbers and four arithmetic operations

#Examples
```julia-repl
julia> maketen(1,1,9,9)
"(1 / 9 + 1) * 9"
```
"""
maketen(xs::Union{Int,Rational}...; tgt=10) = solve!(counter(Node.(xs)), tgt=Rational(tgt))


text:
mutable struct RadauIIA3ConstantCache{F,Tab,Tol,Dt,U,JType} <: OrdinaryDiffEqConstantCache
  uf::F
  tab::Tab
  κ::Tol
  ηold::Tol
  iter::Int
  cont1::U
  cont2::U
  cont3::U
  dtprev::Dt
  W_γdt::Dt
  status::DiffEqBase.NLStatus
  J::JType
end

function alg_cache(alg::RadauIIA3,u,rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},
                   ::Type{tTypeNoUnits},uprev,uprev2,f,t,dt,reltol,p,calck,::Val{false}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
  uf  = UDerivativeWrapper(f, t, p)
  uToltype = constvalue(uBottomEltypeNoUnits)
  tab = RadauIIA3Tableau(uToltype, constvalue(tTypeNoUnits))

  κ = convert(uToltype, 1//100)
  J = false .* _vec(rate_prototype) .* _vec(rate_prototype)'

  RadauIIA3ConstantCache(uf, tab, κ, one(uToltype), 10000, u, u, u, dt, dt, DiffEqBase.Convergence, J)
end

mutable struct RadauIIA3Cache{uType,cuType,uNoUnitsType,rateType,JType,W1Type,UF,JC,F1,Tab,Tol,Dt,rTol,aTol} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  z1::uType
  z2::uType
  w1::uType
  w2::uType
  dw12::cuType
  cubuff::cuType
  cont1::uType
  cont2::uType
  du1::rateType
  fsalfirst::rateType
  k::rateType
  k2::rateType
  fw1::rateType
  fw2::rateType
  J::JType
  W1::W1Type
  uf::UF
  tab::Tab
  κ::Tol
  ηold::Tol
  iter::Int
  tmp::uType
  atmp::uNoUnitsType
  jac_config::JC
  linsolve::F1
  rtol::rTol
  atol::aTol
  dtprev::Dt
  W_γdt::Dt
  status::DiffEqBase.NLStatus
end

function alg_cache(alg::RadauIIA3,u,rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},
                   ::Type{tTypeNoUnits},uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
  uf = UJacobianWrapper(f, t, p)
  uToltype = constvalue(uBottomEltypeNoUnits)
  tab = RadauIIA3Tableau(uToltype, constvalue(tTypeNoUnits))

  κ = alg.κ !== nothing ? convert(uToltype, alg.κ) : convert(uToltype, 1//100)

  z1 = zero(u); z2 = zero(u);
  w1 = zero(u); w2 = zero(u);
  dw12 = similar(u, Complex{eltype(u)})
  cubuff = similar(u, Complex{eltype(u)})
  cont1 = zero(u); cont2 = zero(u);

  fsalfirst = similar(rate_prototype)
  k = similar(rate_prototype); k2 = similar(rate_prototype);
  fw1 = similar(rate_prototype); fw2 = similar(rate_prototype);

  J, W1 = build_J_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits,Val(true))
  W1 = similar(J, Complex{eltype(W1)})

  du1 = similar(rate_prototype)

  tmp = zero(u)
  atmp = similar(u,uEltypeNoUnits)
  jac_config = jac_config = build_jac_config(alg, f, uf, du1, uprev, u, tmp, dw12)

  linprob = LinearProblem(W1,_vec(cubuff); u0=_vec(dw12))
  linsolve = init(linprob,alg.linsolve,alias_A=true,alias_b=true)
                   #Pl = LinearSolve.InvPreconditioner(Diagonal(_vec(weight))),
                   #Pr = Diagonal(_vec(weight)))

  rtol = reltol isa Number ? reltol : similar(reltol)
  atol = reltol isa Number ? reltol : similar(reltol)

  RadauIIA3Cache(u, uprev,
                 z1, z2, w1, w2,
                 dw12, cubuff, cont1, cont2,
                 du1, fsalfirst, k, k2, fw1, fw2,
                 J, W1,
                 uf, tab, κ, one(uToltype), 10000,
                 tmp, atmp, jac_config, linsolve, rtol, atol, dt, dt, DiffEqBase.Convergence)
end

mutable struct RadauIIA5ConstantCache{F,Tab,Tol,Dt,U,JType} <: OrdinaryDiffEqConstantCache
  uf::F
  tab::Tab
  κ::Tol
  ηold::Tol
  iter::Int
  cont1::U
  cont2::U
  cont3::U
  dtprev::Dt
  W_γdt::Dt
  status::DiffEqBase.NLStatus
  J::JType
end

function alg_cache(alg::RadauIIA5,u,rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},
                   ::Type{tTypeNoUnits},uprev,uprev2,f,t,dt,reltol,p,calck,::Val{false}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
  uf  = UDerivativeWrapper(f, t, p)
  uToltype = constvalue(uBottomEltypeNoUnits)
  tab = RadauIIA5Tableau(uToltype, constvalue(tTypeNoUnits))

  κ = alg.κ !== nothing ? convert(uToltype, alg.κ) : convert(uToltype, 1//100)
  J = false .* _vec(rate_prototype) .* _vec(rate_prototype)'

  RadauIIA5ConstantCache(uf, tab, κ, one(uToltype), 10000, u, u, u, dt, dt, DiffEqBase.Convergence, J)
end

mutable struct RadauIIA5Cache{uType,cuType,uNoUnitsType,rateType,JType,W1Type,W2Type,UF,JC,F1,F2,Tab,Tol,Dt,rTol,aTol} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  z1::uType
  z2::uType
  z3::uType
  w1::uType
  w2::uType
  w3::uType
  dw1::uType
  ubuff::uType
  dw23::cuType
  cubuff::cuType
  cont1::uType
  cont2::uType
  cont3::uType
  du1::rateType
  fsalfirst::rateType
  k::rateType
  k2::rateType
  k3::rateType
  fw1::rateType
  fw2::rateType
  fw3::rateType
  J::JType
  W1::W1Type
  W2::W2Type # complex
  uf::UF
  tab::Tab
  κ::Tol
  ηold::Tol
  iter::Int
  tmp::uType
  atmp::uNoUnitsType
  jac_config::JC
  linsolve1::F1
  linsolve2::F2
  rtol::rTol
  atol::aTol
  dtprev::Dt
  W_γdt::Dt
  status::DiffEqBase.NLStatus
end

function alg_cache(alg::RadauIIA5,u,rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},
                   ::Type{tTypeNoUnits},uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
  uf = UJacobianWrapper(f, t, p)
  uToltype = constvalue(uBottomEltypeNoUnits)
  tab = RadauIIA5Tableau(uToltype, constvalue(tTypeNoUnits))

  κ = alg.κ !== nothing ? convert(uToltype, alg.κ) : convert(uToltype, 1//100)

  z1 = zero(u); z2 = zero(u); z3 = zero(u)
  w1 = zero(u); w2 = zero(u); w3 = zero(u)
  dw1 = zero(u); ubuff = zero(u)
  dw23 = similar(u, Complex{eltype(u)}); cubuff = similar(u, Complex{eltype(u)})
  cont1 = zero(u); cont2 = zero(u); cont3 = zero(u)

  fsalfirst = similar(rate_prototype)
  k = similar(rate_prototype); k2 = similar(rate_prototype); k3 = similar(rate_prototype)
  fw1 = similar(rate_prototype); fw2 = similar(rate_prototype); fw3 = similar(rate_prototype)

  J, W1 = build_J_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits,Val(true))
  W2 = similar(J, Complex{eltype(W1)})

  du1 = similar(rate_prototype)

  tmp = zero(u)
  atmp = similar(u,uEltypeNoUnits)
  jac_config = build_jac_config(alg, f, uf, du1, uprev, u, tmp, dw1)

  linprob = LinearProblem(W1,_vec(ubuff); u0=_vec(dw1))
  linsolve1 = init(linprob,alg.linsolve,alias_A=true,alias_b=true)
                   #Pl = LinearSolve.InvPreconditioner(Diagonal(_vec(weight))),
                   #Pr = Diagonal(_vec(weight)))
  linprob = LinearProblem(W2,_vec(cubuff); u0=_vec(dw23))
  linsolve2 = init(linprob,alg.linsolve,alias_A=true,alias_b=true)
                   #Pl = LinearSolve.InvPreconditioner(Diagonal(_vec(weight))),
                   #Pr = Diagonal(_vec(weight)))

  rtol = reltol isa Number ? reltol : similar(reltol)
  atol = reltol isa Number ? reltol : similar(reltol)

  RadauIIA5Cache(u, uprev,
                 z1, z2, z3, w1, w2, w3,
                 dw1, ubuff, dw23, cubuff, cont1, cont2, cont3,
                 du1, fsalfirst, k, k2, k3, fw1, fw2, fw3,
                 J, W1, W2,
                 uf, tab, κ, one(uToltype), 10000,
                 tmp, atmp, jac_config, linsolve1, linsolve2, rtol, atol, dt, dt, DiffEqBase.Convergence)
end

ext: jl


text:
function readFasta(path)
    f = open(path, "r")

    seq = ""
    for line in readlines(f)
        if startswith(line, ">")
            continue
        else
            seq *= line
        end
    end
    close(f)
    return seq
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module Libc
@doc """
Interface to libc, the C standard library.
""" Libc

import Base: transcode, windowserror, show
import Core.Intrinsics: bitcast

export FILE, TmStruct, strftime, strptime, getpid, gethostname, free, malloc, calloc, realloc,
    errno, strerror, flush_cstdio, systemsleep, time, transcode
if Sys.iswindows()
    export GetLastError, FormatMessage
end

include(string(length(Core.ARGS) >= 2 ? Core.ARGS[2] : "", "errno_h.jl"))  # include($BUILDROOT/base/errno_h.jl)

## RawFD ##

# Wrapper for an OS file descriptor (on both Unix and Windows)
"""
    RawFD

Primitive type which wraps the native OS file descriptor.
`RawFD`s can be passed to methods like [`stat`](@ref) to
discover information about the underlying file, and can
also be used to open streams, with the `RawFD` describing
the OS file backing the stream.
"""
primitive type RawFD 32 end
RawFD(fd::Integer) = bitcast(RawFD, Cint(fd))
RawFD(fd::RawFD) = fd
Base.cconvert(::Type{Cint}, fd::RawFD) = bitcast(Cint, fd)

dup(x::RawFD) = ccall((@static Sys.iswindows() ? :_dup : :dup), RawFD, (RawFD,), x)
dup(src::RawFD, target::RawFD) = systemerror("dup", -1 ==
    ccall((@static Sys.iswindows() ? :_dup2 : :dup2), Int32,
                (RawFD, RawFD), src, target))

show(io::IO, fd::RawFD) = print(io, "RawFD(", bitcast(UInt32, fd), ')')  # avoids invalidation via show_default

# Wrapper for an OS file descriptor (for Windows)
if Sys.iswindows()
    primitive type WindowsRawSocket sizeof(Ptr) * 8 end # On Windows file descriptors are HANDLE's and 64-bit on 64-bit Windows
    WindowsRawSocket(handle::Ptr{Cvoid}) = bitcast(WindowsRawSocket, handle)
    WindowsRawSocket(handle::WindowsRawSocket) = handle

    Base.cconvert(::Type{Ptr{Cvoid}}, fd::WindowsRawSocket) = bitcast(Ptr{Cvoid}, fd)
    _get_osfhandle(fd::RawFD) = ccall(:_get_osfhandle, WindowsRawSocket, (RawFD,), fd)
    _get_osfhandle(fd::WindowsRawSocket) = fd
    function dup(src::WindowsRawSocket)
        new_handle = Ref(WindowsRawSocket(Ptr{Cvoid}(-1)))
        my_process = ccall(:GetCurrentProcess, stdcall, Ptr{Cvoid}, ())
        DUPLICATE_SAME_ACCESS = 0x2
        status = ccall(:DuplicateHandle, stdcall, Int32,
            (Ptr{Cvoid}, WindowsRawSocket, Ptr{Cvoid}, Ptr{WindowsRawSocket}, UInt32, Int32, UInt32),
            my_process, src, my_process, new_handle, 0, false, DUPLICATE_SAME_ACCESS)
        windowserror("dup failed", status == 0)
        return new_handle[]
    end
    function dup(src::WindowsRawSocket, target::RawFD)
        fd = ccall(:_open_osfhandle, RawFD, (WindowsRawSocket, Int32), dup(src), 0)
        dup(fd, target)
        ccall(:_close, Int32, (RawFD,), fd)
        nothing
    end

else
    _get_osfhandle(fd::RawFD) = fd
end

## FILE (not auto-finalized) ##

struct FILE
    ptr::Ptr{Cvoid}
end

modestr(s::IO) = modestr(isreadable(s), iswritable(s))
modestr(r::Bool, w::Bool) = r ? (w ? "r+" : "r") : (w ? "w" : throw(ArgumentError("neither readable nor writable")))

function FILE(fd::RawFD, mode)
    FILEp = ccall((@static Sys.iswindows() ? :_fdopen : :fdopen), Ptr{Cvoid}, (Cint, Cstring), fd, mode)
    systemerror("fdopen", FILEp == C_NULL)
    FILE(FILEp)
end

function FILE(s::IO)
    f = FILE(dup(RawFD(fd(s))),modestr(s))
    seek(f, position(s))
    f
end

Base.unsafe_convert(T::Union{Type{Ptr{Cvoid}},Type{Ptr{FILE}}}, f::FILE) = convert(T, f.ptr)
Base.close(f::FILE) = systemerror("fclose", ccall(:fclose, Cint, (Ptr{Cvoid},), f.ptr) != 0)

function Base.seek(h::FILE, offset::Integer)
    systemerror("fseek", ccall(:fseek, Cint, (Ptr{Cvoid}, Clong, Cint),
                               h.ptr, offset, 0) != 0)
    h
end

Base.position(h::FILE) = ccall(:ftell, Clong, (Ptr{Cvoid},), h.ptr)

# flush C stdio output from external libraries

"""
    flush_cstdio()

Flushes the C `stdout` and `stderr` streams (which may have been written to by external C code).
"""
flush_cstdio() = ccall(:jl_flush_cstdio, Cvoid, ())

## time-related functions ##

# TODO: check for usleep errors?
if Sys.isunix()
    systemsleep(s::Real) = ccall(:usleep, Int32, (UInt32,), round(UInt32, s*1e6))
elseif Sys.iswindows()
    function systemsleep(s::Real)
        ccall(:Sleep, stdcall, Cvoid, (UInt32,), round(UInt32, s * 1e3))
        return Int32(0)
    end
else
    error("systemsleep undefined for this OS")
end
"""
    systemsleep(s::Real)

Suspends execution for `s` seconds.
This function does not yield to Julia's scheduler and therefore blocks
the Julia thread that it is running on for the duration of the sleep time.

See also [`sleep`](@ref).
"""
systemsleep

struct TimeVal
   sec::Int64
   usec::Int64
end

function TimeVal()
    tv = Ref{TimeVal}()
    status = ccall(:jl_gettimeofday, Cint, (Ref{TimeVal},), tv)
    status != 0 && error("unable to determine current time: ", status)
    return tv[]
end

"""
    TmStruct([seconds])

Convert a number of seconds since the epoch to broken-down format, with fields `sec`, `min`,
`hour`, `mday`, `month`, `year`, `wday`, `yday`, and `isdst`.
"""
mutable struct TmStruct
    sec::Int32
    min::Int32
    hour::Int32
    mday::Int32
    month::Int32
    year::Int32
    wday::Int32
    yday::Int32
    isdst::Int32
    # on some platforms the struct is 14 words, even though 9 are specified
    _10::Int32
    _11::Int32
    _12::Int32
    _13::Int32
    _14::Int32

    TmStruct(sec, min, hour, mday, month, year, wday, yday, isdst) =
        new(sec, min, hour, mday, month, year, wday, yday, isdst, 0,0,0,0,0)
    TmStruct() = new(0,0,0,0,0,0,0,0,0,0,0,0,0,0)
    function TmStruct(t::Real)
        t = floor(t)
        tm = TmStruct()
        # TODO: add support for UTC via gmtime_r()
        ccall(:localtime_r, Ptr{TmStruct}, (Ref{Int}, Ref{TmStruct}), t, tm)
        return tm
    end
end

"""
    strftime([format], time)

Convert time, given as a number of seconds since the epoch or a `TmStruct`, to a formatted
string using the given format. Supported formats are the same as those in the standard C
library.
"""
strftime(t) = strftime("%c", t)
strftime(fmt::AbstractString, t::Real) = strftime(fmt, TmStruct(t))
# Use wcsftime instead of strftime to support different locales
function strftime(fmt::AbstractString, tm::TmStruct)
    wctimestr = Vector{Cwchar_t}(undef, 128)
    n = ccall(:wcsftime, Csize_t, (Ptr{Cwchar_t}, Csize_t, Cwstring, Ref{TmStruct}),
              wctimestr, length(wctimestr), fmt, tm)
    n == 0 && return ""
    return transcode(String, resize!(wctimestr, n))
end

"""
    strptime([format], timestr)

Parse a formatted time string into a `TmStruct` giving the seconds, minute, hour, date, etc.
Supported formats are the same as those in the standard C library. On some platforms,
timezones will not be parsed correctly. If the result of this function will be passed to
`time` to convert it to seconds since the epoch, the `isdst` field should be filled in
manually. Setting it to `-1` will tell the C library to use the current system settings to
determine the timezone.
"""
strptime(timestr::AbstractString) = strptime("%c", timestr)
function strptime(fmt::AbstractString, timestr::AbstractString)
    tm = TmStruct()
    r = ccall(:strptime, Cstring, (Cstring, Cstring, Ref{TmStruct}), timestr, fmt, tm)
    # the following would tell mktime() that this is a local time, and that
    # it should try to guess the timezone. not sure if/how this should be
    # exposed in the API.
    # tm.isdst = -1
    if r == C_NULL
        # TODO: better error message
        throw(ArgumentError("invalid arguments"))
    end
    @static if Sys.isapple()
        # if we didn't explicitly parse the weekday or year day, use mktime
        # to fill them in automatically.
        if !occursin(r"([^%]|^)%(a|A|j|w|Ow)", fmt)
            ccall(:mktime, Int, (Ref{TmStruct},), tm)
        end
    end
    return tm
end

# system date in seconds

"""
    time(t::TmStruct)

Converts a `TmStruct` struct to a number of seconds since the epoch.
"""
time(tm::TmStruct) = Float64(ccall(:mktime, Int, (Ref{TmStruct},), tm))

"""
    time()

Get the system time in seconds since the epoch, with fairly high (typically, microsecond) resolution.
"""
time() = ccall(:jl_clock_now, Float64, ())

## process-related functions ##

"""
    getpid() -> Int32

Get Julia's process ID.
"""
getpid() = ccall(:jl_getpid, Int32, ())

## network functions ##

"""
    gethostname() -> AbstractString

Get the local machine's host name.
"""
function gethostname()
    hn = Vector{UInt8}(undef, 256)
    err = @static if Sys.iswindows()
        ccall(:gethostname, stdcall, Int32, (Ptr{UInt8}, UInt32), hn, length(hn))
    else
        ccall(:gethostname, Int32, (Ptr{UInt8}, UInt), hn, length(hn))
    end
    systemerror("gethostname", err != 0)
    return GC.@preserve hn unsafe_string(pointer(hn))
end

## system error handling ##

"""
    errno([code])

Get the value of the C library's `errno`. If an argument is specified, it is used to set the
value of `errno`.

The value of `errno` is only valid immediately after a `ccall` to a C library routine that
sets it. Specifically, you cannot call `errno` at the next prompt in a REPL, because lots of
code is executed between prompts.
"""
errno() = ccall(:jl_errno, Cint, ())
errno(e::Integer) = ccall(:jl_set_errno, Cvoid, (Cint,), e)

"""
    strerror(n=errno())

Convert a system call error code to a descriptive string
"""
strerror(e::Integer) = unsafe_string(ccall(:strerror, Cstring, (Int32,), e))
strerror() = strerror(errno())

"""
    GetLastError()

Call the Win32 `GetLastError` function [only available on Windows].
"""
function GetLastError end

"""
    FormatMessage(n=GetLastError())

Convert a Win32 system call error code to a descriptive string [only available on Windows].
"""
function FormatMessage end

if Sys.iswindows()
    GetLastError() = ccall(:GetLastError, stdcall, UInt32, ())

    FormatMessage(e) = FormatMessage(UInt32(e))
    function FormatMessage(e::UInt32=GetLastError())
        FORMAT_MESSAGE_ALLOCATE_BUFFER = UInt32(0x100)
        FORMAT_MESSAGE_FROM_SYSTEM = UInt32(0x1000)
        FORMAT_MESSAGE_IGNORE_INSERTS = UInt32(0x200)
        FORMAT_MESSAGE_MAX_WIDTH_MASK = UInt32(0xFF)
        lpMsgBuf = Ref{Ptr{UInt16}}()
        lpMsgBuf[] = 0
        len = ccall(:FormatMessageW, stdcall, UInt32, (UInt32, Ptr{Cvoid}, UInt32, UInt32, Ptr{Ptr{UInt16}}, UInt32, Ptr{Cvoid}),
                    FORMAT_MESSAGE_ALLOCATE_BUFFER | FORMAT_MESSAGE_FROM_SYSTEM | FORMAT_MESSAGE_IGNORE_INSERTS | FORMAT_MESSAGE_MAX_WIDTH_MASK,
                    C_NULL, e, 0, lpMsgBuf, 0, C_NULL)
        p = lpMsgBuf[]
        len == 0 && return ""
        buf = Vector{UInt16}(undef, len)
        GC.@preserve buf unsafe_copyto!(pointer(buf), p, len)
        ccall(:LocalFree, stdcall, Ptr{Cvoid}, (Ptr{Cvoid},), p)
        return transcode(String, buf)
    end
end

## Memory related ##

"""
    free(addr::Ptr)

Call `free` from the C standard library. Only use this on memory obtained from [`malloc`](@ref), not
on pointers retrieved from other C libraries. [`Ptr`](@ref) objects obtained from C libraries should
be freed by the free functions defined in that library, to avoid assertion failures if
multiple `libc` libraries exist on the system.
"""
free(p::Ptr) = ccall(:free, Cvoid, (Ptr{Cvoid},), p)

"""
    malloc(size::Integer) -> Ptr{Cvoid}

Call `malloc` from the C standard library.
"""
malloc(size::Integer) = ccall(:malloc, Ptr{Cvoid}, (Csize_t,), size)

"""
    realloc(addr::Ptr, size::Integer) -> Ptr{Cvoid}

Call `realloc` from the C standard library.

See warning in the documentation for [`free`](@ref) regarding only using this on memory originally
obtained from [`malloc`](@ref).
"""
realloc(p::Ptr, size::Integer) = ccall(:realloc, Ptr{Cvoid}, (Ptr{Cvoid}, Csize_t), p, size)

"""
    calloc(num::Integer, size::Integer) -> Ptr{Cvoid}

Call `calloc` from the C standard library.
"""
calloc(num::Integer, size::Integer) = ccall(:calloc, Ptr{Cvoid}, (Csize_t, Csize_t), num, size)

free(p::Cstring) = free(convert(Ptr{UInt8}, p))
free(p::Cwstring) = free(convert(Ptr{Cwchar_t}, p))

## Random numbers ##

# To limit dependency on rand functionality implemented in the Random module,
# Libc.rand is used in file.jl, and could be used in error.jl (but it breaks a test)
"""
    rand([T::Type])

Interface to the C `rand()` function. If `T` is provided, generate a value of type `T`
by composing two calls to `rand()`. `T` can be `UInt32` or `Float64`.
"""
rand() = ccall(:rand, Cint, ())
@static if Sys.iswindows()
    # Windows RAND_MAX is 2^15-1
    rand(::Type{UInt32}) = ((rand() % UInt32) << 17) ⊻ ((rand() % UInt32) << 8) ⊻ (rand() % UInt32)
else
    # RAND_MAX is at least 2^15-1 in theory, but we assume 2^16-1
    # on non-Windows systems (in practice, it's 2^31-1)
    rand(::Type{UInt32}) = ((rand() % UInt32) << 16) ⊻ (rand() % UInt32)
end
rand(::Type{Float64}) = rand(UInt32) * 2.0^-32

"""
    srand([seed])

Interface to the C `srand(seed)` function.
"""
srand(seed=Base._make_uint_seed()) = ccall(:srand, Cvoid, (Cuint,), seed)

struct Cpasswd
   username::Cstring
   uid::Culong
   gid::Culong
   shell::Cstring
   homedir::Cstring
   gecos::Cstring
   Cpasswd() = new(C_NULL, typemax(Culong), typemax(Culong), C_NULL, C_NULL, C_NULL)
end
mutable struct Cgroup
    groupname::Cstring # group name
    gid::Culong        # group ID
    mem::Ptr{Cstring}  # group members
    Cgroup() = new(C_NULL, typemax(Culong), C_NULL)
end
struct Passwd
    username::String
    uid::UInt
    gid::UInt
    shell::String
    homedir::String
    gecos::String
end
struct Group
    groupname::String
    gid::UInt
    mem::Vector{String}
end

function getpwuid(uid::Unsigned, throw_error::Bool=true)
    ref_pd = Ref(Cpasswd())
    ret = ccall(:uv_os_get_passwd2, Cint, (Ref{Cpasswd}, Culong), ref_pd, uid)
    if ret != 0
        throw_error && Base.uv_error("getpwuid", ret)
        return
    end
    pd = ref_pd[]
    pd = Passwd(
        pd.username == C_NULL ? "" : unsafe_string(pd.username),
        pd.uid,
        pd.gid,
        pd.shell == C_NULL ? "" : unsafe_string(pd.shell),
        pd.homedir == C_NULL ? "" : unsafe_string(pd.homedir),
        pd.gecos == C_NULL ? "" : unsafe_string(pd.gecos),
    )
    ccall(:uv_os_free_passwd, Cvoid, (Ref{Cpasswd},), ref_pd)
    return pd
end
function getgrgid(gid::Unsigned, throw_error::Bool=true)
    ref_gp = Ref(Cgroup())
    ret = ccall(:uv_os_get_group, Cint, (Ref{Cgroup}, Culong), ref_gp, gid)
    if ret != 0
        throw_error && Base.uv_error("getgrgid", ret)
        return
    end
    gp = ref_gp[]
    members = String[]
    if gp.mem != C_NULL
        while true
            mem = unsafe_load(gp.mem, length(members) + 1)
            mem == C_NULL && break
            push!(members, unsafe_string(mem))
        end
    end
    gp = Group(
         gp.groupname == C_NULL ? "" : unsafe_string(gp.groupname),
         gp.gid,
         members,
    )
    ccall(:uv_os_free_group, Cvoid, (Ref{Cgroup},), ref_gp)
    return gp
end

getuid() = ccall(:jl_getuid, Culong, ())
geteuid() = ccall(:jl_geteuid, Culong, ())

# Include dlopen()/dlpath() code
include("libdl.jl")
using .Libdl

end # module

max_stars_repo_name: lawless-m/julia

max_issues_repo_name: lawless-m/julia

max_forks_repo_name: lawless-m/julia


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

using JuliaFEM, LinearAlgebra, Test

# compare simple 3d heat problem to analytical solution

function calc_3d_heat_model(mesh_name)
    fn = @__DIR__() * "/testdata/rod_short.med"
    mesh = aster_read_mesh(fn, mesh_name)
    problem = Problem(Heat, "rod", 1)
    bc = Problem(Dirichlet, "left support T=100", 1, "temperature")
    problem_elements = create_elements(mesh, "ROD", "FACE2")
    bc_elements = create_elements(mesh, "FACE1")
    update!(problem_elements, "thermal conductivity", 100.0)
    update!(problem_elements, "external temperature", 0.0)
    update!(problem_elements, "heat transfer coefficient", 1000.0)
    update!(bc_elements, "temperature 1", 100.0)
    add_elements!(problem, problem_elements)
    add_elements!(bc, bc_elements)
    analysis = Analysis(Linear)
    add_problems!(analysis, problem, bc)
    run!(analysis)
    time = 0.0
    T = problem("temperature", time)
    T_min = minimum(map(first, values(T)))
    return T_min
end

for model in ["Tet4", "Tet10", "Hex8", "Hex20", "Hex27"]
    Tmin = calc_3d_heat_model(model)
    Tacc = 100/3
    rtol = norm(Tmin-Tacc)/max(Tmin,Tacc)*100.0
    @debug("Temperature for model $model", model, Tmin, Tacc, rtol)
    @test isapprox(Tmin, 100/3)
end


text:
## Computation options

"""
    Bonding

Selection mode for the detection of bonds. The choices are:
- `Input`: use the input bonds. Fail if those are not specified.
- `Guess`: guess bonds using a variant of chemfiles / VMD algorithm.
- `Auto`: if the input specifies bonds, use them unless they look suspicious (too small or
  or too large according to a heuristic). Otherwise, fall back to `Guess`.
"""
module Bonding
    @enum _Bonding begin
        Input
        Guess
        Auto
    end
    """See help for [`Bonding`](@ref)"""
    Input, Guess, Auto
end
import .Bonding

module StructureType
    @enum _StructureType begin
        Auto
        MOF
        Cluster
        Zeolite
        Guess
    end
    """See help for [`StructureType`](@ref)"""
    Auto, MOF, Cluster, Zeolite, Guess
    export MOF, Cluster, Zeolite
end
import .StructureType
import .StructureType: _StructureType

module Clustering
    @enum _Clustering begin
        Auto         = 1
        Input        = 2
        EachVertex   = 3
        SingleNodes  = 4
        AllNodes     = 5
        Standard     = 6
        PE           = 7
        PEM          = 8
    end
    """See help for [`Clustering`](@ref)"""
    Auto, Input, EachVertex, SingleNodes, AllNodes, Standard, PE, PEM
    export EachVertex, SingleNodes, AllNodes, Standard, PE, PEM
end
import .Clustering
import .Clustering: _Clustering

function clustering_from_num(x::Integer)
    if x == 1
        return Clustering.Auto
    elseif x == 2
        return Clustering.Input
    elseif x == 3
        return Clustering.EachVertex
    elseif x == 4
        return Clustering.SingleNodes
    elseif x == 5
        return Clustering.AllNodes
    elseif x == 6
        return Clustering.Standard
    elseif x == 7
        return Clustering.PE
    elseif x == 8
        return Clustering.PEM
    end
    throw(ArgumentError(lazy"No clustering from number $x"))
end

function clustering_from_symb(x::Symbol)
    if x === :Auto
        return Clustering.Auto
    elseif x === :Input
        return Clustering.Input
    elseif x === :EachVertex
        return Clustering.EachVertex
    elseif x === :SingleNodes
        return Clustering.SingleNodes
    elseif x === :AllNodes
        return Clustering.AllNodes
    elseif x === :Standard
        return Clustering.Standard
    elseif x === :PE
        return Clustering.PEM
    elseif x === :PEM
        return Clustering.PEM
    end
    throw(ArgumentError(lazy"No clustering from symbol $x"))
end

function Base.parse(::Type{_Clustering}, s::AbstractString)
    if s == "Auto"
        return Clustering.Auto
    elseif s == "Input"
        return Clustering.Input
    elseif s == "EachVertex"
        return Clustering.EachVertex
    elseif s == "SingleNodes"
        return Clustering.SingleNodes
    elseif s == "AllNodes"
        return Clustering.AllNodes
    elseif s == "Standard"
        return Clustering.Standard
    elseif s == "PE"
        return Clustering.PE
    elseif s == "PEM"
        return Clustering.PEM
    end
    throw(ArgumentError(lazy"No clustering from string $x"))
end

"""
    StructureType

Selection mode for the crystal structure. This choice impacts the bond detection algorithm
as well as the clustering algorithm used.

The choices are:
- `Auto`: No specific structure information. Use Van der Waals radii for bond detection and
  `Input` as [`Clustering`](@ref), or `EachVertex` if the input does not provide residues.
- `MOF`: Use Van der Waals radii for non-metallic atoms and larger radii for metals. Detect
  organic and inorganic clusters and subdivide them according to `AllNodes` and
  `SingleNodes` to identify underlying nets.
- `Cluster`: similar to MOF but metallic atoms are not given a wider radius.
- `Zeolite`: Same as `Auto` but use larger radii for metals (and metalloids) and attempt
  to enforce that each O atom has exactly two neighbours and that they are not O atoms.
- `Guess`: try to identify the clusters as in `Cluster`. If it fails, fall back to `Auto`.
"""
StructureType

"""
    Clustering

The clustering algorithm used to group atoms into vertices.

This choice only affects the creation of a `UnderlyingNets` from a `Crystal`, not the
`Crystal` itself, and in particular not the bond detection algorithm.

The basic choices are:
- `Auto`: determined using the [`StructureType`](@ref).
- `Input`: use the input residues as vertices. Fail if some atom does not belong to a
  residue.
- `EachVertex`: each atom is its own vertex. Vertices with degree 2 or lower are
  iteratively collapsed into edges until all vertices have degree 3 or more.

The next clustering options are designed for MOFs but may target other kinds of frameworks.
In all cases, the clusters are refinements on top of already-defined clusters, such as the
organic and inorganic SBUs defined by the `MOF` structure.
Except for `AllNodes`, infinite clusters (such as the inorganic clusters in a rod MOF) are
split into new finite clusters using heuristics.
- `SingleNodes`: each already-defined cluster is mapped to a vertex.
- `AllNodes`: keep points of extension for organic clusters.
- `Standard`: make each metallic atom its own vertex and do not bond those together if they
  share a common non-metallic neighbour.
- `PE`: stands for Points of Extension. Keep points of extension for organic clusters,
  remove metallic centres and bond their surrounding points of extension.
- `PEM`: stands for Points of Extension and Metals. Keep points of extension for organic
  clusters and each metal centre as a separate vertex.
"""
Clustering

"""
    ClusterKinds(sbus, toclassify=Int[])

Description of the different kinds of SBUs there should be when making clusters.

`sbus` should be a list of set of symbols, each set containing the different elements
acceptable in this SBU (an empty set designates all remaining elements). All elements of
the same category of the periodic table can be grouped together by putting the name of the
category.
For example, `ClusterKinds([[:Au, :halogen, :nonmetal], [:metal, :metalloid], []])` means
that there are three kinds of SBUs:
- the first kind can only hold halogens, non-metals and Au atoms
- the second kind can only hold metalloids and metals (except Au)
- the third kind can hold all the other elements.

The list of possible categories is: :actinide, :noble (for noble gas), :halogen,
:lanthanide, :metal, :metalloid and :nonmetal.

`toclassify` contains the list of SBUs which are not actual SBUs but only groups of
atoms waiting to be merged to a neighboring SBU. The neighboring SBU is chosen
by order in the `sbus` list.

The cluster kinds used by default are
`CrystalNets.ClusterKinds([[:metal, :actinide, :lanthanide], [:C,], [:P, :S],
                           [:nonmetal, :metalloid, :halogen], [:noble]], [3, 4])`.
This means that all atoms that are either metals, actinides or lanthanides are assigned to
class 1 and all C atoms in SBUs of class 2.
Afterwards, each group of adjacent P or S atoms is assigned either class 1 if any of its
neighbor is of class 1, or class 2 otherwise if any of its neighbor is of class 2.
If no such neighbor exist, it is assigned to class 1.
Finally, each group of adjacent nonmetals, metalloids and halogens is assigned class 1 or 2
following the same rule as for P and S atoms.

At the end of the procedure, all atoms are thus given a class between `1` and `length(sbus)`
which is not in `toclassify`. See also [`find_sbus`](@ref) for the implementation of this
procedure.

To determine which SBU kind corresponds to a given atom, use `getindex`:
```jldoctest
julia> sbu_kinds = CrystalNets.ClusterKinds([[:nonmetal, :halogen], [:metal, :F]]);

julia> sbu_kinds[:O] # nonmetal
1

julia> sbu_kinds[:Au] # metal
2

julia> sbu_kinds[:F] # specifically F
2

julia> sbu_kinds[:Ne] # no given SBU kind
0
```
If no empty set has been explicitly added to `sbus` and an element falls outside of the
included categories, the returned SBU kind is 0.

An exception is made for nonmetals which are part of an aromatic heterocycle: those will be
treated separately and put in the SBU of the corresponding carbons.
"""
struct ClusterKinds
    dict::Dict{Symbol,Int}
    default::Int
    tomerge::Vector{Int}
    len::Int

    function ClusterKinds(sbus, tomerge=Int[])
        dict = Dict{Symbol,Int}()
        default = 0
        for (i, sbu) in enumerate(sbus)
            if isempty(sbu)
                if default != 0
                    throw(ArgumentError("ClusterKinds cannot contain multiple empty sets."))
                end
                default = i
                continue
            end
            for x in sbu
                old = get!(dict, x, i)
                if old != i
                    throw(ArgumentError(lazy"Element \"$x\" appears in at least two SBU kinds ($old and $i)."))
                end
            end
        end
        return new(dict, default, tomerge, length(sbus))
    end
end

function Base.getindex(sbus::ClusterKinds, i::Int)
    get(sbus.dict, element_categories[i], sbus.default)
end
function Base.getindex(sbus::ClusterKinds, x::Symbol)
    ret = get(sbus.dict, x, 0)
    ret == 0 || return ret
    num = get(atomic_numbers, x, 0)
    num == 0 && throw(MissingAtomInformation(lazy"unknown atom name: $x."))
    return sbus[num]
end
function getmetal(sbus::ClusterKinds)
    ret = get(sbus.dict, :metal, nothing)
    ret isa Int && return ret
    ret = get(sbus.dict, :Fe, nothing)
    ret isa Int && return ret
    ret = get(sbus.dict, :Cu, nothing)
    ret isa Int && return ret
    for (sym, at) in atomic_numbers
        ismetal[at] || continue
        ret = get(sbus.dict, sym, nothing)
        ret isa Int && return ret
    end
    return 0
end

Base.length(sbus::ClusterKinds) = sbus.len

const default_sbus = ClusterKinds([
    [:metal, :actinide, :lanthanide], [:C, :Pc, :Ss], [:P, :S],
    [:nonmetal, :metalloid, :halogen], [:noble]
], [3, 4])
# Note: Pc (respectively Ss) is the internal designation of an organic P (resp. S) atom

function ifbooltempdirorempty(x)::String
    if x isa Bool
        x ? tempdir() : ""
    else
        x
    end
end

"""
    Options

Different options, passed as keyword arguments.

## Basic options
- `name`: a name for the structure.
- `bonding`: one of the [`Bonding`](@ref) options. Default is `Bonding.Auto`.
- `structure`: one of the [`StructureType`](@ref) options. Default is `StructureType.Auto`.
- `clusterings`: a list of [`Clustering`](@ref) options. Default is `[Clustering.Auto]`.

## Exports
For each export option, the accepted values are either a string, indicating the path to
the directory in which to store the export, or a boolean, specifying whether or not to do
the export. If the value is `true`, a path will be automatically determined. An empty
string is equivalent to `false`.
- `export_input`: the parsed structure, as a .vtf
- `export_trimmed`: the parsed structure after iteratively removing all atoms having only one
  neighbour, as a .vtf
- `export_attributions`: the attribution of vertices into SBUs, as a .pdb. Only relevant for
  the `MOF` [`StructureType`](@ref).
- `export_clusters`: the clustering of vertices, as a .vtf
- `export_net`: the overall extracted net on which the topology is computed, as a .vtf.
- `export_subnets`: each connected component of the overall net as a separate .vtf file.
  These subnets are defined after grouping vertices according to their [`Clustering`](@ref).

## Other options
- `ignore_atoms`: set of atom symbols to ignore (for instance `[:C,:H]` will
  remove carbohydrate solvent residues).
- `ignore_types`: disregard atom types to compute the topology, making pcu and
  pcu-b identical for example (default is true)
- `cutoff_coeff`: coefficient used to detect bonds. Default is 0.75, higher
  values will include bonds that were considered too long before.
- `skip_minimize`: assume that the cell is already the unit cell (default is false).
- `dimensions`: the set of crystal net dimensions to consider. For instance, putting
  `Set(3)` will ensure that only 3-dimensional nets are considered.
  Default is `Set([1,2,3])`.
- `cluster_kinds`: a [`ClusterKinds`](@ref). Default separates organic and inorganic SBUs.
- `ignore_homoatomic_bonds`: a `Set{Symbol}` such that all X-X bonds of the net are removed
  if X is an atom whose type is in `ignore_homoatomic_bonds`.
- `max_polyhedron_radius`: an integer specifying the maximum number of bonds between two
  corners of the coordination polyhedron built for the [`Clustering.PE`](@ref Clustering) option.
  Default is 4.

## Miscellaneous options
These boolean options have a default value that may be determined by [`Bonding`](@ref),
[`StructureType`](@ref) and [`Clustering`](@ref). They can be directly overriden here.
- `bond_adjacent_sbus`: bond together SBUs which are only separated by a single C atom.
- `authorize_pruning`: remove colliding atoms in the input. Default is true.
- `wider_metallic_bonds`: for bond detections, metals have a radius equal to 1.5× their Van
  der Waals radius. Default is false, unless [`StructureType`](@ref) is `MOF` or `Zeolite`.
- `ignore_homometallic_bonds`: do not bond two metallic atoms of the same type if they are
  up to third neighbours anyway. Default is false, unless [`StructureType`](@ref) is `MOF`.
- `ignore_metal_cluster_bonds`: do not bond two metallic clusters together if they share at
  least one non-metallic neighbour. Default is false.
- `ignore_low_occupancy`: atoms with occupancy lower than 0.5 are ignored. Default is false.
- `detect_paddlewheels`: detect paddle-wheel pattern and group them into an inorganic vertex.
  Default is true.
- `detect_organiccycles`: detect organic cycles and collapse all belonging C atoms into a new
  vertex. Default is true.
- `detect_pe`: detect organic points-of-extension (organic atoms bonded to another SBU) and
  transform them into vertices. Default is true.
- `cluster_simple_pe`: cluster adjacent points-of-extension if they are not part of a cycle.
  Default is true.
- `separate_metals`: separate each metal atom into its own vertex (instead of grouping them
  to form metallic clusters if they are adjacent or bonded by an oxygen). Default is false,
  unless [`Clustering`](@ref) is `Standard` or `PEM`.
- `split_O_vertex`: if a vertex is composed of a single O, remove it and bond together all of
  its neighbors. Default is true.
- `unify_sbu_decomposition`: apply the same rule to decompose both periodic and finite SBUs.
  Default is false.
- `force_warn`: force printing warning and information even during `..._dataset` and
  `..._topologies` function calls. Default is false.

## Internal fields
These fields are for internal use and should not be modified by the user:
- `dryrun`: store information on possible options to try (for `guess_topology`).
- `_pos`: the positions of the centre of the clusters collapsed into vertices.
- `error`: store the first error that occured when building the net.
- `throw_error`: if set, throw the error instead of storing it in the `error` field.
"""
struct Options
    name::String # used for exports

    # Input options
    bonding::Bonding._Bonding
    cutoff_coeff::Float64
    structure::StructureType._StructureType
    authorize_pruning::Bool
    wider_metallic_bonds::Bool
    ignore_atoms::Set{Symbol}
    ignore_homoatomic_bonds::Set{Symbol}
    ignore_homometallic_bonds::Bool
    ignore_low_occupancy::Bool
    export_input::String
    export_trimmed::String
    force_warn::Bool

    # Clustering options
    clusterings::Vector{_Clustering}
    bond_adjacent_sbus::Bool
    ignore_metal_cluster_bonds::Union{Nothing,Bool}
    cluster_kinds::ClusterKinds
    detect_paddlewheels::Bool
    detect_organiccycles::Bool
    detect_pe::Bool
    cluster_simple_pe::Bool
    split_O_vertex::Bool
    unify_sbu_decomposition::Bool
    separate_metals::Union{Nothing,Bool}
    max_polyhedron_radius::Int
    export_attributions::String
    export_clusters::String

    # Topology computation options
    skip_minimize::Bool
    dimensions::Set{Int}
    ignore_types::Bool
    export_net::String
    export_subnets::String

    # Internal
    _pos::Vector{SVector{3,Float64}}
    dryrun::Union{Nothing,Dict{Symbol,Union{Nothing,Set{Symbol}}}}
    error::String
    throw_error::Bool

    function Options(; name="unnamed",
                       bonding=Bonding.Auto,
                       structure=StructureType.Auto,
                       cutoff_coeff=0.75,
                       wider_metallic_bonds=nothing,
                       authorize_pruning=true,
                       ignore_atoms=Set{Symbol}(),
                       ignore_homoatomic_bonds=Set{Symbol}(),
                       ignore_homometallic_bonds=nothing,
                       ignore_low_occupancy=false,
                       export_input=DOEXPORT[],
                       export_trimmed=false,
                       force_warn=false,
                       clusterings=_Clustering[Clustering.Auto],
                       bond_adjacent_sbus=false,
                       ignore_metal_cluster_bonds=nothing,
                       cluster_kinds=default_sbus,
                       detect_paddlewheels=true,
                       detect_organiccycles=true,
                       detect_pe=true,
                       cluster_simple_pe=true,
                       split_O_vertex=true,
                       unify_sbu_decomposition=false,
                       separate_metals=nothing,
                       max_polyhedron_radius=4,
                       export_attributions="",
                       export_clusters="",
                       skip_minimize=false,
                       dimensions=Set{Int}([1,2,3]),
                       ignore_types=true,
                       export_net=false,
                       export_subnets=DOEXPORT[],
                       _pos=SVector{3,Float64}[],
                       dryrun=nothing,
                       error="",
                       throw_error=false,
                    )

        _export_input = ifbooltempdirorempty(export_input)
        _export_trimmed = ifbooltempdirorempty(export_trimmed)
        _export_attributions = ifbooltempdirorempty(export_attributions)
        _export_clusters = ifbooltempdirorempty(export_clusters)
        _export_net = ifbooltempdirorempty(export_net)
        _export_subnets = ifbooltempdirorempty(export_subnets)

        _ignore_homometallic_bonds = if ignore_homometallic_bonds === nothing
            structure == StructureType.MOF
        else
            ignore_homometallic_bonds
        end

        _wider_metallic_bonds = if wider_metallic_bonds === nothing
            structure == StructureType.MOF || structure == StructureType.Zeolite
        else
            wider_metallic_bonds
        end

        new(
            name,
            bonding,
            cutoff_coeff,
            structure,
            _wider_metallic_bonds,
            authorize_pruning,
            Set{Symbol}(ignore_atoms),
            Set{Symbol}(ignore_homoatomic_bonds),
            _ignore_homometallic_bonds,
            ignore_low_occupancy,
            _export_input,
            _export_trimmed,
            force_warn,
            clusterings,
            bond_adjacent_sbus,
            ignore_metal_cluster_bonds,
            cluster_kinds,
            detect_paddlewheels,
            detect_organiccycles,
            detect_pe,
            cluster_simple_pe,
            split_O_vertex,
            unify_sbu_decomposition,
            separate_metals,
            max_polyhedron_radius,
            _export_attributions,
            _export_clusters,
            skip_minimize,
            Set{Int}(dimensions),
            ignore_types,
            _export_net,
            _export_subnets,
            _pos,
            dryrun,
            error,
            throw_error,
        )
    end
end

function Options(options::Options; kwargs...)
    isempty(kwargs) && return options
    base = Dict{Symbol,Any}([x => getfield(options, x) for x in fieldnames(Options)])
    for (kwarg, val) in kwargs
        T = fieldtype(Options, kwarg)
        val = if isconcretetype(T) && !(T <: Enum)
            if T <: Set
                union(base[kwarg], T(val))
            elseif T === String
                if kwarg === :error
                    ifelse(isempty(options.error), val::String, options.error)
                else
                    val::Union{String,Bool}
                end
            else
                T(val)
            end
        else
            val
        end
        base[kwarg] = val
    end
    return Options(; base...)
end


text:
export krann

import ..RANNModel

using mlpack._Internal.io

import mlpack_jll
const krannLibrary = mlpack_jll.libmlpack_julia_krann

# Call the C binding of the mlpack krann binding.
function krann_mlpackMain()
  success = ccall((:krann, krannLibrary), Bool, ())
  if !success
    # Throw an exception---false means there was a C++ exception.
    throw(ErrorException("mlpack binding error; see output"))
  end
end

" Internal module to hold utility functions. "
module krann_internal
  import ..krannLibrary

import ...RANNModel

# Get the value of a model pointer parameter of type RANNModel.
function IOGetParamRANNModel(paramName::String)::RANNModel
  RANNModel(ccall((:IO_GetParamRANNModelPtr, krannLibrary), Ptr{Nothing}, (Cstring,), paramName))
end

# Set the value of a model pointer parameter of type RANNModel.
function IOSetParamRANNModel(paramName::String, model::RANNModel)
  ccall((:IO_SetParamRANNModelPtr, krannLibrary), Nothing, (Cstring, Ptr{Nothing}), paramName, model.ptr)
end

# Serialize a model to the given stream.
function serializeRANNModel(stream::IO, model::RANNModel)
  buf_len = UInt[0]
  buf_ptr = ccall((:SerializeRANNModelPtr, krannLibrary), Ptr{UInt8}, (Ptr{Nothing}, Ptr{UInt}), model.ptr, Base.pointer(buf_len))
  buf = Base.unsafe_wrap(Vector{UInt8}, buf_ptr, buf_len[1]; own=true)
  write(stream, buf)
end
# Deserialize a model from the given stream.
function deserializeRANNModel(stream::IO)::RANNModel
  buffer = read(stream)
  RANNModel(ccall((:DeserializeRANNModelPtr, krannLibrary), Ptr{Nothing}, (Ptr{UInt8}, UInt), Base.pointer(buffer), length(buffer)))
end
end # module

"""
    krann(; [alpha, first_leaf_exact, input_model, k, leaf_size, naive, query, random_basis, reference, sample_at_leaves, seed, single_mode, single_sample_limit, tau, tree_type, verbose])

This program will calculate the k rank-approximate-nearest-neighbors of a set of
points. You may specify a separate set of reference points and query points, or
just a reference set which will be used as both the reference and query set. You
must specify the rank approximation (in %) (and optionally the success
probability).

For example, the following will return 5 neighbors from the top 0.1% of the data
(with probability 0.95) for each point in `input` and store the distances in
`distances` and the neighbors in `neighbors.csv`:

```julia
julia> using CSV
julia> input = CSV.read("input.csv")
julia> distances, neighbors, _ = krann(k=5, reference=input,
            tau=0.1)
```

Note that tau must be set such that the number of points in the corresponding
percentile of the data is greater than k.  Thus, if we choose tau = 0.1 with a
dataset of 1000 points and k = 5, then we are attempting to choose 5 nearest
neighbors out of the closest 1 point -- this is invalid and the program will
terminate with an error message.

The output matrices are organized such that row i and column j in the neighbors
output file corresponds to the index of the point in the reference set which is
the i'th nearest neighbor from the point in the query set with index j.  Row i
and column j in the distances output file corresponds to the distance between
those two points.

# Arguments

 - `alpha::Float64`: The desired success probability.  Default value
      `0.95`.

 - `first_leaf_exact::Bool`: The flag to trigger sampling only after
      exactly exploring the first leaf.  Default value `false`.

 - `input_model::RANNModel`: Pre-trained kNN model.
 - `k::Int`: Number of nearest neighbors to find.  Default value `0`.

 - `leaf_size::Int`: Leaf size for tree building (used for kd-trees, UB
      trees, R trees, R* trees, X trees, Hilbert R trees, R+ trees, R++ trees,
      and octrees).  Default value `20`.

 - `naive::Bool`: If true, sampling will be done without using a tree.
      Default value `false`.

 - `query::Array{Float64, 2}`: Matrix containing query points (optional).
 - `random_basis::Bool`: Before tree-building, project the data onto a
      random orthogonal basis.  Default value `false`.

 - `reference::Array{Float64, 2}`: Matrix containing the reference
      dataset.
 - `sample_at_leaves::Bool`: The flag to trigger sampling at leaves.
      Default value `false`.

 - `seed::Int`: Random seed (if 0, std::time(NULL) is used).  Default
      value `0`.

 - `single_mode::Bool`: If true, single-tree search is used (as opposed to
      dual-tree search.  Default value `false`.

 - `single_sample_limit::Int`: The limit on the maximum number of samples
      (and hence the largest node you can approximate).  Default value `20`.

 - `tau::Float64`: The allowed rank-error in terms of the percentile of
      the data.  Default value `5`.

 - `tree_type::String`: Type of tree to use: 'kd', 'ub', 'cover', 'r',
      'x', 'r-star', 'hilbert-r', 'r-plus', 'r-plus-plus', 'oct'.  Default value
      `kd`.

 - `verbose::Bool`: Display informational messages and the full list of
      parameters and timers at the end of execution.  Default value `false`.

# Return values

 - `distances::Array{Float64, 2}`: Matrix to output distances into.
 - `neighbors::Array{Int, 2}`: Matrix to output neighbors into.
 - `output_model::RANNModel`: If specified, the kNN model will be output
      here.

"""
function krann(;
               alpha::Union{Float64, Missing} = missing,
               first_leaf_exact::Union{Bool, Missing} = missing,
               input_model::Union{RANNModel, Missing} = missing,
               k::Union{Int, Missing} = missing,
               leaf_size::Union{Int, Missing} = missing,
               naive::Union{Bool, Missing} = missing,
               query = missing,
               random_basis::Union{Bool, Missing} = missing,
               reference = missing,
               sample_at_leaves::Union{Bool, Missing} = missing,
               seed::Union{Int, Missing} = missing,
               single_mode::Union{Bool, Missing} = missing,
               single_sample_limit::Union{Int, Missing} = missing,
               tau::Union{Float64, Missing} = missing,
               tree_type::Union{String, Missing} = missing,
               verbose::Union{Bool, Missing} = missing,
               points_are_rows::Bool = true)
  # Force the symbols to load.
  ccall((:loadSymbols, krannLibrary), Nothing, ());

  IORestoreSettings("K-Rank-Approximate-Nearest-Neighbors (kRANN)")

  # Process each input argument before calling mlpackMain().
  if !ismissing(alpha)
    IOSetParam("alpha", convert(Float64, alpha))
  end
  if !ismissing(first_leaf_exact)
    IOSetParam("first_leaf_exact", convert(Bool, first_leaf_exact))
  end
  if !ismissing(input_model)
    krann_internal.IOSetParamRANNModel("input_model", convert(RANNModel, input_model))
  end
  if !ismissing(k)
    IOSetParam("k", convert(Int, k))
  end
  if !ismissing(leaf_size)
    IOSetParam("leaf_size", convert(Int, leaf_size))
  end
  if !ismissing(naive)
    IOSetParam("naive", convert(Bool, naive))
  end
  if !ismissing(query)
    IOSetParamMat("query", query, points_are_rows)
  end
  if !ismissing(random_basis)
    IOSetParam("random_basis", convert(Bool, random_basis))
  end
  if !ismissing(reference)
    IOSetParamMat("reference", reference, points_are_rows)
  end
  if !ismissing(sample_at_leaves)
    IOSetParam("sample_at_leaves", convert(Bool, sample_at_leaves))
  end
  if !ismissing(seed)
    IOSetParam("seed", convert(Int, seed))
  end
  if !ismissing(single_mode)
    IOSetParam("single_mode", convert(Bool, single_mode))
  end
  if !ismissing(single_sample_limit)
    IOSetParam("single_sample_limit", convert(Int, single_sample_limit))
  end
  if !ismissing(tau)
    IOSetParam("tau", convert(Float64, tau))
  end
  if !ismissing(tree_type)
    IOSetParam("tree_type", convert(String, tree_type))
  end
  if verbose !== nothing && verbose === true
    IOEnableVerbose()
  else
    IODisableVerbose()
  end

  IOSetPassed("distances")
  IOSetPassed("neighbors")
  IOSetPassed("output_model")
  # Call the program.
  krann_mlpackMain()

  return IOGetParamMat("distances", points_are_rows),
         IOGetParamUMat("neighbors", points_are_rows),
         krann_internal.IOGetParamRANNModel("output_model")
end


text:
using Test
include("../src/IntCodeMachine.jl")
using .IntCodeMachine

@testset "Day 2 Example - Add, Mul, End" begin
    program = [1, 9, 10, 3, 2, 3, 11, 0, 99, 30, 40, 50]
    expected = [3500, 9, 10, 70, 2, 3, 11, 0, 99, 30, 40, 50]
    @test IntCodeMachine.execute!(IntCodeMachine.IntCode(program)).memory == expected
end

@testset "Day 5 Example - Immediate/Position mode" begin
    program = [1002, 4, 3, 4, 33]
    expected = [1002, 4, 3, 4, 99]
    @test IntCodeMachine.execute!(IntCodeMachine.IntCode(program)).memory == expected
end

# "Day 5 Example - I/O" begin
# program = [3, 0, 4, 0, 99]
# IntCodeMachine.execute!(IntCodeMachine.IntCode(program))

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
#=
Created on Monday 13 Jan 2020
Last update: Monday 20 Jan 2020

@author: Michiel Stock
michielfmstock@gmail.com

Just some plots of the testing functions
=#

using STMO, Plots, LaTeXStrings

using STMO.TestFuns: branin, rosenbrock, rastrigine, flower, booth, fquadr, fnonquadr

cmap = :speed

# NOTE: some contours are log-scaled to visualize the minimum better

# branin
contour(-5:0.02:20, -3:0.02:23, (x1, x2) -> log(branin([x1, x2])), xlabel="\$ x_1 \$",
                                    ylabel="\$ x_2 \$", color=cmap)
savefig("figures/brainin.png")

# rosenbrock
contour(-2:0.01:2, -2:0.01:2, (x1, x2) -> log(rosenbrock([x1, x2])), xlabel="\$ x_1 \$",
                                    ylabel="\$ x_2 \$", color=cmap)
savefig("figures/rosenbrock.png")

# rastrigine
contourf(-5.2:0.1:5.2, -5.2:0.1:5.2, (x1, x2) -> rastrigine([x1, x2]),
                    xlabel="\$ x_1 \$", ylabel="\$ x_2 \$", color=cmap)
savefig("figures/rastrigine.png")

# flower
contourf(-3:0.1:3, -3:0.1:3, (x1, x2) -> flower([x1, x2]),
                        xlabel="\$ x_1 \$", ylabel="\$ x_2 \$", color=cmap)
savefig("figures/flower.png")

# booth
contour(-10:0.1:10, -10:0.1:10, (x1, x2) -> booth([x1, x2]) |> log,
                        xlabel="\$ x_1 \$", ylabel="\$ x_2 \$", color=cmap)
savefig("figures/booth.png")

# quadratic
contour(-10:0.1:10, -5:0.1:5, (x1, x2) -> fquadr((x1, x2)), xlabel="\$ x_1 \$",
                ylabel="\$ x_2 \$", title="quadratic", fill=true)
savefig("figures/quadratic.png")

# non-quadratic
contour(-2:0.1:2, -1:0.1:1, (x1, x2) -> fnonquadr((x1, x2)), xlabel="\$ x_1 \$",
                ylabel="\$ x_2 \$", title="non-quadratic", fill=true)
savefig("figures/nonquadratic.png")


text:
#Normalized Difference Vegetation Index
function NDVI(img::Landsat7ETMImage)
    (img.b4 - img.b3) / (img.b4 + img.b3)
end

function NDVI(img::Landsat8OLIImage)
    (img.b5 - img.b4) / (img.b5 + img.b4)
end

function NDVI(img::AbstractSentinelImage)
    (img.b8 - img.b4) / (img.b8 + img.b4);
end

function TassledCap(img::Landsat8OLIImage)
    brightness = img.b1 * 0.3029 + img.b2 * 0.2786 + img.b3 * 0.4733 + img.b4 * 0.5599 + img.b5 * 0.508 + img.b6 * 0.1872
    greenness =  img.b1 * (-0.2941) + img.b2 * (-0.243) + img.b3 * (-0.5424) + img.b4 * 0.7276 + img.b5 * 0.0713 + img.b6 * (-0.1608)
    wetness = img.b1 * 0.1511 + img.b2 * 0.1973 + img.b3 * 0.3283 + img.b4 * 0.3407 + img.b5 * (-0.7117) + img.b6 * (-0.4559)
end

function TassledCap(img::Landsat7ETMImage)
   brightness =  img.b1 * 0.3561 + img.b2*  0.3972 + img.b3 * 0.3904 + img.b4 * 0.6966 + img.b5 * 0.2286 + img.b7 0.1596
   greenness=  img.b1 * -0.3344 + img.b2 * -0.3544 + img.b3 * -0.4556 + img.b4 *  0.6966+-0.0242+-0.2630
   wetness =  img.b1 * 0.2626 + img.b2 * 0.2141 + img.b3 * 0.0926 + img.b4 * 0.0656 + img.b5 * -0.7629 + img.b7 * -0.5388
end

#---TO DO ---
#Enhanced Vegetation Index
#function EVI end

#Modified Soil Adjusted Vegetation Index
#function MSAVI end

#Normalized Burn Ratio Index
#function NBRI end

#Normalized Difference Water Index
#function NDWI end


text:
module SincFun

using SpecialFunctions

import Base: +,-,*,/,\
import SpecialFunctions: sinint, cosint

include("helper.jl")
include("Domains.jl")

export sincfun, hilbert, domain

#
# A sincfun represents a function by
#
# f(x) = ω(ϕ^{-1}(x))/ϕ'(ϕ^{-1}(x))
#       * Σ_{j=-n}^{+n} (-1)^j f(ϕ(jh))*ϕ'(jh)/(ϕ^{-1}(x)-jh)
#       / Σ_{j=-n}^{+n} (-1)^j ω(jh)/(ϕ^{-1}(x)-jh)
#
# where ϕ is a double exponential variable transformation
# and ω is a double exponential weight.
#

type sincfun{D,T}
    n::Integer
    h::T
    fϕv::Vector{T}
    ϕpv::Vector{T}
    ωv::Vector{T}
    ωscale::T
    ωβ::T
    jh::Vector{T}
    domain::D
end

function sincfun{T<:Number}(f::Function,domain::Domain{T})
    n=2^2
    intold,intnew = Inf,-one(T)
    while n < 2^14 && abs(intold-intnew) > abs(intnew)*eps(T)^(1/3)
        n *= 2
        jh = h(T,n)*collect(-n:n)
        sinhv,coshv = sinh.(jh)*π/2,cosh.(jh)*π/2
        fϕv = T[f(ψ(domain,z)) for z in sinhv]
        ϕpv = ψp(domain,sinhv).*coshv
        test = abs.(fϕv.*ϕpv)
        cutoff = .!(.!isinf.(test).*.!isnan.(test))
        fϕv[cutoff],ϕpv[cutoff] = zeros(T,sum(cutoff)),zeros(T,sum(cutoff))
        intold,intnew = intnew,h(T,n)*sum(fϕv.^2.*ϕpv)
    end
    jh = h(T,n)*(collect(-n:n-1)+one(T)/2)
    sinhv,coshv = sinh.(jh)*π/2,cosh.(jh)*π/2
    fϕv = interlace([fϕv;T[f(ψ(domain,z)) for z in sinhv]])
    ϕpv = interlace([ϕpv;ψp(domain,sinhv).*coshv])
    test = abs.(fϕv.*ϕpv)
    cutoff = .!(.!isinf.(test).*.!isnan.(test))
    fϕv[cutoff],ϕpv[cutoff] = zeros(T,sum(cutoff)),zeros(T,sum(cutoff))
    ωscale = maximum(test)
    tru = div(findfirst(reverse(interlace2(h(T,n)/2*fϕv.^2.*ϕpv)) .> ωscale*eps(T)^3),2)
    fϕv,ϕpv,jh = fϕv[tru+1:4n-tru+1],ϕpv[tru+1:4n-tru+1],h(T,n)/2*collect(-2n+tru:2n-tru)
    ωβ = -2log(eps(T))/(cosh(jh[end])-one(T))
    ωv = ωscale*ω.(ωβ,jh)
    return sincfun{typeof(domain),T}(length(fϕv),h(T,n)/2,fϕv,ϕpv,ωv,ωscale,ωβ,jh,domain)
end
sincfun(f::Function) = sincfun(f,Finite())

Base.length(sf::sincfun) = sf.n
Base.eltype{D,T}(sf::sincfun{D,T}) = T
domain(sf::sincfun) = sf.domain

## Evaluation by the barycentric formula

function Base.getindex{D<:Domain,S<:Number,T<:Number}(sf::sincfun{D,S},x::T)
    xc = convert(promote_type(S,T),x)
    z = ψinv(sf.domain,xc)
    barycentric(sf,xc)*singularities(sf.domain,z)
end
Base.getindex{D<:Domain,S<:Number,T<:Number}(sf::sincfun{D,S},x::AbstractVector{T}) = promote_type(S,T)[sf[xk] for xk in x]
Base.getindex{D<:Domain,S<:Number,T<:Number,N}(sf::sincfun{D,S},x::AbstractArray{T,N}) = reshape(sf[vec(x)],size(x))

function barycentric{D<:Domain,T<:Number}(sf::sincfun{D,T},x::T)
    t = asinh(2ψinv(sf.domain,x)/π)
    if t < sf.jh[1]
        return envelope(sf,sf.jh[1])*sf.fϕv[1]*sf.ϕpv[1]/sf.ωv[1]
    elseif t > sf.jh[end]
        return envelope(sf,sf.jh[end])*sf.fϕv[end]*sf.ϕpv[end]/sf.ωv[end]
    else
        idx = findfirst(sf.jh,t)
        if idx == 0
            common = t-sf.jh
            num,den = sf.fϕv.*sf.ϕpv./common,sf.ωv./common
            valN,valD = sum(num[1:2:end])-sum(num[2:2:end]),sum(den[1:2:end])-sum(den[2:2:end])
            return envelope(sf,t)*valN/valD
        else
            return envelope(sf,t)*sf.fϕv[idx]*sf.ϕpv[idx]/sf.ωv[idx]
        end
    end
end
envelope{D<:Domain,T<:Number}(sf::sincfun{D,T},t::T) = sf.ωscale*ω(sf.ωβ,t)/ψp(sf.domain,sinh(t)*π/2)/(cosh(t)*π/2)

# Algebra

for op in (:+,:-,:*)
    @eval begin
        function $op{D<:Domain,T<:Number}(c::Number,sf::sincfun{D,T})
            sf1 = deepcopy(sf)
            sf1.fϕv = $op(convert(T,c),sf.fϕv)
            return sf1
        end
        $op{D<:Domain,T<:Number}(sf::sincfun{D,T},c::Number) = $op(c,sf)

        function $op{D<:Domain,T<:Number}(sf1::sincfun{D,T},sf2::sincfun{D,T})
            #TODO: sf.domain = $op(sf1.domain,sf2.domain)
            sincfun(x->$op(sf1[x],sf2[x]),sf1.domain)
        end
    end
end

# real, imag, conj, abs2

for op in (:(Base.real),:(Base.imag),:(Base.conj),:(Base.abs2))
    @eval begin
        $op{D<:Domain,T<:Real}(sf::sincfun{D,T}) = deepcopy(sf)
    end
end

for op in (:(Base.real),:(Base.imag))
    @eval begin
        function $op{D<:Domain,T<:Real}(sf::sincfun{D,Complex{T}})
            sincfun{typeof(sf.domain),T}(sf.n,real(sf.h),$op(sf.fϕv),real(sf.ϕpv),real(sf.ωv),real(sf.ωscale),real(sf.ωβ),real(sf.jh),sf.domain)
        end
    end
end

function Base.conj{D<:Domain,T<:Real}(sf::sincfun{D,Complex{T}})
    sf1 = deepcopy(sf)
    sf1.fϕv = conj(sf1.fϕv)
    return sf1
end

function Base.abs2{D<:Domain,T<:Real}(sf::sincfun{D,Complex{T}})
    rsf,isf = real(sf),imag(sf)
    return rsf*rsf + isf*isf
end

# sum, dot, cumsum, norm, and diff.

function Base.sum{D<:Domain,T<:Number}(sf::sincfun{D,T})
    sinhv = sinh.(sf.jh)*π/2
    singv = singularities(sf.domain,sinhv)
    sf.h*sum(sf.fϕv.*sf.ϕpv.*singv)
end
Base.dot{D<:Domain,T<:Number}(sf1::sincfun{D,T},sf2::sincfun{D,T}) = sum(conj(sf1)*sf2)

function Base.cumsum{D<:Domain,T<:Number}(sf::sincfun{D,T})
    SM = sinc.(-1,one(T)*collect(-sf.n+1:sf.n-1))
    sf1 = deepcopy(sf)
    temp = sf.h*sf.fϕv.*sf.ϕpv
    [sf1.fϕv[i] = sum(SM[i+sf.n-1:-1:i].*temp) for i=1:sf.n]
    return sf1
end

function Base.norm{D<:Domain,T<:Number}(sf::sincfun{D,T},p::Int)
    if isodd(p) error("Only even integer norms supported.") end
    sinhv = sinh.(sf.jh)*π/2
    singv = singularities(sf.domain,sinhv)
    abs(sf.h*sum((sf.fϕv.*singv).^p.*sf.ϕpv))^(one(T)/p)
end
Base.norm{D<:Domain,T<:Number}(sf::sincfun{D,T}) = norm(sf,2)

function hilbert{D<:Domain,T<:Number}(sf::sincfun{D,T})
    sinhv = sinh.(sf.jh)*π/2
    SM = (sinc.(0,one(T)/2/sf.h*(sf.jh.-sf.jh')).*(sf.jh.-sf.jh')).^2./(ψ(sf.domain,sinhv).-ψ(sf.domain,sinhv)')
    [SM[i,i] = zero(T) for i=1:sf.n]
    sf1 = deepcopy(sf)
    singv = singularities(sf.domain,sinhv)
    temp = -sf.fϕv.*sf.ϕpv.*singv*π/2/sf.h
    sf1.fϕv = SM*temp
    return sf1
end

#=
function Base.diff{D<:Domain,T<:Number}(sf::sincfun{D,T})
    SM = sinc(1,one(T)*[-sf.n+1:sf.n-1])
    #SM = sinc(1,one(T)*[-(sf.n-1)/2:(sf.n-1)/2].-[-(sf.n-1)/2:(sf.n-1)/2]')
    sf1 = deepcopy(sf)
    temp = sf.fϕv./(sf.h*sf.ϕpv)
    [sf1.fϕv[i] = sum(SM[i+sf.n-1:-1:i].*temp) for i=1:sf.n]
    sf1.fϕv .*= sf.ϕpv
    #sf1.fϕv = SM*temp
    return sf1
end
=#
#=
function Base.diff{D<:Domain,T<:Number}(sf::sincfun{D,T})
    SM = sinc(1,one(T)*[-sf.n:sf.n])
    sinhv,coshv = convert(T,π)/2*sinh(sf.jh),convert(T,π)/2*cosh(sf.jh)
    temp = ψp(sf.domain,sinhv).*coshv
    phippoverphip = sinhv./coshv - 2tanh(sinhv).*coshv
    sf1 = deepcopy(sf)
    [sf1.fϕpv[i] = sum(SM[i:i+sf.n-1].*(sf.fϕpv/sf.h)) for i=1:sf.n]
    sf1.fϕpv = (sf1.fϕpv-sf.fϕpv.*phippoverphip)#./temp
    #sf1.fϕv = sf1.fϕv./(sf.ϕpv*sf.h)
    #sf1.ϕpv = ones(T,length(sf.ϕpv))
    return sf1
end
=#

include("roots.jl")
include("fftBigFloat.jl")
include("KrylovMethods.jl")
include("steig.jl")
include("sinc.jl")
include("ConformalMap.jl")
include("SincMatrix.jl")
include("BandMatrix.jl")

end #module


text:
using Pretty
using Test

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


text:
mutable struct CTPF <: TopicModel
	"Collaborative topic Poisson factorization."

	K::Int
	M::Int
	V::Int
	U::Int
	N::Vector{Int}
	C::Vector{Int}
	R::Vector{Int}
	corp::Corpus
	topics::VectorList{Int}
	scores::Matrix{Float64}
	libs::VectorList{Int}
	drecs::VectorList{Int}
	urecs::VectorList{Int}
	a::Float64
	b::Float64
	c::Float64
	d::Float64
	e::Float64
	f::Float64
	g::Float64
	h::Float64
	alef::Matrix{Float64}
	alef_old::Matrix{Float64}
	alef_temp::Matrix{Float64}
	he::Matrix{Float64}
	he_old::Matrix{Float64}
	he_temp::Matrix{Float64}
	bet::Vector{Float64}
	bet_old::Vector{Float64}
	vav::Vector{Float64}
	vav_old::Vector{Float64}
	gimel::VectorList{Float64}
	gimel_old::VectorList{Float64}
	zayin::VectorList{Float64}
	zayin_old::VectorList{Float64}
	dalet::Vector{Float64}
	dalet_old::Vector{Float64}
	het::Vector{Float64}
	het_old::Vector{Float64}
	phi::MatrixList{Float64}
	xi::MatrixList{Float64}
	elbo::Float64

	function CTPF(corp::Corpus, K::Integer)
		check_corp(corp)
		K > 0 || throw(ArgumentError("Number of topics must be a positive integer."))

		M, V, U = size(corp)
		N = [length(doc) for doc in corp]
		C = [size(doc) for doc in corp]
		R = [length(doc.readers) for doc in corp]

		topics = [collect(1:V) for _ in 1:K]
		scores = zeros(M, U)

		libs = [Int[] for _ in 1:U]
		for d in 1:M, u in corp[d].readers
			push!(libs[u], d)
		end

		urecs = VectorList{Int}(undef, U)
		for u in 1:U
			ur = trues(M)
			ur[libs[u]] .= false
			urecs[u] = findall(ur)
		end

		drecs = VectorList{Int}(undef, M)
		for d in 1:M
			nr = trues(U)
			nr[corp[d].readers] .= false
			drecs[d] = findall(nr)
		end

		a, b, c, d, e, f, g, h = fill(0.1, 8)

		alef = exp.(rand(Dirichlet(V, 1.0), K)' .- 0.5)
		alef_old = copy(alef)
		alef_temp = fill(a, K, V)
		he = ones(K, U)
		he_old = copy(he)
		he_temp = fill(e, K, U)
		bet = ones(K)
		bet_old = copy(bet)
		vav = ones(K)
		vav_old = copy(vav)
		gimel = [ones(K) for _ in 1:M]
		gimel_old = deepcopy(gimel)
		zayin = [ones(K) for _ in 1:M]
		zayin_old = deepcopy(zayin)
		dalet = ones(K)
		dalet_old = copy(dalet)
		het = ones(K)
		het_old = copy(het)
		phi = [ones(K, N[d]) / K for d in 1:min(M, 1)]
		xi = [ones(2K, R[d]) / 2K for d in 1:min(M, 1)]
		elbo = 0

		model = new(K, M, V, U, N, C, R, copy(corp), topics, scores, libs, drecs, urecs, a, b, c, d, e, f, g, h, alef, alef_old, alef_temp, he, he_old, he_temp, bet, bet_old, vav, vav_old, gimel, gimel_old, zayin, zayin_old, dalet, dalet_old, het, het_old, phi, xi, elbo)
		return model
	end
end

function Elogpya(model::CTPF, d::Int)
	"Compute E_q[log(P(ya))]."

	x = -dot(model.gimel[d] ./ (model.dalet .* model.vav), sum(model.he, dims=2))
	readers, ratings = model.corp[d].readers, model.corp[d].ratings
	for (u, (re, ra)) in enumerate(zip(readers, ratings)), i in 1:model.K
		binom = Binomial(ra, model.xi[1][i,u])
		x += (ra * model.xi[1][i,u] * (digamma(model.gimel[d][i]) - log(model.dalet[i]) + digamma(model.he[i,re]) - log(model.vav[i])) - sum([pdf(binom, y) * loggamma(y + 1) for y in 0:ra]))
	end
	return x
end

function Elogpyb(model::CTPF, d::Int)
	"Compute E_q[log(P(yb))]."

	x = -dot(model.zayin[d] ./ (model.het .* model.vav), sum(model.he, dims=2))
	readers, ratings = model.corp[d].readers, model.corp[d].ratings
	for (u, (re, ra)) in enumerate(zip(readers, ratings)), i in 1:model.K
		binom = Binomial(ra, model.xi[1][model.K+i,u])
		x += (ra * model.xi[1][model.K+i,u] * (digamma(model.zayin[d][i]) - log(model.het[i]) + digamma(model.he[i,re]) - log(model.vav[i])) - sum([pdf(binom, y) * loggamma(y + 1) for y in 0:ra]))
	end
	return x
end

function Elogpz(model::CTPF, d::Int)
	"Compute E_q[log(P(z))]."

	x = -dot(model.gimel[d] ./ (model.dalet .* model.bet), sum(model.alef, dims=2))
	terms, counts = model.corp[d].terms, model.corp[d].counts
	for (n, (j, c)) in enumerate(zip(terms, counts)), i in 1:model.K
		binom = Binomial(c, model.phi[1][i,n])
		x += (c * model.phi[1][i,n] * (digamma(model.gimel[d][i]) - log(model.dalet[i]) + digamma(model.alef[i,j]) - log(model.bet[i])) - sum([pdf(binom, z) * loggamma(z + 1) for z in 0:c]))
	end
	return x
end

function Elogpbeta(model::CTPF)
	"Compute E_q[log(P(beta))]."

	x = model.V * model.K * (model.a * log(model.b) - loggamma(model.a))
	for j in 1:model.V, i in 1:model.K
		x += (model.a - 1) * (digamma(model.alef[i,j]) - log(model.bet[i])) - model.b * model.alef[i,j] / model.bet[i]
	end
	return x
end

function Elogptheta(model::CTPF, d::Int)
	"Compute E_q[log(P(theta))]."

	x = model.K * (model.c * log(model.d) - loggamma(model.c))
	for i in 1:model.K
		x += (model.c - 1) * (digamma(model.gimel[d][i]) - log(model.dalet[i])) - model.d * model.gimel[d][i] / model.dalet[i]
	end
	return x
end

function Elogpeta(model::CTPF)
	"Compute E_q[log(P(eta))]."

	x = model.U * model.K * (model.e * log(model.f) - loggamma(model.e))
	for u in 1:model.U, i in 1:model.K
		x += (model.e - 1) * (digamma(model.he[i,u]) - log(model.vav[i])) - model.f * model.he[i,u] / model.vav[i]
	end
	return x
end

function Elogpepsilon(model::CTPF, d::Int)
	"Compute E_q[log(P(epsilon))]."

	x = model.K * (model.g * log(model.h) - loggamma(model.g))
	for i in 1:model.K
		x += (model.g - 1) * (digamma(model.zayin[d][i]) - log(model.het[i])) - model.h * model.zayin[d][i] / model.het[i]
	end
	return x
end

function Elogqy(model::CTPF, d::Int)
	"Compute E_q[log(q(y))]."

	x = 0
	for (u, ra) in enumerate(model.corp[d].ratings)
		x -= entropy(Multinomial(ra, model.xi[1][:,u]))
	end
	return x
end

function Elogqz(model::CTPF, d::Int)
	"Compute E_q[log(q(z))]."

	x = 0
	for (n, c) in enumerate(model.corp[d].counts)
		x -= entropy(Multinomial(c, model.phi[1][:,n]))
	end
	return x
end

function Elogqbeta(model::CTPF)
	"Compute E_q[log(q(beta))]."

	x = 0
	for j in 1:model.V, i in 1:model.K
		x -= entropy(Gamma(model.alef[i,j], 1 / model.bet[i]))
	end
	return x
end

function Elogqtheta(model::CTPF, d::Int)
	"Compute E_q[log(q(theta))]."

	x = 0
	for i in 1:model.K
		x -= entropy(Gamma(model.gimel[d][i], 1 / model.dalet[i]))
	end
	return x
end

function Elogqeta(model::CTPF)
	"Compute E_q[log(q(eta))]."

	x = 0
	for u in 1:model.U, i in 1:model.K
		x -= entropy(Gamma(model.he[i,u], 1 / model.vav[i]))
	end
	return x
end

function Elogqepsilon(model::CTPF, d::Int)
	"Compute E_q[log(q(epsilon))]."

	x = 0
	for i in 1:model.K
		x -= entropy(Gamma(model.zayin[d][i], 1 / model.het[i]))
	end
	return x
end

function update_elbo!(model::CTPF)
	"Update the evidence lower bound."

	model.elbo = Elogpbeta(model) + Elogpeta(model) - Elogqbeta(model) - Elogqeta(model)
	for d in 1:model.M
		terms = model.corp[d].terms
		readers = model.corp[d].readers

		model.phi[1] = additive_logistic(digamma.(model.gimel_old[d]) - log.(model.dalet_old) - log.(model.bet_old) .+ digamma.(model.alef_old[:,terms]), dims=1)
		model.xi[1] = additive_logistic(vcat(digamma.(model.gimel_old[d]) - log.(model.dalet_old) - log.(model.vav_old) .+ digamma.(model.he_old[:,readers]), digamma.(model.zayin_old[d]) - log.(model.het_old) - log.(model.vav_old) .+ digamma.(model.he_old[:,readers])), dims=1)

		model.elbo += Elogpya(model, d) + Elogpyb(model, d) + Elogpz(model, d) + Elogptheta(model, d) + Elogpepsilon(model, d) - Elogqy(model, d) - Elogqz(model, d) - Elogqtheta(model, d) - Elogqepsilon(model, d)
	end

	return model.elbo
end

function update_alef!(model::CTPF)
	"Update alef."
	"Analytic."

	model.alef_old = model.alef
	model.alef = model.alef_temp
	model.alef_temp = fill(model.a, model.K, model.V)
end

function update_alef!(model::CTPF, d::Int)
	"Update alef."
	"Analytic."

	terms, counts = model.corp[d].terms, model.corp[d].counts
	model.alef_temp[:,terms] += model.phi[1] .* counts'
end

function update_he!(model::CTPF)
	"Update he."
	"Analytic."

	model.he_old = model.he
	model.he = model.he_temp
	model.he_temp = fill(model.e, model.K, model.U)
end

function update_he!(model::CTPF, d::Int)
	"Update he."
	"Analytic."

	readers, ratings = model.corp[d].readers, model.corp[d].ratings
	model.he_temp[:,readers] += (model.xi[1][1:model.K,:] + model.xi[1][model.K+1:end,:]) .* ratings'
end

function update_bet!(model::CTPF)
	"Update bet."
	"Analytic"

	model.bet_old = model.bet
	model.bet = model.b .+ sum(model.gimel) ./ model.dalet
end

function update_vav!(model::CTPF)
	"Update vav."
	"Analytic."

	model.vav_old = model.vav
	model.vav = model.f .+ sum(model.gimel) ./ model.dalet + sum(model.zayin) ./ model.het
end

function update_dalet!(model::CTPF)
	"Update dalet."
	"Analytic."

	model.dalet_old = model.dalet
	model.dalet = model.d .+ vec(sum(model.alef, dims=2)) ./ model.bet + vec(sum(model.he, dims=2)) ./ model.vav
end

function update_het!(model::CTPF)
	"Update het."
	"Analytic"

	model.het_old = model.het
	model.het = model.h .+ vec(sum(model.he, dims=2)) ./ model.vav
end

function update_gimel!(model::CTPF, d::Int)
	"Update gimel."
	"Analytic."

	model.gimel_old[d] = model.gimel[d]

	counts, ratings = model.corp[d].counts, model.corp[d].ratings
	model.gimel[d] = model.c .+ model.phi[1] * counts + model.xi[1][1:model.K,:] * ratings
end

function update_zayin!(model::CTPF, d::Int)
	"Update zayin."
	"Analytic"

	model.zayin_old[d] = model.zayin[d]

	ratings = model.corp[d].ratings
	model.zayin[d] = model.g .+ model.xi[1][model.K+1:end,:] * ratings
end

function update_phi!(model::CTPF, d::Int)
	"Update phi."
	"Analytic."

	terms = model.corp[d].terms
	model.phi[1] = additive_logistic(digamma.(model.gimel[d]) - log.(model.dalet) - log.(model.bet) .+ digamma.(model.alef[:,terms]), dims=1)
end

function update_xi!(model::CTPF, d::Int)
	"Update xi."
	"Analytic."

	readers = model.corp[d].readers
	model.xi[1] = additive_logistic(vcat(digamma.(model.gimel[d]) - log.(model.dalet) - log.(model.vav) .+ digamma.(model.he[:,readers]), digamma.(model.zayin[d]) - log.(model.het) - log.(model.vav) .+ digamma.(model.he[:,readers])), dims=1)
end

function train!(model::CTPF; iter::Integer=150, tol::Real=1.0, viter::Integer=10, vtol::Real=1/model.K^2, checkelbo::Real=1, printelbo::Bool=true)
	"Coordinate ascent optimization procedure for collaborative topic Poisson factorization variational Bayes algorithm."

	check_model(model)
	all([tol, vtol] .>= 0)												|| throw(ArgumentError("Tolerance parameters must be nonnegative."))
	all([iter, viter] .>= 0)											|| throw(ArgumentError("Iteration parameters must be nonnegative."))
	(isa(checkelbo, Integer) & (checkelbo > 0)) | (checkelbo == Inf)	|| throw(ArgumentError("checkelbo parameter must be a positive integer or Inf."))
	all([isempty(doc) for doc in model.corp]) && (iter = 0)
	(checkelbo <= iter) && update_elbo!(model)

	for k in 1:iter
		for d in 1:model.M
			for _ in 1:viter
				update_xi!(model, d)
				update_phi!(model, d)
				update_zayin!(model, d)
				update_gimel!(model, d)
				if norm(model.gimel[d] - model.gimel_old[d]) < vtol
					break
				end
			end
			update_he!(model, d)
			update_alef!(model, d)
		end
		update_he!(model)
		update_alef!(model)
		update_dalet!(model)
		update_het!(model)
		update_bet!(model)
		update_vav!(model)

		if check_elbo!(model, checkelbo, printelbo, k, tol)
			break
		end
	end

	Ebeta = model.alef ./ model.bet
	model.topics = [reverse(sortperm(vec(Ebeta[i,:]))) for i in 1:model.K]

	Eeta = model.he ./ model.vav
	for d in 1:model.M
		Etheta = model.gimel[d] ./ model.dalet
		Eepsilon = model.zayin[d] ./ model.het
		model.scores[d,:] = sum(Eeta .* (Etheta + Eepsilon), dims=1)
	end

	model.urecs = VectorList{Int}(undef, model.U)
	for u in 1:model.U
		ur = trues(model.M)
    	ur[model.libs[u]] .= false
		model.urecs[u] = findall(ur)[reverse(sortperm(model.scores[ur,u]))]
	end

	model.drecs = VectorList{Int}(undef, model.M)
	for d in 1:model.M
    	nr = trues(model.U)
    	nr[model.corp[d].readers] .= false
    	model.drecs[d] = findall(nr)[reverse(sortperm(model.scores[d,nr]))]
	end
	nothing
end


text:
function test_text_hist(parallelization::HistogramParallelization)
    @testset "1D text" begin
        strings = split(lowercase(read(joinpath(@__DIR__, "pride_and_prejudice.txt"), String)))
        expected = Dict([s => count(isequal(s), strings) for s in unique(strings)])

        h = create_fast_histogram(UnboundedWidth(), HashFunction(), parallelization)

        increment_bins!(h, strings)
        @test counts(h) == expected

        zero!(h)
        increment_bins!(h, strings)
        @test counts(h) == expected

        zero!(h)
        @test isempty(counts(h))
    end
end


text:
module AcousticNewAssemble
using FEniCS

c = 5000
#problem variables
dt = 0.000004;
global t = 0;
T = 0.004;

mesh = RectangleMesh(Point([-2., -2.]),Point([2., 2.]),40,40)
V = FunctionSpace(mesh,"Lagrange",1)

# Previous and current solution
u1= interpolate(Constant(0.0), V)
u0= interpolate(Constant(0.0), V)

# Variational problem at each time
u = TrialFunction(V)
v = TestFunction(V)

a = u*v*dx + dt*dt*c*c*inner(grad(u), grad(v))*dx
L = 2*u1*v*dx-u0*v*dx

delta = Expression("sin(c*10*t)",degree=2,c=c,t=t)

#Define boundary conditions
top ="on_boundary && near(x[1], 2)"
BCT = DirichletBC(V,Constant(0.0),top)

down ="on_boundary && near(x[1], -2)"
BCD = DirichletBC(V,Constant(0.0),down)

left = "on_boundary && near(x[0],-2)"
BCL = DirichletBC(V,delta,left)

right = "on_boundary && near(x[0],2)"
BCR = DirichletBC(V,Constant(0.0),right)

bcs_dir = [BCL,BCD,BCT,BCR]
bcs_neu = [BCL,BCD,BCT]

u=FeFunction(V)

while t <= T
    delta.pyobject.t = t
    lvsolve(a,L,u,bcs_dir) #linear variational solver
    assign(u0,u1)
    assign(u1,u)
    global t +=dt
end
end#module


text:
# using Atom
# include("types.jl")
using FileIO
using Images
using ImageMagick # for load rgb
# using Images
"""
ids start from 0
"""
mutable struct KlgReader <: RgbdReader
    frame_num::Int64
    ids::Vector{Int64}
    timestamps::Vector{Int64}
    dw::Int64
    dh::Int64
    num_pixels::Int64
    f_idx::Int64
    klg_file_name::String
    klg_file::IOStream
    # store the position of each frame
    frame_position::Dict{Int64, Int64}

    depth_buffer::Array{UInt8}
    rgb_buffer::Array{UInt8}

    depth_uncompressed::Array{UInt16}
    #codec::ZlibDecompressor
    cam::CameraModel.RgbdCamParams
end

function destruct_klg!(klg_data::KlgReader)
    Base.close(klg_data.klg_file)
end

function KlgReader(
    klg_file_name::String;
    cam::CamT
) where CamT <: CameraModel.AbstractCamParams
    dw = cam.width
    dh = cam.height
    num_pixels = dw * dh
    klg_file::IOStream = open(klg_file_name,"r")
    FrameNum::Int64 = Base.read(klg_file, Int32)
    frame_position = Dict{Int64, Int64}()
    # FrameNum = Base.read(klg_file, Int32)
    # frame_position[0] = 4

    depth_buffer::Array{UInt8} = Array{UInt8}(undef, num_pixels * 2);
    rgb_buffer::Array{UInt8} = Array{UInt8}(undef, num_pixels * 3);

    depth_uncompressed =  Array{UInt8}(undef, num_pixels);

    #codec = ZlibDecompressor()
    #TranscodingStreams.initialize(codec)
    # finalizer(k, k->print("finalized!"))

    ids = collect(0:FrameNum-1)
    timestamps = ids .+ 1

    klg_data = KlgReader(
        FrameNum, ids, timestamps, dw, dh, num_pixels, -1,
        klg_file_name, klg_file,
        frame_position, depth_buffer, rgb_buffer, depth_uncompressed,
        cam)
    finalizer(destruct_klg!, klg_data)
    return klg_data
end

"""
`isrgb`: true for rgb order, false for bgr order
"""
function readRGBD_base(klg_data::KlgReader; isrgb::Bool=true)
    time_stamp = Base.read(klg_data.klg_file, Int64)
    depth_size = Base.read(klg_data.klg_file, Int32)
    rgb_size = Base.read(klg_data.klg_file, Int32)

#     println(time_stamp)
    # @show depth_size
    # @show rgb_size

    readbytes!(klg_data.klg_file, klg_data.depth_buffer, depth_size)
    readbytes!(klg_data.klg_file, klg_data.rgb_buffer, rgb_size)

    # # C++ function decompress
    # decompress(klg_data.depth_buffer[1:depth_size], klg_data.depth_uncompressed)

    if klg_data.num_pixels * 2 != depth_size
        # use ccall to perform decompression
        decompLength::UInt64 = size(klg_data.depth_uncompressed, 1) * 2;
        depth_buffer_ptr = Base.unsafe_convert(Ptr{Cuchar}, klg_data.depth_buffer)
        depth_uncompressed_ptr = Base.unsafe_convert(Ptr{Cuchar}, klg_data.depth_uncompressed)
        # //typedef unsigned char Byte;
        # ZEXTERN int ZEXPORT uncompress OF((Bytef *dest,   uLongf *destLen,
        #                                    const Bytef *source, uLong sourceLen));
        res = ccall((:uncompress, :libz),
                    Cint,
                    (Ptr{Cuchar}, Ptr{Culonglong}, Ptr{Cuchar}, Culonglong),
                    depth_uncompressed_ptr, Ref(decompLength), depth_buffer_ptr, depth_size)

        depth = permutedims(reshape(klg_data.depth_uncompressed, (klg_data.dw, klg_data.dh)), [2, 1])
    else
        depth = permutedims(reshape(reinterpret(UInt16, klg_data.depth_buffer), (klg_data.dw, klg_data.dh)), [2, 1])
    end

    if klg_data.num_pixels * 3 != rgb_size
        img = readblob(klg_data.rgb_buffer[1:rgb_size])
    else
        rgb_mat = permutedims(reshape(klg_data.rgb_buffer[1:rgb_size],
                                     (3, klg_data.dw, klg_data.dh)), [1, 3, 2])
        if isrgb
            img = RGB{Normed{UInt8,8}}.(rgb_mat[1, :, :]./255, rgb_mat[2, :, :]./255, rgb_mat[3, :, :]./255)
        else
            img = RGB{Normed{UInt8,8}}.(rgb_mat[3, :, :]./255, rgb_mat[2, :, :]./255, rgb_mat[1, :, :]./255)
        end
            # img = copy(dropdims(reinterpret(RGB{Normed{UInt8,8}}, rgb_mat), dims=1))
    end
    # @show size(img)
    return  img, depth;#,depth_uncompressed;#
end

# end

# klg_reader = io.KlgReader("/home/gzhang8/project/cavedrone/data/belize/2017-06-07.01.klg", 640, 480, 640*480);
# rgb, depth = io.readOneFrame(klg_reader);
# end

function readRGBDFrame(klg_reader::KlgReader; id::Int64, isrgb::Bool=true)
    rgb, depth = readRGBD(klg_reader, id=id, isrgb=isrgb);
    return RGBDFrame(rgb, depth)

end

"""
read data for frame id
id is the index for klg_data, starting from 0
    `isrgb::Bool` true for rgb order, false for bgr order
"""
function readRGBD(klg_data::KlgReader; id::Int64, isrgb::Bool=true)
    @assert(klg_data.frame_num > id)
    old_position = mark(klg_data.klg_file)
    # FrameNum = Base.read(klg_file, Int32)
    current_fid::Int64 = 0
    current_pos::Int64 = 4
    if length(klg_data.frame_position) > 0
        max_fid = maximum(keys(klg_data.frame_position))
        if max_fid >= id
            current_fid = id
            current_pos = klg_data.frame_position[id]
        else
            current_fid = max_fid
            current_pos = klg_data.frame_position[max_fid]
        end

    end
    seek(klg_data.klg_file, current_pos)

    # fast forward to the frame we want to read
    while current_fid < id
        # time_stamp = Base.read(klg_data.klg_file, Int64)
        skip(klg_data.klg_file, 8)
        klg_data.frame_position[current_fid] = current_pos
        # 4 bytes for depth_size
        depth_size = Base.read(klg_data.klg_file, Int32)
        # 4 bytes for rgb_size
        rgb_size = Base.read(klg_data.klg_file, Int32)
        skip(klg_data.klg_file, depth_size + rgb_size)
        current_fid += 1
        current_pos += (16 + depth_size + rgb_size)

    end
    # add the last one in the loop
    klg_data.frame_position[current_fid] = current_pos
    img, depth = readRGBD_base(klg_data, isrgb=isrgb)
    # reset file location
    reset(klg_data.klg_file)
    unmark(klg_data.klg_file)
    return img, depth
end

function view_depth(depth)
    return depth_view = Gray.(depth/4000);
end


text:
using TransformVariables, LogDensityProblems, DynamicHMC, MCMCDiagnostics, Parameters,
    Distributions, Statistics, StatsFuns, ForwardDiff

struct NUTS_result{T}
    chain::Vector{NUTS_Transition{Vector{Float64},Float64}}
    transformation
    samples::Vector{T}
    tuning
end

Base.show(io::IO, n::NUTS_result) = begin
    println(io, "NUTS_result with samples:")
    println(IOContext(io, :limit => true, :compact => true), n.samples)
end

export nuts

function nuts(m :: Model; kwargs...)
    ℓ = makeLogdensity(m)

    result = NUTS_result{}
    t = xform(m; data...)
    P = TransformedLogDensity(t, ℓ)
    ∇P = ADgradient(:ForwardDiff, P)
    chain, tuning = NUTS_init_tune_mcmc(∇P, 1000);
    samples = TransformVariables.transform.(Ref(parent(∇P).transformation), get_position.(chain));
    NUTS_result(chain, t, samples, tuning)
end

function nuts(m :: Model, data :: NamedTuple; kwargs...)
    f1 = makeLogdensity(m)
    ℓ(pars) = f1(merge(data,pars))

    result = NUTS_result{}
    t = xform(m; data...)
    P = TransformedLogDensity(t, ℓ)
    ∇P = ADgradient(:ForwardDiff, P)
    chain, tuning = NUTS_init_tune_mcmc(∇P, 1000);
    samples = TransformVariables.transform.(Ref(parent(∇P).transformation), get_position.(chain));
    NUTS_result(chain, t, samples, tuning)
end


text:
## Code for RENDEREING algorithms #############################################

## On-Off render algorithm ####################################################
struct OnOff_Renderer <: Renderer
    world::World
    on_color::RGB
    off_color::RGB
    OnOff_Renderer(world::World, on_color::RGB=RGB(1.f0,1.f0,1.f0), off_color::RGB=RGB(0.f0,0.f0,0.f0) ) = new(world, on_color, off_color)
end

function OnOff(ray::Ray, renderer::OnOff_Renderer)
    if ray_intersection(renderer.world, ray) === nothing
        return renderer.off_color
    else
        return renderer.on_color
    end
end

## Flat render algorithm ######################################################
struct Flat_Renderer <: Renderer
    world::World
    background_color::RGB

    Flat_Renderer(world::World, background_color::RGB=RGB(0.f0,0.f0,0.f0) ) = new(world, background_color)
end

function Flat(ray::Ray, renderer::Flat_Renderer)
    hit_record = ray_intersection(renderer.world, ray)
    if hit_record === nothing
        return renderer.background_color
    else
        material = hit_record.shape.material
        return get_color(material.brdf.pigment, hit_record.surface_point) + get_color(material.emitted_radiance, hit_record.surface_point)
    end
end

## Path tracer render algorithm ###############################################
struct PathTracer_Renderer <: Renderer
    world::World
    background_color::RGB
    pcg::PCG
    num_of_rays::Int
    max_depth::Int
    russian_roulette_limit::Int

    PathTracer_Renderer(world::World;
                        background_color::RGB = RGB(0.f0,0.f0,0.f0),
                        pcg::PCG = PCG(),
                        num_of_rays::Int = 10,
                        max_depth::Int = 2,
                        russian_roulette_limit::Int = 3 ) =
                        new(world, background_color, pcg, num_of_rays, max_depth, russian_roulette_limit)
end

function PathTracer(ray::Ray, rend::PathTracer_Renderer)
    if ray.depth > rend.max_depth
        return RGB(0.f0, 0.f0, 0.f0)
    end

    hit_record = ray_intersection(rend.world, ray)
    if hit_record === nothing
        return rend.background_color
    end

    hit_material = hit_record.shape.material
    hit_color = get_color(hit_material.brdf.pigment, hit_record.surface_point)
    if hit_material.is_light_source == true
        emitted_radiance = get_color(hit_material.emitted_radiance, hit_record.surface_point) * hit_material.emission_intensity
    else
        emitted_radiance = RGB(0.f0, 0.f0, 0.f0)
    end
    hit_color_lum = max(hit_color.r, hit_color.g, hit_color.b)

    # Russian roulette
    if ray.depth >= rend.russian_roulette_limit
        if pcg_randf(rend.pcg) > hit_color_lum
            # Keep the recursion going, but compensate for other potentially discarded rays
            hit_color *= Float32(1. / (1. - hit_color_lum))
        else
            # Terminate prematurely
            return emitted_radiance
        end
    end

    cum_radiance = RGB(0.f0, 0.f0, 0.f0)
    if hit_color_lum > 0.f0  # Only do costly recursions if it's worth it
        for ray_index ∈ 1:rend.num_of_rays
            new_ray = scatter_ray(hit_material.brdf, rend.pcg, hit_record.ray.dir, hit_record.world_point, hit_record.normal, ray.depth + 1)
            # Recursive call
            new_radiance = PathTracer(new_ray, rend)
            cum_radiance += hit_color * new_radiance
        end
    end

    return (emitted_radiance + cum_radiance * (1.f0 / rend.num_of_rays))
end

## Point-light render algorithm ###############################################
struct PointLight_Renderer <: Renderer
    world::World
    background_color::RGB

    PointLight_Renderer(world::World, background_color::RGB=RGB(0.1f0,0.1f0,0.1f0) ) = new(world, background_color)
end

function PointLight(ray::Ray, renderer::PointLight_Renderer)
    hit_record = ray_intersection(renderer.world, ray)
    if hit_record === nothing
        return renderer.background_color
    else
        hit_point=hit_record.world_point
        radiance = RGB(0.,0.,0.)
        min_luminosity= typemax(Float32)

            for i ∈ 1:length(renderer.world.shapes)

            if isa(renderer.world.shapes[i], LightPoint)
                new_ray = Ray(hit_point, renderer.world.shapes[i].point-hit_point)
                is_illuminated = ray_intersection(renderer.world, new_ray)
                min_luminosity = min(min_luminosity, renderer.world.shapes[i].material.emission_intensity)
                if is_illuminated === nothing
                    hit_color = get_color(hit_record.shape.material.brdf.pigment, hit_record.surface_point)
                    emitted_radiance = get_color(hit_record.shape.material.emitted_radiance, hit_record.surface_point)
                    illumination = get_color(renderer.world.shapes[i].material.emitted_radiance, hit_record.surface_point) * renderer.world.shapes[i].material.emission_intensity
                    radiance += (hit_color + emitted_radiance + illumination) * 0.33

                elseif norm(is_illuminated.world_point-new_ray.origin) > norm(renderer.world.shapes[i].point-hit_point)
                    hit_color = get_color(hit_record.shape.material.brdf.pigment, hit_record.surface_point)
                    emitted_radiance = get_color(hit_record.shape.material.emitted_radiance, hit_record.surface_point)
                    illumination = get_color(renderer.world.shapes[i].material.emitted_radiance, hit_record.surface_point) * renderer.world.shapes[i].material.emission_intensity
                    radiance += (hit_color + emitted_radiance + illumination) * 0.33
                end
            end

        end

        if radiance == RGB(0.,0.,0.)
            hit_color = get_color(hit_record.shape.material.brdf.pigment, hit_record.surface_point)
            emitted_radiance = get_color(hit_record.shape.material.emitted_radiance, hit_record.surface_point)
            radiance += (hit_color + emitted_radiance) *  min_luminosity * 0.25
        end
        return radiance
    end
end


text:
export NeuralNetwork, NNPath

# const NNEnsemble{F, B} = FWEnsemble{F, B}

function _randominit(prob::Problem{L, B, CS}, C::F, nnodes::Int)  where {L, B<:AbSignedUniDirBasis, CS, F<:AbstractFloat}
    ens = FWEnsemble(prob)
    β = ones(nnodes)*√C/nnodes #projectl1ball(randn(nnodes), C/2.)
    for j in 1:nnodes
        u = project(prob, randn(sizep(prob)))
        sgn = rand([-1, 1])
        push!(ens, B(u, sgn), 0.*β[j])
    end
    ens
end

NNPath(prob::Problem{L, B, CS}, Cs::Vector{F}, nnodes::Int;
       kwargs...) where {F, L, B<:AbBasis{F}, CS} =
        [NeuralNetwork(prob, c, nnodes; kwargs...) for c in Cs]

function NeuralNetwork(prob::Problem{L, B, CS},
              C::F,
              nnodesorens::Union{Int,FWEnsemble};
              stepsize=1e-2, βstep=:safe,
              maxiter=5000) where {L, B, CS, F<:AbstractFloat}
    if isa(nnodesorens, Int)
        nnodes = nnodesorens
        ens = _randominit(prob, C, nnodes)
    else
        ens = deepcopy(nnodesorens)
        nnodes = length(ens)
    end
    ∇L_u = Array{F}(sizep(prob), nnodes)
    ∇L_β = similar(ens.β)
    currloss = trainerr(ens)
    iter = 0
    # llog(txt="") = printfmtln("{:3d}) {:6.6f} {:.2f}({}) {}", iter, currloss, sum(abs, ens.β), countnz(ens.β), txt)
    llog(t="") = nothing
    llog("Initial")
    errs = fill(NaN, 2maxiter+1)
    errs[1] = currloss
    lim = NaN
    projected = false
    σ1sqrd = NaN
    if βstep === :safe
        z=getZ(ens)
        ZZ=z'*z
        σ1sqrd = eig(ZZ)[1][end]
    elseif βstep === :limited
        lim = 16C/maxiter/nnodes
    end

    for iter in 1:maxiter
        prevloss = currloss
        step = stepsize/√iter
        ∇L_z = derv(L(), prob.y, ens.fit)
        for j in 1:nnodes
            ∇L_u[:, j] = ens.β[j] * jacobian(ens.ωs[j], prob.X)' * ∇L_z
        end
        stepalong!(ens, ∇L_u, step)
        # sanitycheck(ens)
        currloss = trainerr(ens)
        llog("new {ω}")
        errs[2iter] = currloss
        if βstep === :lasso
            Z = getZ(ens)
            path = glmnet!(Z, _gety(prob), getdist(L), intercept=false, standardize=false)
            β = getβfromC(path.betas, C)
        else
            ∇L_β = [ens.zs[j]' * ∇L_z for j in 1:nnodes]
            if βstep === :safe
                Δβ = -∇L_β/σ1sqrd
            elseif βstep === :verysafe
                z=getZ(ens);  ZZ=z'*z;  σ1sqrd = eig(ZZ)[1][end]
                Δβ = -∇L_β/σ1sqrd
            elseif βstep === :newton
                z=getZ(ens);  ZZ=z'*z
                Δβ = -ZZ\∇L_β
            elseif βstep === :limited
                ∇L_β /= maximum(abs, ∇L_β) / lim
                Δβ = -∇L_β
                println(sum(Δβ.*ens.β .> 0))
            elseif βstep===:none
                Δβ = -step*∇L_β
            else
                error("βstep not understood: ", βstep)
            end
            β = ens.β + Δβ
            projected = sum(abs, β) > C
            false && println(round.(Int, 1*β), sum(abs, β))
            β = projectl1ball(β, C)
            false && println(round.(Int, 1*β), sum(abs, β))
        end
        updateβ!(ens, β, false)
        # sanitycheck(ens)
        currloss = trainerr(ens)
        llog(format("new β |{}| |{:.1f}| {}", countnz(β), sum(abs, β), projected?"!":""))
        errs[2iter+1] = currloss
        (βstep===:lasso) && (abs(1-currloss/prevloss) < 1e-5) && (println("loss converged"); break)
    end
    ens #, errs
end