text:
==′(x::T, y::T) where {T} = x == y
!=′(x::T, y::T) where {T} = x != y

assertion_enabled() = false
# assertion_enabled() = true

function enable_assertion()
    prev = assertion_enabled()
    @eval assertion_enabled() = true
    return prev
end

@noinline unreachable() = error("unreachable reached")

@noinline unexpected(x) = error("unexpected value: $x")

@noinline static_error(::Val{msg}) where {msg} = error(string(msg))

macro static_error(msg::AbstractString)
    sym = Symbol(msg)
    :(static_error(Val{$(QuoteNode(sym))}()))
end

function noinline(f)
    @noinline g() = f()
    return g()
end

"""
    is_pointerfree_type(T::Type) :: Bool

Return `true` if any instances of `T` do not contain boxed Julia objects.
"""
is_pointerfree_type(::Type{T}) where {T} = isconcretetype(T) && Base.datatype_pointerfree(T)

function ceillog2(n::Integer)
    n > 0 || throw(DomainError(n))
    i = trailing_zeros(n)
    j = 8 * sizeof(n) - leading_zeros(n) - 1
    if i == j
        return i
    else
        return j + 1
    end
end

@inline function uint_for(::Type{T}) where {T}
    if sizeof(T) <= 1
        return UInt8
    elseif sizeof(T) <= 2
        return UInt16
    elseif sizeof(T) <= 4
        return UInt32
    elseif sizeof(T) <= 8
        return UInt64
    elseif sizeof(T) <= 16
        return UInt128
    else
        @static_error("size of type is larger than 16 bytes")
    end
end

# TODO: use memcpy?
function uint_from(x::T) where {T}
    UIType = uint_for(T)
    ref = Ref(zero(UIType))
    GC.@preserve ref begin
        unsafe_store!(Ptr{T}(pointer_from_objref(ref)), x)
    end
    return ref[]
end

mutable struct PaddedRef{T,Pad}
    x::T
    pad::NTuple{Pad,UInt8}
    PaddedRef{T,Pad}() where {T,Pad} = new{T,Pad}()
end

@inline paddedref(::Type{T}, ::Type{Desired}) where {T,Desired} =
    PaddedRef{T,sizeof(Desired) - sizeof(T)}()

function from_bytes(::Type{T}, uint::UIntType) where {T,UIntType}
    Base.allocatedinline(T) || @static_error("expected isbits or union of them")
    let ref = paddedref(T, UIntType)
        GC.@preserve ref begin
            unsafe_store!(Ptr{UIntType}(pointer_from_objref(ref)), uint)
        end
        return ref.x
    end
end

function UnsafeAtomics.load(p::Ptr{T}, ord::Val) where {T}
    q = Ptr{uint_for(T)}(p)
    uint = UnsafeAtomics.load(q, ord)
    return from_bytes(T, uint)
end

function UnsafeAtomics.store!(p::Ptr{T}, v::T, ord::Val) where {T}
    uint = uint_from(v)
    q = Ptr{typeof(uint)}(p)
    UnsafeAtomics.store!(q, uint, ord)
end

function UnsafeAtomics.cas!(p::Ptr{T}, cmp::T, new::T, so::Val, fo::Val) where {T}
    ci = uint_from(cmp)
    ni = uint_from(new)
    q = Ptr{typeof(ci)}(p)
    oi = UnsafeAtomics.cas!(q, ci, ni, so, fo)
    return from_bytes(T, oi)
end

# Read /sys/devices/system/cpu/cpu0/cache/index0/coherency_line_size?
const CACHELINE_SIZE = 64

# Workaround
# https://github.com/JuliaconCurrent/ConcurrentCollections.jl/issues/35
# https://github.com/JuliaLang/julia/issues/42326
#=
primitive type PadAfter64 448 end
=#
const PadAfter64 = NTuple{CACHELINE_SIZE - sizeof(UInt64),UInt8}
_PadAfter64() = Ref{PadAfter64}()[]

mutable struct CheckPadAfter64
    a::UInt64
    pad::PadAfter64
    b::UInt64
end
if fieldoffset(CheckPadAfter64, 3) != CACHELINE_SIZE
    @error "Unexpected padding with `PadAfter64`" fieldoffset(CheckPadAfter64, 3)
end

#=
const PadAfter32 = PadAfter64
=#
const PadAfter32 = NTuple{CACHELINE_SIZE - sizeof(UInt32),UInt8}
_PadAfter32() = Ref{PadAfter32}()[]

mutable struct CheckPadAfter32
    a::UInt32
    pad::PadAfter32
    b::UInt32
end
if fieldoffset(CheckPadAfter32, 3) != CACHELINE_SIZE
    @error "Unexpected padding with `PadAfter32`" fieldoffset(CheckPadAfter32, 3)
end

function cacheline_padded_vector(::Type{T}, n::Integer) where {T}
    cacheline = cld(sizeof(T), CACHELINE_SIZE)
    xs = Vector{T}(undef, cacheline * (n + 1))
    ys = view(xs, cacheline+1:cacheline:length(xs))
    @assert length(ys) == n
    return ys
end

function threaded_foreach(f, xs)
    y = iterate(xs)
    y === nothing && return
    x1, st = y
    y = iterate(xs, st)
    if y === nothing
        f(x1)
        return
    end
    tasks = Task[]
    while true
        local x
        x, st = y
        t = Threads.@spawn f(x)
        push!(tasks, t)
        y = iterate(xs, st)
        y === nothing && break
    end
    f(x1)
    foreach(wait, tasks)
    return
end

function is_singleton_collection(xs)
    y = iterate(xs)
    y === nothing && return false
    _, st = y
    y = iterate(xs, st)
    y === nothing && return true
    return false
end

function threaded_typed_mapreduce(f, ::Type{T}, op, xs; kw...) where {T}
    is_singleton_collection(xs) && return f(first(xs))
    refs = [Ref{T}() for _ in xs]
    threaded_foreach(zip(refs, xs)) do (y, x)
        y[] = f(x)
    end
    return mapreduce(getindex, op, refs; kw...)
end

function zip_strict(a, args...)
    lens = map(length, args)
    all(==(length(a)), lens) || noinline() do
        error("collections have non-identical length: ", length(a), ", ", join(lens, ", "))
    end
    return zip(a, args...)
end

function define_docstrings()
    docstrings = [:ConcurrentCollections => joinpath(dirname(@__DIR__), "README.md")]
    docsdir = joinpath(@__DIR__, "docs")
    for filename in readdir(docsdir)
        stem, ext = splitext(filename)
        ext == ".md" || continue
        name = Symbol(stem)
        name in names(ConcurrentCollections, all = true) || continue
        push!(docstrings, name => joinpath(docsdir, filename))
    end
    for (name, path) in docstrings
        include_dependency(path)
        doc = read(path, String)
        doc = replace(doc, r"^```julia"m => "```jldoctest $name")
        doc = replace(doc, "<kbd>TAB</kbd>" => "_TAB_")
        @eval ConcurrentCollections $Base.@doc $doc $name
    end
end


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons, Josef Kircher
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

using Plots

# our simulation setup
t_start = 0.0
t_stop = 8.0

# load the FMU container
myFMU = fmiLoad("SpringPendulum1D", ENV["EXPORTINGTOOL"], ENV["EXPORTINGVERSION"])

# print some useful FMU-information into the REPL
fmiInfo(myFMU)

# make an instance from the FMU
fmiInstantiate!(myFMU)

recordValues = ["mass.s", "mass.v"]
solutionME = fmiSimulateME(myFMU, t_start, t_stop; recordValues=recordValues)
solutionCS = fmiSimulateCS(myFMU, t_start, t_stop; recordValues=recordValues)

# plot the results
fig = fmiPlot(solutionME)

fig = Plots.plot()
fmiPlot!(fig, solutionME)

fig = fmiPlot(solutionCS)

fig = Plots.plot()
fmiPlot!(fig, solutionCS)

# unload the FMU, remove unpacked data on disc ("clean up")
fmiUnload(myFMU)

ext: jl


text:
using MicrochipRTDM
using FTD2XX
#=
38400
57600
115200
128000
224000
256000
460800
921600
=#
projectdirectory = "c:/Users/Chris/MPLABXProjects/lbcmcp/PWM_RTDM.x"
mapdict = microchip_parsemap(projectdirectory,"production")
uartconfig1 = UARTConfiguration(115200,8,1,"n",readtimeout = 1000,writetimeout = 1000)
io = open(FT_DeviceIndex(0),uartconfig1)
setlatencytimer(io,2); #Should speed things up a little
rtdminterface = RTDMInterface(mapdict,io)
const r = 1
const address = mapdict[:offsetP3]
const buffer = Array(UInt16,1)
const buffer2 = Array(UInt8,128)
const buffer3 = Array(UInt8,9)
i = UInt16(35)  # 35 is bad value

rtdm_write(rtdminterface, i, :DTRB, retry = r)

close(rtdminterface)


text:
using LinearAlgebra, StatsFuns, Nonconvex
using LinearAlgebra: norm
Nonconvex.@load Ipopt
# Nonconvex.@load NLopt
# Nonconvex.@load Multistart

# Point load in the middle, kN
P = -100.0
# half-span, meter
h_span = 2.0

# O------C-------S
# \    / | \    /
#  \  /  |  \  /
#   \/   |   \/
#   X--------X'
# O = (0,0), X = (x,y) the variables
# norm(O,S) = 2*h_span
# P is a vertical point load applied at C
# X' node is mirrored

x0 = [1.0, -1.0]

function obj(X)
    x = X[1:2]
    # x = (x,y) of the variable node
    # compute the load path as a proxy for necessary structural weight
    L = [h_span, norm(x), norm(x-[h_span,0.0]), abs(h_span-x[1])*2]
    cosθ, sinθ = x/L[2]
    cosϕ, sinϕ = (x-[h_span,0.0])/L[3]
    F2 = 0.5*P/sinθ
    F3 = -0.5*P/sinϕ
    F = [-F2*cosθ, F2, F3, F2*cosθ+F3*cosϕ]
    FL = abs.(F).*L
    return 2*sum(FL[1:3]) + FL[4]
end

xstar = [1.0, -1.5]
shift = 1.0
power = 2.0
function deflation_constr(X)
    return 1.0/norm(X[1:2]-xstar, power)^power + shift - X[3]
end

model = Model(obj)
# X = (x1,y1,y), where y is the deflation auxilary variable
addvar!(model, [0.0, -3.0], [h_span, 3.0])
# add_ineq_constraint!(model, deflation_constr)

alg = IpoptAlg()
options = IpoptOptions(tol = 1e-4)
r1 = optimize(model, alg, x0, options = options)
@show(r1.minimizer, r1.minimum)

xstar = r1.minimizer
addvar!(model, [-1e3], [1e6])
add_ineq_constraint!(model, deflation_constr)
r2 = optimize(model, alg, vcat(x0, 1.0), options = options)
@show(r2.minimizer, r2.minimum)

###########################

# alg = NLoptAlg(:LD_MMA)
# options = NLoptOptions()
# r = optimize(model, alg, vcat(x0, 1.0), options = options)
# println(r.minimizer, r.minimum)

###########################

# options = DeflatedOptions(
#     ndeflations = 10, sub_options = IpoptOptions(tol = 1e-4), radius = 2.0,
# )
# alg = DeflatedAlg(IpoptAlg())
# @time r = Nonconvex.optimize(model, alg, x0, options = options)


text:
# compute ranges of consistency bars
# multiple bins
function consistency_ranges(
    probabilities::DataStructures.OrderedDict{Int,<:AbstractVector{<:Real}}, bars
)
    results = map(values(probabilities)) do probs
        return consistency_ranges(probs, bars)
    end
    return results
end
# single bin
function consistency_ranges(probabilities::AbstractVector{<:Real}, bars::ConsistencyBars)
    probs = similar(probabilities)
    targets = Vector{Bool}(undef, length(probabilities))
    samples = StructArrays.StructArray((probs, targets))
    spl = Random.Sampler(Random.GLOBAL_RNG, ConsistencyResampling.Consistent(probabilities))

    frequencies = map(1:(bars.samples)) do _
        Random.rand!(Random.GLOBAL_RNG, samples, spl)
        return Statistics.mean(targets)
    end

    return StatsBase.quantile(frequencies, bars.quantiles)
end

# no consistency bars
function consistency_ranges(
    probabilities::DataStructures.OrderedDict{Int,<:AbstractVector{<:Real}}, ::Nothing
)
    return fill((NaN, NaN), length(probabilities))
end

# compute mean predictions, frequencies and ranges of consistency bars
function means_and_bars(
    probabilities::AbstractVector{<:Real},
    outcomes::AbstractVector{Bool};
    binning,
    consistencybars,
)
    length(probabilities) == length(outcomes) ||
        error("number of probabilities and outcomes must be equal")

    # obtain bins
    bins = ReliabilityDiagrams.bins(probabilities, outcomes, binning)
    maxbins = max_nonempty_bins(bins, length(probabilities))

    # compute dictionaries of binned probabilities and mean frequencies of outcomes
    T = eltype(probabilities)
    dict_probabilities = DataStructures.OrderedDict{Int,Vector{T}}()
    dict_meanfrequencies = DataStructures.OrderedDict{Int,Float64}()
    sizehint!(dict_probabilities, maxbins)
    sizehint!(dict_meanfrequencies, maxbins)
    for (probability, outcome) in zip(probabilities, outcomes)
        # compute index of bin
        index = binindex(bins, probability)

        # append to the array of probabilities in the current bin
        bin = get!(dict_probabilities, index, Vector{T}(undef, 0))
        push!(bin, probability)

        # update mean frequency for the current bin
        n = length(bin)
        meanfrequency = get(dict_meanfrequencies, index, 0.0)
        dict_meanfrequencies[index] = meanfrequency + (outcome - meanfrequency) / n
    end

    # compute mean probabilities
    meanprobabilities = map(Statistics.mean, values(dict_probabilities))

    # extract mean frequencies (has the same order!)
    meanfrequencies = collect(values(dict_meanfrequencies))

    # compute consistency ranges (in the same order)
    low_high = consistency_ranges(dict_probabilities, consistencybars)

    # sort results according to mean probabilities
    perm = sortperm(meanprobabilities)
    sorted_meanprobabilities = meanprobabilities[perm]
    sorted_meanfrequencies = meanfrequencies[perm]
    sorted_low_high = low_high[perm]

    return sorted_meanprobabilities, sorted_meanfrequencies, sorted_low_high
end


text:
### BAM APIs

### bam header
function bam_hdr_init()
    @fncall(:bam_hdr_init, Ptr{Bam_hdr_t},())
end

function bam_hdr_read(fp::Ptr{BGZF})
    @fncall(:bam_hdr_read,Ptr{Bam_hdr_t},(Ptr{BGZF},),fp)
end

function bam_hdr_write(fp::Ptr{BGZF}, h::Ptr{Bam_hdr_t})
    @fncall(:bam_hdr_write, Int32, (Ptr{BGZF}, Ptr{Bam_hdr_t}), fp, h)
end

function bam_hdr_destroy(h::Ptr{Bam_hdr_t})
    @fncall(:bam_hdr_destroy, Void, (Ptr{Bam_hdr_t},), h)
end

function bam_name2id(h::Ptr{Bam_hdr_t}, ref::Ptr{Char})
    @fncall(:bam_name2id, Int32, (Ptr{Bam_hdr_t},Ptr{Char}), h, ref)
end

function bam_hdr_dup(h0::Ptr{Bam_hdr_t})
    @fncall(:bam_hdr_dup, Ptr{Bam_hdr_t}, (Ptr{Bam_hdr_t},), h0)
end

### bgzf

function bgzf_is_bgzf(fn::Ptr{Cchar})
    flag = @fncall(:bgzf_is_bgzf, Int32, (Ptr{Cchar},), fn)
    if flag == 0
        error("fn is not BGZF, OR I/O error")
    end
    true
end

function bgzf_open(path::Ptr{Cchar}, mode::Ptr{Cchar})
    @fncall(:bgzf_open, Ptr{BGZF}, (Ptr{Cchar},Ptr{Cchar}), path, mode)
end

function bgzf_close(fp::Ptr{BGZF})
    @fncall(:bgzf_close, Int32, (Ptr{BGZF},), fp)
end

function bgzf_eof(fp::Ptr{BGZF})
    error("return flag error")
    flag = @fncall(:bgzf_check_EOF, Int32, (Ptr{BGZF},), fp)
    info("bgzf_eof $flag")
    if flag == 2
        error("bgzf isn't seekable")
    end
    if flag == -1
        error("on error")
    end
    if flag == 0
        return false
    end
    if flag == 1
        return true
    end
end

function bgzf_getline!(fp::Ptr{BGZF}, delim::Int32, str::Ptr{KStr})
    @fncall(:bgzf_getline, Int32, (Ptr{BGZF}, Int32, Ptr{KStr}), fp, delim, str)
end

### bam records
function bam_init1()
    @fncall(:bam_init1, Ptr{Bam1_t},())
end

function bam_destroy1(b::Ptr{Bam1_t})
    @fncall(:bam_destroy1, Void, (Ptr{Bam1_t},), b)
end

function bam_read1!(fp::Ptr{BGZF},b::Ptr{Bam1_t})
    @fncall(:bam_read1, Int32, (Ptr{BGZF},Ptr{Bam1_t}), fp, b)
end

function bam_write1(fp::Ptr{BGZF},b::Ptr{Bam1_t})
    @fncall(:bam_write1, Int32, (Ptr{BGZF},Ptr{Bam1_t}), fp, b)
end

function bam_copy1(bdst::Ptr{Bam1_t}, bsrc::Ptr{Bam1_t})
    @fncall(:bam_copy1, Ptr{Bam1_t}, (Ptr{Bam1_t},Ptr{Bam1_t}), bdst, bsrc)
end

function bam_dup1(bsrc::Ptr{Bam1_t})
    @fncall(:bam_dup1, Ptr{Bam1_t},(Ptr{Bam1_t},),bsrc)
end

function bam_cigar2qlen(n_cigar::Int32, cigar::Ptr{UInt32})
    @fncall(:bam_cigar2qlen, Int32, (Int32, Ptr{UInt32}), n_cigar,cigar)
end

function bam_cigar2rlen(n_cigar::Int32, cigar::Ptr{UInt32})
    @fncall(:bam_cigar2rlen, Int32, (Int32, Ptr{UInt32}), n_cigar,cigar)
end

function bam_endpos(b::Ptr{Bam1_t})
    @fncall(:bam_endpos, Int32, (Ptr{Bam1_t},), b)
end

function bam_str2flag(str::Ptr{Cchar})
    @fncall(:bam_str2flag, Int32, (Ptr{Cchar},), str)
end

function bam_flag2str(flag::Int32)
    @fncall(:bam_flag2str, Ptr{Cchar}, (Int32,), flag)
end

function bam_itr_destroy(iter::Ptr{Hts_itr_t})
    @fncall(:hts_itr_destroy,Void,(Ptr{Hts_itr_t},), iter)
end

function bam_itr_queryi(idx::Ptr{Hts_idx_t},tid::Cint,bg::Cint,ed::Cint)
    @fncall(:sam_itr_queryi, Ptr{Hts_itr_t}, (Ptr{Hts_idx_t}, Cint, Cint, Cint), idx, tid, bg, ed)
end

function bam_itr_querys(idx::Ptr{Hts_idx_t}, hdr::Ptr{Bam_hdr_t}, region::Ptr{Cchar})
    @fncall(:sam_itr_querys, Ptr{HTs_itr_t}, (Ptr{Hts_idx_t},Ptr{Bam_hdr_t},Ptr{Cchar}), idx, hdr, region)
end

function bam_itr_next(htsfp::Ptr{Htsfile}, itr::Ptr{Hts_itr_t}, r::Ptr{Void})
    bgzf = unsafe_load(unsafe_load(htsfp).fp).bgzf
    @fncall(:hts_itr_next, Cint, (Ptr{BGZF}, Ptr{Hts_itr_t}, Ptr{Void}, Ptr{Void}), bgzf, itr, r, C_NULL)
end

function sam_index_load(fp::Ptr{HtsFile}, fn::Ptr{Cchar})
    @fncall(:sam_index_load, Ptr{Hts_idx_t}, (Ptr{HtsFile}, Ptr{Cchar}), fp, fn)
end

function sam_index_load2(fp::Ptr{HtsFile}, fn::Ptr{Cchar}, fnidx::Ptr{Cchar})
    @fncall(:sam_index_load2, Ptr{Hts_idx_t}, (Ptr{HtsFile}, Ptr{Cchar}, Ptr{Cchar}), fp, fn, fnidx)
end

function sam_index_build(fn::Ptr{Cchar}, min_shift::Cint)
    @fncall(:sam_index_build, Cint, (Ptr{Cchar}, Cint), fn, min_shift)
end

function sam_index_build2(fn::Ptr{Cchar}, fnidx::Ptr{Cchar}, min_shift::Cint)
    @fncall(:sam_index_build2, Cint, (Ptr{Cchar}, Ptr{Cchar}, Cint), fn, fnidx, min_shift)
end

function sam_itr_destroy(iter)
    bam_itr_destroy(iter)
end

function sam_itr_queryi(idx::Ptr{Hts_idx_t}, tid::Cint, bg::Cint, ed::Cint)
    @fncall(:sam_itr_queryi, Ptr{Hts_itr_t}, (Ptr{Hts_idx_t}, Cint, Cint, Cint), idx, tid, bg, ed)
end

function sam_itr_querys(idx::Ptr{Hts_idx_t},hdr::Ptr{Bam_hdr_t},region::Ptr{Cchar})
    @fncall(:sam_itr_querys, Ptr{Hts_itr_t}, (Ptr{Hts_idx_t}, Ptr{Bam_hdr_t}, Ptr{Cchar}), idx, hdr, region))
end

function sam_itr_next(htsfp::Ptr{HtsFile},itr::Ptr{Hts_itr_t},r::Cint)
    bam_itr_next(htsfp, itr, r)
end

function sam_open(fn::Ptr{Cchar}, mode::Ptr{Cchar})
    @fncall(:hts_open, Ptr{HtsFile}, Ptr{Cchar}, Ptr{Cchar}), fn, mode)
end

function sam_open_format(fn::Ptr{Cchar}, mode::Ptr{Cchar}, format::Ptr{Cchar})
    @fncall(:hts_open_format, Ptr{Cchar}, (Ptr{Char}),Ptr{Char},Ptr{Char}), fn, mode, format)
end

function sam_close(fp::Ptr{HtsFile})
    @fncall(:hts_close, Cint, (Ptr{HtsFile},), fp)
end

function sam_open_mode(mode::Ptr{Cchar}, fn::Ptr{Cchar}, format::Ptr{Char})
    @fncall(:sam_open_mode, Cint, (Ptr{Cchar},Ptr{Cchar},Ptr{Cchar}),mode,fn,format)
end

function sam_hdr_parse(l_text::Cint, text::Ptr{Cchar})
    @fncall(:sam_hdr_parse, Ptr{Bam_hdr_t},(Cint, Ptr{Char}), l_text, text)
end

function sam_hdr_read(fp::Ptr{SamFile})
    @fncall(:sam_hdr_read, Ptr{Bam_hdr_t}, (Ptr{SamFile},), fp)
end

function sam_hdr_write(fp::Ptr{SamFile}, h::Ptr{Bam_hdr_t})
    @fncall(:sam_hdr_write, Cint, (Ptr{SamFile}, Ptr{Bam_hdr_t}), fp, h)
end

function sam_parse1(s::Ptr{KStr}, h::Ptr{Bam_hdr_t}, b::Ptr{Bam1_t})
    @fncall(:sam_parse1, Cint, (Ptr{KStr}, Ptr{Bam_hdr_t}, Ptr{Bam1_t}), s, h, b)
end

function sam_format1!(h::Ptr{Bam_hdr_t}, b::Ptr{Bam1_t}, kstr::Ptr{KStr})
    @fncall(:sam_format1,Cint,(Ptr{Bam_hdr_t},Ptr{Bam1_t},Ptr{KStr}), h, b, kstr)
end

function sam_read1(fp::Ptr{SamFile}, h::Ptr{Bam_hdr_t}, b::Ptr{Bam1_t})
    @fncall(:sam_read1, Cint, (Ptr{SamFile}, Ptr{Bam_hdr_t}, Ptr{Bam1_t}), fp, h, b)
end

function sam_write1(fp::Ptr{SamFile}, h::Ptr{Bam_hdr_t}, b::Ptr{Bam1_t})
    @fncall(:sam_write1, Cint, (Ptr{SamFile}, Ptr{Bam_hdr_t}, Ptr{Bam1_t}), fp, h, b)
end

### Manipulating auxilary fields

function bam_aux_get(b::Ptr{Bam1_t})
    tag = Array{Cchar,1}(2)
    ptag = pointer(tag)
    @fncall(:bam_aux_get, Ptr{UInt8}, (Ptr{Bam1_t}, Ptr{Cchar}), b, ptag)
end

function bam_aux2i(s::Ptr{UInt8})
    @fncall(:bam_aux2i, Int32, (Ptr{UInt8},), s)
end

function bam_aux2f(s::Ptr{UInt8})
    @fncall(:bam_aux2f, Cdouble, (Ptr{UInt8},), s)
end

function bam_aux2A(s::Ptr{UInt8})
    @fncall(:bam_aux2A, Cchar, (Ptr{UInt8},), s)
end

function bam_aux2Z(s::Ptr{UInt8})
    @fncall(:bam_aux2Z, Ptr{Cchar}, (Ptr{UInt8},), s)
end

function bam_aux_append(b::Ptr{Bam1_t}, tag::Array{Cchar,1}, typ::Cchar, len::Cint, data::Ptr{UInt8})
    @assert length(tag) == 2
    ptag = pointer(tag)
    @fncall(:bam_aux_append, Void, (Ptr{Bam1_t}, Ptr{Cchar}, Cchar, Cint, Ptr{UInt8}), b, ptag, typ, len, data)
end

function bam_aux_del(b::Ptr{Bam1_t}, s::Ptr{UInt8})
    @fncall(:bam_aux_del, Cint, (Ptr{Bam1_t}, Ptr{UInt8}), b, s)
end

### Pileup and Mpileup
#typedef int (*bam_plp_auto_f)(void *data, bam1_t *b)

function bam_plp_init(func::Function, data::Ptr{Void} )
    @fncall(:bam_plp_init, Bam_plp_t, )
end


text:
struct AnalysisParams
    # disables caching of native remarks (that may speed up profiling time)
    filter_native_remarks::Bool

    function AnalysisParams(; filter_native_remarks::Bool = true,
                              # dummy kwargs so that kwargs for other functions can be passed on
                              __kwargs...)
        return new(filter_native_remarks,
                   )
    end
end

const LocalCache = Dict{Vector{Any},Vector{InferenceErrorReportCache}}

struct JETInterpreter <: AbstractInterpreter
    #= native =#

    native::NativeInterpreter
    optimize::Bool
    compress::Bool
    discard_trees::Bool

    #= JET.jl specific =#

    # local report cache for constant analysis
    cache::LocalCache

    # for sequential assignment of virtual global variables
    id::Symbol

    # reports found so far
    reports::Vector{InferenceErrorReport}

    # stashes `UncaughtExceptionReport`s that are not caught so far
    uncaught_exceptions::Vector{UncaughtExceptionReport}

    # stashes `NativeRemark`s
    native_remarks::Vector{NativeRemark}

    # toplevel profiling (skip inference on actually interpreted statements)
    concretized::BitVector

    # configurations for analysis performed by this interpreter
    analysis_params::AnalysisParams

    # debugging
    depth::Ref{Int}

    function JETInterpreter(world               = get_world_counter();
                            inf_params          = gen_inf_params(),
                            opt_params          = gen_opt_params(),
                            optimize            = true,
                            compress            = false,
                            discard_trees       = false,
                            id                  = gensym(:JETInterpreterID),
                            reports             = InferenceErrorReport[],
                            uncaught_exceptions = UncaughtExceptionReport[],
                            native_remarks      = NativeRemark[],
                            concretized         = BitVector(),
                            analysis_params     = AnalysisParams(),
                            # dummy kwargs so that kwargs for other functions can be passed on
                            __kwargs...)
        @assert !opt_params.inlining "inlining should be disabled for JETInterpreter analysis"

        native = NativeInterpreter(world; inf_params, opt_params)
        return new(native,
                   optimize,
                   compress,
                   discard_trees,
                   LocalCache(),
                   id,
                   reports,
                   uncaught_exceptions,
                   native_remarks,
                   concretized,
                   analysis_params,
                   Ref(0),
                   )
    end
end

# AbstractInterpreter API
# -----------------------

CC.InferenceParams(interp::JETInterpreter) = InferenceParams(interp.native)
CC.OptimizationParams(interp::JETInterpreter) = OptimizationParams(interp.native)
CC.get_world_counter(interp::JETInterpreter) = get_world_counter(interp.native)
CC.get_inference_cache(interp::JETInterpreter) = get_inference_cache(interp.native)

# JET only works for runtime inference
CC.lock_mi_inference(::JETInterpreter, ::MethodInstance) = nothing
CC.unlock_mi_inference(::JETInterpreter, ::MethodInstance) = nothing

function CC.add_remark!(interp::JETInterpreter, sv::InferenceState, s::String)
    AnalysisParams(interp).filter_native_remarks && return
    push!(interp.native_remarks, NativeRemark(interp, sv, s))
    return
end

CC.may_optimize(interp::JETInterpreter) = interp.optimize
CC.may_compress(interp::JETInterpreter) = interp.compress
CC.may_discard_trees(interp::JETInterpreter) = interp.discard_trees

# JETInterpreter specific
# -----------------------

AnalysisParams(interp::JETInterpreter) = interp.analysis_params

function gen_inf_params(; # more constant prop, more correct reports ?
                          aggressive_constant_propagation::Bool = true,
                          # turn this off to get profiles on `throw` blocks, this might be good to default
                          # to `true` since `throw` calls themselves will be reported anyway
                          unoptimize_throw_blocks::Bool = true,
                          # dummy kwargs so that kwargs for other functions can be passed on
                          __kwargs...)
    return @static VERSION ≥ v"1.6.0-DEV.837" ?
           InferenceParams(; aggressive_constant_propagation,
                             unoptimize_throw_blocks,
                             ) :
           InferenceParams(; aggressive_constant_propagation,
                             )
end

function gen_opt_params(; # inlining should be disabled for `JETInterpreter`, otherwise virtual stack frame
                          # traversing will fail for frames after optimizer runs on
                          inlining = false,
                          # dummy kwargs so that kwargs for other functions can be passed on
                          __kwargs...)
    return OptimizationParams(; inlining,
                                )
end

get_id(interp::JETInterpreter) = interp.id

# TODO do report filtering or something configured by `AnalysisParams(interp)`
function report!(interp::JETInterpreter, report::InferenceErrorReport)
    push!(interp.reports, report)
    return
end


text:
include("helpers/doc.jl")
include("helpers/page.jl")

typeexample(code, output=eval(parse(code))) =
  vbox(
     codemirror(code, readonly=true, linenumbers=false) |> fontcolor("#777"),
     output |> pad([left], 2em))

titles = vbox(
   intersperse(vskip(2em),
       map(n -> typeexample("title($n, \"Title $n\")"), 4:-1:1)))

headings = vbox(
   intersperse(vskip(2em),
       map(n -> typeexample("heading($n, \"Heading $n\")"), 1:4)))

function main(window)
    push!(window.assets, "codemirror")

    md"""
$(title(3, "Typography"))
$(vskip(1em))
Escher provides primitives that directly map to CSS font styling properties, as well as higher-level functions which form a standard typographic scale you can use to give your documents a consistent, pleasant look. By default, Escher uses the [*Source Sans Pro*](http://www.google.com/fonts/specimen/Source+Sans+Pro) (sans-serif) and [*Source Code Pro*](http://www.google.com/fonts/specimen/Source+Code+Pro) (monospaced) font families for great-looking and legible type.

# High-level Functions
$(vskip(1em))

$(h2("Titles"))

$(titles |> pad([left], 4em) |> pad([top, bottom], 1em))

$(h2("Headings"))

$(headings |> pad([left], 4em) |> pad([top, bottom], 1em))

$(vskip(1em))

    """ |> docpage
end


text:
export AdjointOperator

struct AdjointOperator{T <: SymSemiseparableMatrix} <: SymSemiseparableMatrix
	A::T
	function AdjointOperator(A::T) where {T<:SymSemiseparableMatrix}
		new{T}(A)
	end
end

# Constructor
AdjointOperator(L::AdjointOperator) = AdjointOperator(L)


text:
# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .jl
#       format_name: percent
#       format_version: '1.3'
#       jupytext_version: 1.4.2
#   kernelspec:
#     display_name: Julia 1.4.1
#     language: julia
#     name: julia-1.4
# ---

# %%
# Ensure your working directory is: https://github.com/henriquebecker91/phd/tree/master/latex/revised_PPG2KP
include("notebook_setup.jl")

# %%
# Read the data, show nothing.
lp_method_csv_path = "./data/lp_method.csv"
lp_method_df = DataFrame(CSV.File(lp_method_csv_path))
showtable(lp_method_df)
#nothing

# %%
# Clean the data a little. Try to keep this safe to re-apply to a cleaned dataframe, if possible.
lp_method_df = let
    # Keep only the instance name (not the path).
    @with(lp_method_df, :instance_name .= basename.(:instance_name))
    @with(lp_method_df, :datafile .= basename.(:datafile))
    # Replace parameter numeric values with more descriptive strings.
    lp_method_code2name = Dict{Int, String}(1 => "simplex", 2 => "barrier")
    lp_switch_code2name = Dict{Int, String}(-2 => "only", 1 => "switch")
    if eltype(lp_method_df.lp_method) == Int
        lp_method_df.lp_method = getindex.((lp_method_code2name,), lp_method_df.lp_method)
    end
    if eltype(lp_method_df.lp_method_switch) == Int
        lp_method_df.lp_method_switch = getindex.((lp_switch_code2name,), lp_method_df.lp_method_switch)
    end
    # Change the not-a-number values in final_root_time to zero. They happen
    lp_method_df
end
nothing

# %%
let lp_method_df = deepcopy(lp_method_df)
    finished = @linq lp_method_df |> select(:instance_name, :solution_profit, :finished)
    finished = @linq filter!(:finished => identity, finished) |> select!(Not(:finished))
    finished = @linq finished |> groupby(:instance_name) |> based_on(;
        solution_profit = maximum(:solution_profit), qt_distinct = length(unique(:solution_profit))
    )
    @assert all(isone, finished[!, "qt_distinct"])
    finished = select!(finished, Not(:qt_distinct))
    opt59 = DataFrame(CSV.File("./data/opt59_by_thomo.csv"))
    joined = innerjoin(finished, opt59; on = :instance_name)
    wrong = filter(:solution_profit => isequal(:best), joined)
end

# %%
# Shows the cleaned data.
println(names(lp_method_df))
showtable(lp_method_df)

# %%
# Just checking which instances ended by time limit or memory limit.
@linq lp_method_df |>
    where(.!:finished) |>
    select(:instance_name, :lp_method, :lp_method_switch, :pricing_method, :total_instance_time)

# %%
lp_params_df = let lp_method_df = deepcopy(lp_method_df)
    @show length(join.(zip(
        lp_method_df.lp_method, lp_method_df.lp_method_switch, lp_method_df.pricing_method
    ), "_"))
    @show length(lp_method_df.lp_method)
    lp_method_df[!, "params"] = join.(zip(
        lp_method_df.lp_method, lp_method_df.lp_method_switch, lp_method_df.pricing_method
    ), "_")
    pretty_params = Dict{String, String}(
        "barrier_only_none" => "NP B",
        "barrier_switch_furini" => "FP DS+B",
        "simplex_only_furini" => "FP DS",
        "barrier_only_furini" => "FP B",
        "simplex_only_none" => "NP DS"
    )
    lp_method_df[!, "params"] = getindex.((pretty_params,), lp_method_df.params)
    lp_method_df
end

# %%
lttt = let lp_params_df = lp_params_df
    sel_columns = [
        :instance_name, :params, :build_and_solve_time, :final_root_time
    ]
    # long (opposed to wide) table of total time
    lttt = select(lp_params_df, sel_columns)
    percents = @linq lttt |>
        where((:params .== "NP B") .| (:params .== "NP DS")) |>
        transform(value = (:final_root_time ./ :build_and_solve_time) .* 100.0) |>
        transform(variable = ifelse.(:params .== "NP B", "% B RNR", "% DS RNR")) |>
        select!(:instance_name, :variable, :value)

    select!(lttt, :instance_name, :params => :variable, :build_and_solve_time => :value)
    lttt = vcat(lttt, percents; cols = :setequal)

    lttt.value = (x -> @sprintf("%.2f", x)).(lttt.value)
    lttt.value = ifelse.(
        lttt.value .== "NaN", "--", string.(lttt.value)
    )

    # Table with instance, params, and time, in long format.
    lttt
end
nothing
showtable(lttt)

# %%
# Table with the build+solve time for each instance (rows) and parameter combination (columns).
# Note this counts the theoretically unnecessary root node solving of the priced models.
# wide (opposed to long) table of total time
wttt = unstack(lttt)

let
    latex_table = deepcopy(wttt)
    rename!(latex_table, :instance_name => :Name)
    desired_col_order = ["NP DS", "\\% DS RNR", "FP DS", "NP B", "\\% B RNR", "FP B", "FP DS+B"]
    rename!(latex_table, "% DS RNR" => "\\% DS RNR", "% B RNR" => "\\% B RNR")
    select!(latex_table, :Name, desired_col_order)
    highlight_best_values!(latex_table; columns = ["FP B", "FP DS", "FP DS+B", "NP B", "NP DS"])
    pretty_table(
        latex_table; backend = :latex, nosubheader = true, alignment = [:l; repeat([:r], 7)]
    )
end

# %%
# Now let us do the same but without the root node time.
lcttt = let
    sel_columns = [
        :instance_name, :lp_method, :lp_method_switch, :pricing_method, :build_and_solve_time,
        :final_root_time
    ]
    # long (opposed to wide) corrected table of total time
    lcttt = select(lp_method_df, sel_columns)
    lcttt.corrected_time = lcttt.build_and_solve_time .- lcttt.final_root_time
    lcttt.both_timings = tuple.(lcttt.build_and_solve_time, lcttt.corrected_time)
    lcttt.params = join.(zip(lcttt.lp_method, lcttt.lp_method_switch, lcttt.pricing_method), "_")
    pretty_params = Dict{String, String}(
        "barrier_only_none" => "NP B",
        "barrier_switch_furini" => "FP DS+B",
        "simplex_only_furini" => "FP DS",
        "barrier_only_furini" => "FP B",
        "simplex_only_none" => "NP DS"
    )
    lcttt.params = getindex.((pretty_params,), lcttt.params)
    lcttt = select(lcttt, :instance_name, :params, :both_timings)
    lcttt
end
showtable(lcttt)

# %%
# Table with both the build+solve time and this time minus the root node solving for each instance (rows)
# and parameter combination (columns).
# wide (opposed to long) corrected table of total time
wcttt = unstack(lcttt, :instance_name, :params, :both_timings)
clean_col(df, col) = setproperty!(df, Symbol(col), (both -> round.(both; digits = 2)).(getproperty(df, Symbol(col))))
clean_col.((wcttt,), ["FP B", "FP DS", "FP DS+B", "NP B", "NP DS"])

pretty_table(wcttt; backend = :latex)
#showtable(wcttt)

# %%
?setproperty!


text:
isolate_name(name) = split(string(name), '.')[end]

for struct_type in (:Struct, :Mutable, :CustomStruct, :OrderedStruct, :AbstractType, :DictType, :ArrayType, :StringType, :NumberType, :BoolType, :NullType)
    @eval begin
        """
            @$(isolate_name($struct_type))(expr::Expr)
            @$(isolate_name($struct_type))(expr::Symbol)

        If `expr` is a struct definition, sets the `StructType` of the defined struct to
        `$(isolate_name($struct_type))()`. If `expr` is the name of a `Type`, sets the `StructType` of that
        type to `$(isolate_name($struct_type))()`.

        # Examples
        ```julia
        @$(isolate_name($struct_type)) MyStruct
        ```
        is equivalent to
        ```julia
        StructTypes.StructType(::Type{MyStruct}) = StructType.Struct()
        ```
        and
        ```julia
        @$(isolate_name($struct_type)) struct MyStruct
            val::Int
        end
        ```
        is equivalent to
        ```julia
        struct MyStruct
            val::Int
        end
        StructTypes.StructType(::Type{MyStruct}) = StructType.Struct()
        ```
        """ macro $struct_type(expr)
            return macro_constructor(expr, $struct_type)
        end
    end
end

function macro_constructor(expr::Expr, structtype)
    if expr.head == :struct
        return esc(quote
            $expr
            StructTypes.StructType(::Type{$(expr.args[2])}) = $structtype()
        end)
    else
        return :(throw(ArgumentError("StructType macros can only be used with a struct definition or a Type")))
    end
end
function macro_constructor(expr::Symbol, structtype)
    return esc(quote
        if $(expr) isa Type
            StructTypes.StructType(::Type{$(expr)}) = $structtype()
        else
            throw(ArgumentError("StructType macros can only be used with a struct definition or a Type"))
        end
    end)
end


text:
using MixedModelsGadfly
using Test

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


text:
 using DSP, StatsBase

function ContinuousRP(ts1, ts2)

    # center the data
    ts1 = ts1 .- mean(ts1)
    ts2 = ts2 .- mean(ts2)

    # hilbert transform
    h1 = DSP.hilbert(ts1) |> imag
    h2 = DSP.hilbert(ts2) |> imag
    num = h2.*ts1 - ts2.*h1
    denom = ts2.*ts1 + h2.*h1

    # time series of CRP in radians
    radians = atan.(num,denom)

    # convert to degrees
    pRP = rad2deg.(radians)

    # circular stats in radians
    wgh = ones(size(radians))
    wr = wgh'*exp.(1im*radians)
    meanRP = angle(wr) |> rad2deg
    rvRP = abs(wr)/sum(wgh)
    sdRP = sqrt(2*(1-rvRP))

    (meanRP = meanRP, rvRP = rvRP, sdRP = sdRP, pRP = pRP)

    # Dict("meanRP" => meanRP,
    #     "rvRP" => rvRP,
    #     "sdRP" => sdRP,
    #     "pRP" => pRP)

end


text:
const YAHOO_URL = "https://query1.finance.yahoo.com/v7/finance/download"  # for querying yahoo's servers
const YAHOO_TMP = "https://ca.finance.yahoo.com/quote/^GSPC/history?p=^GSPC"  # for getting the cookies and crumbs

import HTTP
import Dates

function yahoo_get_crumb()::Tuple{SubString{String}, Dict{String, Set{HTTP.Cookies.Cookie}}}
    response = HTTP.request("GET", YAHOO_TMP)
    m = match(r"\"user\":{\"crumb\":\"(.*?)\"", String(response.body)).captures[1]
    h = response.headers
    # manually grab and parse the cookie
    function substring_after_equal(str;offset::Int64=1)
        idx = collect(findfirst("=",str))[1]
        return String(str[idx+offset:end])
    end
    function return_cookie_string(head)
        for x in head
            if occursin("cookie",lowercase(x[1]))
                return String(x[2])
            end
        end
        return nothing
    end
    cookie_str = return_cookie_string(h)
    split_h = split(cookie_str,";")
    n = String(split_h[1][1:1])
    v = String(split_h[1][3:end])
    expire = substring_after_equal(split_h[2])
    pth = substring_after_equal(split_h[3])
    dom = substring_after_equal(split_h[4],offset=2)
    c = HTTP.Cookies.Cookie(n,v,domain=dom,path=pth,expires=DateTime(expire,"e, d-u-y H:M:S G\\MT"),unparsed=[cookie_str])
    c_dict = (m,Dict(m=>Set([c])))
    return c_dict
end

"""
Download stock price data from Yahoo! Finance into a TS object.

`yahoo(symb::String; from::String="1900-01-01", thru::String=string(Dates.today()), freq::String="d", event::String="history", crumb_tuple::Tuple{SubString{String}, Dict{String, HTTP.Cookie}}=yahoo_get_crumb())::TS`

# Arguments
- `symb` ticker symbol of the stock
- `from` starting date of the historical data request (string formatted as yyyy-mm-dd)
- `thru` ending date of the historical data request (string formatted as yyyy-mm-dd)
- `freq` frequency interval of the requested dowload (valid options are \"d\" for daily, \"wk\" for weekly, and \"mo\" for monthly)
- `event` type of data download to request (valid options are \"history\" for standard historical price data, \"div\" for dividend payments, and \"split\" for stock splits)
- `crumb_tuple` workaround to provide crumbs/cookies for the new Yahoo Finance portal (which requires such data to fulfill the requests)

# Example

```
julia> yahoo("AAPL", from="2010-06-09", thru=string(Dates.today()), freq="wk")
356x6 Temporal.TS{Float64,Date}: 2010-06-09 to 2017-03-27
Index       Open    High    Low     Close   Volume      AdjClose
2010-06-09  251.47  253.86  242.2   253.51  1.813954e8  32.8446
2010-06-14  255.96  275.0   254.01  274.07  1.814594e8  35.5084
2010-06-21  277.69  279.01  265.81  266.7   1.763214e8  34.5535
2010-06-28  266.93  269.75  243.2   246.94  2.087241e8  31.9934
2010-07-06  251.0   262.9   246.16  259.62  1.525786e8  33.6362
⋮
2017-02-27  137.14  140.28  136.28  139.78  2.54267e7   139.78
2017-03-06  139.37  139.98  137.05  139.14  1.97315e7   139.14
2017-03-13  138.85  141.02  138.82  139.99  2.41057e7   139.99
2017-03-20  140.4   142.8   139.73  140.64  2.54857e7   140.64
2017-03-27  139.39  144.49  138.62  144.12  2.86449e7   144.12
```
"""
function yahoo(symb::String;
               from::String="1900-01-01",
               thru::String=string(Dates.today()),
               freq::String="d",
               event::String="history",
               crumb_tuple::Tuple{SubString{String}, Dict{String, Set{HTTP.Cookies.Cookie}}}=yahoo_get_crumb())::TS
    @assert freq in ["d","wk","mo"] "Argument `freq` must be either \"d\" (daily), \"wk\" (weekly), or \"mo\" (monthly)."
    @assert event in ["history","div","split"] "Argument `event` must be either \"history\", \"div\", or \"split\"."
    @assert from[5] == '-' && from[8] == '-' "Argument `from` has invalid date format."
    @assert thru[5] == '-' && thru[8] == '-' "Argument `thru` has invalid date format."
    period1 = Int(floor(Dates.datetime2unix(Dates.DateTime(from))))
    period2 = Int(floor(Dates.datetime2unix(Dates.DateTime(thru))))
    urlstr = "$(YAHOO_URL)/$(symb)?period1=$(period1)&period2=$(period2)&interval=1$(freq)&events=$(event)&crumb=$(crumb_tuple[1])"
    response = HTTP.request("POST",HTTP.URIs.URI(urlstr), cookies=true, cookiejar=crumb_tuple[2])
    indata = Temporal.csvresp(response)
    data = TS(indata[1], indata[2], indata[3][2:end])
    rename!(data, Symbol("Adj Close")=>:AdjClose)
    return data
end

function yahoo(syms::Vector{String};
               from::String="1900-01-01",
               thru::String=string(Dates.today()),
               freq::String="d",
               event::String="history",
               crumb_tuple::Tuple{SubString{String}, Dict{String, HTTP.Cookie}}=yahoo_get_crumb())::Dict{String,TS}
    out = Dict()
    @inbounds for s = syms
        out[s] .= yahoo(s, from=from, thru=thru, freq=freq, event=event, crumb_tuple=crumb_tuple)
    end
    return out
end


text:
# modules are about how you structure your program
# packages are about how you distribute your program `pkg.add(package)`

# modules - partioning your code into individual functions
# Julia allows you to include modules by `#include <modules>`
# modules provide a namespace for our code

# defining the start and end of module
module modules

export f

function f(x::Int)
    return x + 1
end

function g(x:Int)
    return x - 1
end

end

max_stars_repo_name: adeeconometrics/Algorithms

max_issues_repo_name: adeeconometrics/Algorithms

max_forks_repo_name: adeeconometrics/Algorithms


text:
# Script to demonstrate convergence analysis of Ecut values for DFTK.jl

using Atomistic
using AtomsBase
using DFTK
using NBodySimulator
using Unitful
using UnitfulAtomic

setup_threading(n_blas = 4)

N = 8
element = :Ar
box_size = 1.5u"nm" # this number was chosen arbitrarily
reference_temp = 94.4u"K"
thermostat_prob = 0.1 # this number was chosen arbitrarily
Δt = 1e-2u"ps"

pspkey = list_psp(:Ar, functional = "lda")[1].identifier
initial_bodies = generate_bodies_in_cell_nodes(N, element, box_size, reference_temp)
for body ∈ initial_bodies
    body.data[:pseudopotential] = pspkey
end
initial_system = FlexibleSystem(initial_bodies, CubicPeriodicBoundaryConditions(austrip(box_size)))

eq_steps = 20000
eq_thermostat = AndersenThermostat(austrip(reference_temp), thermostat_prob / austrip(Δt))
eq_simulator = NBSimulator(Δt, eq_steps, thermostat = eq_thermostat)
potential = LennardJonesParameters(1.657e-21u"J", 0.34u"nm", 0.765u"nm")

eq_result = @time simulate(initial_system, eq_simulator, potential)

display(@time plot_rdf(eq_result, potential.σ, 0.5))

dftk_potential = DFTKPotential(5u"hartree", [1, 1, 1]; damping = 0.7)

display(@time analyze_convergence(get_system(eq_result), dftk_potential, (5:2.5:25)u"hartree"))

;


text:
# Julia wrapper for header: /home/astukalov/.julia/v0.5/Sundials/deps/usr/include/cvode/cvode.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function __CVodeCreate(lmm::Cint,iter::Cint)
    ccall((:CVodeCreate,libsundials_cvode),CVODEMemPtr,(Cint,Cint),lmm,iter)
end

function CVodeCreate(lmm,iter)
    __CVodeCreate(convert(Cint,lmm),convert(Cint,iter))
end

function __CVodeSetErrHandlerFn(cvode_mem::CVODEMemPtr,ehfun::CVErrHandlerFn,eh_data::Ptr{Void})
    ccall((:CVodeSetErrHandlerFn,libsundials_cvode),Cint,(CVODEMemPtr,CVErrHandlerFn,Ptr{Void}),cvode_mem,ehfun,eh_data)
end

function CVodeSetErrHandlerFn(cvode_mem,ehfun,eh_data)
    __CVodeSetErrHandlerFn(convert(CVODEMemPtr,cvode_mem),ehfun,pointer(eh_data))
end

function __CVodeSetErrFile(cvode_mem::CVODEMemPtr,errfp::Ptr{FILE})
    ccall((:CVodeSetErrFile,libsundials_cvode),Cint,(CVODEMemPtr,Ptr{FILE}),cvode_mem,errfp)
end

function CVodeSetErrFile(cvode_mem,errfp)
    __CVodeSetErrFile(convert(CVODEMemPtr,cvode_mem),errfp)
end

function __CVodeSetUserData(cvode_mem::CVODEMemPtr,user_data::Any)
    ccall((:CVodeSetUserData,libsundials_cvode),Cint,(CVODEMemPtr,Any),cvode_mem,user_data)
end

function CVodeSetUserData(cvode_mem,user_data)
    __CVodeSetUserData(convert(CVODEMemPtr,cvode_mem),user_data)
end

function __CVodeSetMaxOrd(cvode_mem::CVODEMemPtr,maxord::Cint)
    ccall((:CVodeSetMaxOrd,libsundials_cvode),Cint,(CVODEMemPtr,Cint),cvode_mem,maxord)
end

function CVodeSetMaxOrd(cvode_mem,maxord)
    __CVodeSetMaxOrd(convert(CVODEMemPtr,cvode_mem),convert(Cint,maxord))
end

function __CVodeSetMaxNumSteps(cvode_mem::CVODEMemPtr,mxsteps::Clong)
    ccall((:CVodeSetMaxNumSteps,libsundials_cvode),Cint,(CVODEMemPtr,Clong),cvode_mem,mxsteps)
end

function CVodeSetMaxNumSteps(cvode_mem,mxsteps)
    __CVodeSetMaxNumSteps(convert(CVODEMemPtr,cvode_mem),convert(Clong,mxsteps))
end

function __CVodeSetMaxHnilWarns(cvode_mem::CVODEMemPtr,mxhnil::Cint)
    ccall((:CVodeSetMaxHnilWarns,libsundials_cvode),Cint,(CVODEMemPtr,Cint),cvode_mem,mxhnil)
end

function CVodeSetMaxHnilWarns(cvode_mem,mxhnil)
    __CVodeSetMaxHnilWarns(convert(CVODEMemPtr,cvode_mem),convert(Cint,mxhnil))
end

function __CVodeSetStabLimDet(cvode_mem::CVODEMemPtr,stldet::Cint)
    ccall((:CVodeSetStabLimDet,libsundials_cvode),Cint,(CVODEMemPtr,Cint),cvode_mem,stldet)
end

function CVodeSetStabLimDet(cvode_mem,stldet)
    __CVodeSetStabLimDet(convert(CVODEMemPtr,cvode_mem),convert(Cint,stldet))
end

function __CVodeSetInitStep(cvode_mem::CVODEMemPtr,hin::realtype)
    ccall((:CVodeSetInitStep,libsundials_cvode),Cint,(CVODEMemPtr,realtype),cvode_mem,hin)
end

function CVodeSetInitStep(cvode_mem,hin)
    __CVodeSetInitStep(convert(CVODEMemPtr,cvode_mem),hin)
end

function __CVodeSetMinStep(cvode_mem::CVODEMemPtr,hmin::realtype)
    ccall((:CVodeSetMinStep,libsundials_cvode),Cint,(CVODEMemPtr,realtype),cvode_mem,hmin)
end

function CVodeSetMinStep(cvode_mem,hmin)
    __CVodeSetMinStep(convert(CVODEMemPtr,cvode_mem),hmin)
end

function __CVodeSetMaxStep(cvode_mem::CVODEMemPtr,hmax::realtype)
    ccall((:CVodeSetMaxStep,libsundials_cvode),Cint,(CVODEMemPtr,realtype),cvode_mem,hmax)
end

function CVodeSetMaxStep(cvode_mem,hmax)
    __CVodeSetMaxStep(convert(CVODEMemPtr,cvode_mem),hmax)
end

function __CVodeSetStopTime(cvode_mem::CVODEMemPtr,tstop::realtype)
    ccall((:CVodeSetStopTime,libsundials_cvode),Cint,(CVODEMemPtr,realtype),cvode_mem,tstop)
end

function CVodeSetStopTime(cvode_mem,tstop)
    __CVodeSetStopTime(convert(CVODEMemPtr,cvode_mem),tstop)
end

function __CVodeSetMaxErrTestFails(cvode_mem::CVODEMemPtr,maxnef::Cint)
    ccall((:CVodeSetMaxErrTestFails,libsundials_cvode),Cint,(CVODEMemPtr,Cint),cvode_mem,maxnef)
end

function CVodeSetMaxErrTestFails(cvode_mem,maxnef)
    __CVodeSetMaxErrTestFails(convert(CVODEMemPtr,cvode_mem),convert(Cint,maxnef))
end

function __CVodeSetMaxNonlinIters(cvode_mem::CVODEMemPtr,maxcor::Cint)
    ccall((:CVodeSetMaxNonlinIters,libsundials_cvode),Cint,(CVODEMemPtr,Cint),cvode_mem,maxcor)
end

function CVodeSetMaxNonlinIters(cvode_mem,maxcor)
    __CVodeSetMaxNonlinIters(convert(CVODEMemPtr,cvode_mem),convert(Cint,maxcor))
end

function __CVodeSetMaxConvFails(cvode_mem::CVODEMemPtr,maxncf::Cint)
    ccall((:CVodeSetMaxConvFails,libsundials_cvode),Cint,(CVODEMemPtr,Cint),cvode_mem,maxncf)
end

function CVodeSetMaxConvFails(cvode_mem,maxncf)
    __CVodeSetMaxConvFails(convert(CVODEMemPtr,cvode_mem),convert(Cint,maxncf))
end

function __CVodeSetNonlinConvCoef(cvode_mem::CVODEMemPtr,nlscoef::realtype)
    ccall((:CVodeSetNonlinConvCoef,libsundials_cvode),Cint,(CVODEMemPtr,realtype),cvode_mem,nlscoef)
end

function CVodeSetNonlinConvCoef(cvode_mem,nlscoef)
    __CVodeSetNonlinConvCoef(convert(CVODEMemPtr,cvode_mem),nlscoef)
end

function __CVodeSetIterType(cvode_mem::CVODEMemPtr,iter::Cint)
    ccall((:CVodeSetIterType,libsundials_cvode),Cint,(CVODEMemPtr,Cint),cvode_mem,iter)
end

function CVodeSetIterType(cvode_mem,iter)
    __CVodeSetIterType(convert(CVODEMemPtr,cvode_mem),convert(Cint,iter))
end

function __CVodeSetRootDirection(cvode_mem::CVODEMemPtr,rootdir::Ptr{Cint})
    ccall((:CVodeSetRootDirection,libsundials_cvode),Cint,(CVODEMemPtr,Ptr{Cint}),cvode_mem,rootdir)
end

function CVodeSetRootDirection(cvode_mem,rootdir)
    __CVodeSetRootDirection(convert(CVODEMemPtr,cvode_mem),pointer(rootdir))
end

function __CVodeSetNoInactiveRootWarn(cvode_mem::CVODEMemPtr)
    ccall((:CVodeSetNoInactiveRootWarn,libsundials_cvode),Cint,(CVODEMemPtr,),cvode_mem)
end

function CVodeSetNoInactiveRootWarn(cvode_mem)
    __CVodeSetNoInactiveRootWarn(convert(CVODEMemPtr,cvode_mem))
end

function __CVodeInit(cvode_mem::CVODEMemPtr,f::CVRhsFn,t0::realtype,y0::N_Vector)
    ccall((:CVodeInit,libsundials_cvode),Cint,(CVODEMemPtr,CVRhsFn,realtype,N_Vector),cvode_mem,f,t0,y0)
end

function CVodeInit(cvode_mem,f,t0,y0)
    __y0 = convert(NVector,y0)
    __CVodeInit(convert(CVODEMemPtr,cvode_mem),CVRhsFn_wrapper(f),t0,convert(N_Vector,__y0))
end

function __CVodeReInit(cvode_mem::CVODEMemPtr,t0::realtype,y0::N_Vector)
    ccall((:CVodeReInit,libsundials_cvode),Cint,(CVODEMemPtr,realtype,N_Vector),cvode_mem,t0,y0)
end

function CVodeReInit(cvode_mem,t0,y0)
    __y0 = convert(NVector,y0)
    __CVodeReInit(convert(CVODEMemPtr,cvode_mem),t0,convert(N_Vector,__y0))
end

function __CVodeSStolerances(cvode_mem::CVODEMemPtr,reltol::realtype,abstol::realtype)
    ccall((:CVodeSStolerances,libsundials_cvode),Cint,(CVODEMemPtr,realtype,realtype),cvode_mem,reltol,abstol)
end

function CVodeSStolerances(cvode_mem,reltol,abstol)
    __CVodeSStolerances(convert(CVODEMemPtr,cvode_mem),reltol,abstol)
end

function __CVodeSVtolerances(cvode_mem::CVODEMemPtr,reltol::realtype,abstol::N_Vector)
    ccall((:CVodeSVtolerances,libsundials_cvode),Cint,(CVODEMemPtr,realtype,N_Vector),cvode_mem,reltol,abstol)
end

function CVodeSVtolerances(cvode_mem,reltol,abstol)
    __abstol = convert(NVector,abstol)
    __CVodeSVtolerances(convert(CVODEMemPtr,cvode_mem),reltol,convert(N_Vector,__abstol))
end

function __CVodeWFtolerances(cvode_mem::CVODEMemPtr,efun::CVEwtFn)
    ccall((:CVodeWFtolerances,libsundials_cvode),Cint,(CVODEMemPtr,CVEwtFn),cvode_mem,efun)
end

function CVodeWFtolerances(cvode_mem,efun)
    __CVodeWFtolerances(convert(CVODEMemPtr,cvode_mem),efun)
end

function __CVodeRootInit(cvode_mem::CVODEMemPtr,nrtfn::Cint,g::CVRootFn)
    ccall((:CVodeRootInit,libsundials_cvode),Cint,(CVODEMemPtr,Cint,CVRootFn),cvode_mem,nrtfn,g)
end

function CVodeRootInit(cvode_mem,nrtfn,g)
    __CVodeRootInit(convert(CVODEMemPtr,cvode_mem),convert(Cint,nrtfn),CVRootFn_wrapper(g))
end

function __CVode(cvode_mem::CVODEMemPtr,tout::realtype,yout::N_Vector,tret::Ptr{realtype},itask::Cint)
    ccall((:CVode,libsundials_cvode),Cint,(CVODEMemPtr,realtype,N_Vector,Ptr{realtype},Cint),cvode_mem,tout,yout,tret,itask)
end

function CVode(cvode_mem,tout,yout,tret,itask)
    __yout = convert(NVector,yout)
    __CVode(convert(CVODEMemPtr,cvode_mem),tout,convert(N_Vector,__yout),pointer(tret),convert(Cint,itask))
end

function __CVodeGetDky(cvode_mem::CVODEMemPtr,t::realtype,k::Cint,dky::N_Vector)
    ccall((:CVodeGetDky,libsundials_cvode),Cint,(CVODEMemPtr,realtype,Cint,N_Vector),cvode_mem,t,k,dky)
end

function CVodeGetDky(cvode_mem,t,k,dky)
    __dky = convert(NVector,dky)
    __CVodeGetDky(convert(CVODEMemPtr,cvode_mem),t,convert(Cint,k),convert(N_Vector,__dky))
end

function __CVodeGetWorkSpace(cvode_mem::CVODEMemPtr,lenrw::Ptr{Clong},leniw::Ptr{Clong})
    ccall((:CVodeGetWorkSpace,libsundials_cvode),Cint,(CVODEMemPtr,Ptr{Clong},Ptr{Clong}),cvode_mem,lenrw,leniw)
end

function CVodeGetWorkSpace(cvode_mem,lenrw,leniw)
    __CVodeGetWorkSpace(convert(CVODEMemPtr,cvode_mem),pointer(lenrw),pointer(leniw))
end

function __CVodeGetNumSteps(cvode_mem::CVODEMemPtr,nsteps::Ptr{Clong})
    ccall((:CVodeGetNumSteps,libsundials_cvode),Cint,(CVODEMemPtr,Ptr{Clong}),cvode_mem,nsteps)
end

function CVodeGetNumSteps(cvode_mem,nsteps)
    __CVodeGetNumSteps(convert(CVODEMemPtr,cvode_mem),pointer(nsteps))
end

function __CVodeGetNumRhsEvals(cvode_mem::CVODEMemPtr,nfevals::Ptr{Clong})
    ccall((:CVodeGetNumRhsEvals,libsundials_cvode),Cint,(CVODEMemPtr,Ptr{Clong}),cvode_mem,nfevals)
end

function CVodeGetNumRhsEvals(cvode_mem,nfevals)
    __CVodeGetNumRhsEvals(convert(CVODEMemPtr,cvode_mem),pointer(nfevals))
end

function __CVodeGetNumLinSolvSetups(cvode_mem::CVODEMemPtr,nlinsetups::Ptr{Clong})
    ccall((:CVodeGetNumLinSolvSetups,libsundials_cvode),Cint,(CVODEMemPtr,Ptr{Clong}),cvode_mem,nlinsetups)
end

function CVodeGetNumLinSolvSetups(cvode_mem,nlinsetups)
    __CVodeGetNumLinSolvSetups(convert(CVODEMemPtr,cvode_mem),pointer(nlinsetups))
end

function __CVodeGetNumErrTestFails(cvode_mem::CVODEMemPtr,netfails::Ptr{Clong})
    ccall((:CVodeGetNumErrTestFails,libsundials_cvode),Cint,(CVODEMemPtr,Ptr{Clong}),cvode_mem,netfails)
end

function CVodeGetNumErrTestFails(cvode_mem,netfails)
    __CVodeGetNumErrTestFails(convert(CVODEMemPtr,cvode_mem),pointer(netfails))
end

function __CVodeGetLastOrder(cvode_mem::CVODEMemPtr,qlast::Ptr{Cint})
    ccall((:CVodeGetLastOrder,libsundials_cvode),Cint,(CVODEMemPtr,Ptr{Cint}),cvode_mem,qlast)
end

function CVodeGetLastOrder(cvode_mem,qlast)
    __CVodeGetLastOrder(convert(CVODEMemPtr,cvode_mem),pointer(qlast))
end

function __CVodeGetCurrentOrder(cvode_mem::CVODEMemPtr,qcur::Ptr{Cint})
    ccall((:CVodeGetCurrentOrder,libsundials_cvode),Cint,(CVODEMemPtr,Ptr{Cint}),cvode_mem,qcur)
end

function CVodeGetCurrentOrder(cvode_mem,qcur)
    __CVodeGetCurrentOrder(convert(CVODEMemPtr,cvode_mem),pointer(qcur))
end

function __CVodeGetNumStabLimOrderReds(cvode_mem::CVODEMemPtr,nslred::Ptr{Clong})
    ccall((:CVodeGetNumStabLimOrderReds,libsundials_cvode),Cint,(CVODEMemPtr,Ptr{Clong}),cvode_mem,nslred)
end

function CVodeGetNumStabLimOrderReds(cvode_mem,nslred)
    __CVodeGetNumStabLimOrderReds(convert(CVODEMemPtr,cvode_mem),pointer(nslred))
end

function __CVodeGetActualInitStep(cvode_mem::CVODEMemPtr,hinused::Ptr{realtype})
    ccall((:CVodeGetActualInitStep,libsundials_cvode),Cint,(CVODEMemPtr,Ptr{realtype}),cvode_mem,hinused)
end

function CVodeGetActualInitStep(cvode_mem,hinused)
    __CVodeGetActualInitStep(convert(CVODEMemPtr,cvode_mem),pointer(hinused))
end

function __CVodeGetLastStep(cvode_mem::CVODEMemPtr,hlast::Ptr{realtype})
    ccall((:CVodeGetLastStep,libsundials_cvode),Cint,(CVODEMemPtr,Ptr{realtype}),cvode_mem,hlast)
end

function CVodeGetLastStep(cvode_mem,hlast)
    __CVodeGetLastStep(convert(CVODEMemPtr,cvode_mem),pointer(hlast))
end

function __CVodeGetCurrentStep(cvode_mem::CVODEMemPtr,hcur::Ptr{realtype})
    ccall((:CVodeGetCurrentStep,libsundials_cvode),Cint,(CVODEMemPtr,Ptr{realtype}),cvode_mem,hcur)
end

function CVodeGetCurrentStep(cvode_mem,hcur)
    __CVodeGetCurrentStep(convert(CVODEMemPtr,cvode_mem),pointer(hcur))
end

function __CVodeGetCurrentTime(cvode_mem::CVODEMemPtr,tcur::Ptr{realtype})
    ccall((:CVodeGetCurrentTime,libsundials_cvode),Cint,(CVODEMemPtr,Ptr{realtype}),cvode_mem,tcur)
end

function CVodeGetCurrentTime(cvode_mem,tcur)
    __CVodeGetCurrentTime(convert(CVODEMemPtr,cvode_mem),pointer(tcur))
end

function __CVodeGetTolScaleFactor(cvode_mem::CVODEMemPtr,tolsfac::Ptr{realtype})
    ccall((:CVodeGetTolScaleFactor,libsundials_cvode),Cint,(CVODEMemPtr,Ptr{realtype}),cvode_mem,tolsfac)
end

function CVodeGetTolScaleFactor(cvode_mem,tolsfac)
    __CVodeGetTolScaleFactor(convert(CVODEMemPtr,cvode_mem),pointer(tolsfac))
end

function __CVodeGetErrWeights(cvode_mem::CVODEMemPtr,eweight::N_Vector)
    ccall((:CVodeGetErrWeights,libsundials_cvode),Cint,(CVODEMemPtr,N_Vector),cvode_mem,eweight)
end

function CVodeGetErrWeights(cvode_mem,eweight)
    __eweight = convert(NVector,eweight)
    __CVodeGetErrWeights(convert(CVODEMemPtr,cvode_mem),convert(N_Vector,__eweight))
end

function __CVodeGetEstLocalErrors(cvode_mem::CVODEMemPtr,ele::N_Vector)
    ccall((:CVodeGetEstLocalErrors,libsundials_cvode),Cint,(CVODEMemPtr,N_Vector),cvode_mem,ele)
end

function CVodeGetEstLocalErrors(cvode_mem,ele)
    __ele = convert(NVector,ele)
    __CVodeGetEstLocalErrors(convert(CVODEMemPtr,cvode_mem),convert(N_Vector,__ele))
end

function __CVodeGetNumGEvals(cvode_mem::CVODEMemPtr,ngevals::Ptr{Clong})
    ccall((:CVodeGetNumGEvals,libsundials_cvode),Cint,(CVODEMemPtr,Ptr{Clong}),cvode_mem,ngevals)
end

function CVodeGetNumGEvals(cvode_mem,ngevals)
    __CVodeGetNumGEvals(convert(CVODEMemPtr,cvode_mem),pointer(ngevals))
end

function __CVodeGetRootInfo(cvode_mem::CVODEMemPtr,rootsfound::Ptr{Cint})
    ccall((:CVodeGetRootInfo,libsundials_cvode),Cint,(CVODEMemPtr,Ptr{Cint}),cvode_mem,rootsfound)
end

function CVodeGetRootInfo(cvode_mem,rootsfound)
    __CVodeGetRootInfo(convert(CVODEMemPtr,cvode_mem),pointer(rootsfound))
end

function __CVodeGetIntegratorStats(cvode_mem::CVODEMemPtr,nsteps::Ptr{Clong},nfevals::Ptr{Clong},nlinsetups::Ptr{Clong},netfails::Ptr{Clong},qlast::Ptr{Cint},qcur::Ptr{Cint},hinused::Ptr{realtype},hlast::Ptr{realtype},hcur::Ptr{realtype},tcur::Ptr{realtype})
    ccall((:CVodeGetIntegratorStats,libsundials_cvode),Cint,(CVODEMemPtr,Ptr{Clong},Ptr{Clong},Ptr{Clong},Ptr{Clong},Ptr{Cint},Ptr{Cint},Ptr{realtype},Ptr{realtype},Ptr{realtype},Ptr{realtype}),cvode_mem,nsteps,nfevals,nlinsetups,netfails,qlast,qcur,hinused,hlast,hcur,tcur)
end

function CVodeGetIntegratorStats(cvode_mem,nsteps,nfevals,nlinsetups,netfails,qlast,qcur,hinused,hlast,hcur,tcur)
    __CVodeGetIntegratorStats(convert(CVODEMemPtr,cvode_mem),pointer(nsteps),pointer(nfevals),pointer(nlinsetups),pointer(netfails),pointer(qlast),pointer(qcur),pointer(hinused),pointer(hlast),pointer(hcur),pointer(tcur))
end

function __CVodeGetNumNonlinSolvIters(cvode_mem::CVODEMemPtr,nniters::Ptr{Clong})
    ccall((:CVodeGetNumNonlinSolvIters,libsundials_cvode),Cint,(CVODEMemPtr,Ptr{Clong}),cvode_mem,nniters)
end

function CVodeGetNumNonlinSolvIters(cvode_mem,nniters)
    __CVodeGetNumNonlinSolvIters(convert(CVODEMemPtr,cvode_mem),pointer(nniters))
end

function __CVodeGetNumNonlinSolvConvFails(cvode_mem::CVODEMemPtr,nncfails::Ptr{Clong})
    ccall((:CVodeGetNumNonlinSolvConvFails,libsundials_cvode),Cint,(CVODEMemPtr,Ptr{Clong}),cvode_mem,nncfails)
end

function CVodeGetNumNonlinSolvConvFails(cvode_mem,nncfails)
    __CVodeGetNumNonlinSolvConvFails(convert(CVODEMemPtr,cvode_mem),pointer(nncfails))
end

function __CVodeGetNonlinSolvStats(cvode_mem::CVODEMemPtr,nniters::Ptr{Clong},nncfails::Ptr{Clong})
    ccall((:CVodeGetNonlinSolvStats,libsundials_cvode),Cint,(CVODEMemPtr,Ptr{Clong},Ptr{Clong}),cvode_mem,nniters,nncfails)
end

function CVodeGetNonlinSolvStats(cvode_mem,nniters,nncfails)
    __CVodeGetNonlinSolvStats(convert(CVODEMemPtr,cvode_mem),pointer(nniters),pointer(nncfails))
end

function __CVodeGetReturnFlagName(flag::Clong)
    ccall((:CVodeGetReturnFlagName,libsundials_cvode),Ptr{UInt8},(Clong,),flag)
end

function CVodeGetReturnFlagName(flag)
    __CVodeGetReturnFlagName(convert(Clong,flag))
end

function CVodeFree(cvode_mem::Ref{CVODEMemPtr})
    ccall((:CVodeFree,libsundials_cvode),Void,(Ref{CVODEMemPtr},),cvode_mem)
end
# Julia wrapper for header: /home/astukalov/.julia/v0.5/Sundials/deps/usr/include/cvode/cvode_direct.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function __CVDlsSetDenseJacFn(cvode_mem::CVODEMemPtr,jac::CVDlsDenseJacFn)
    ccall((:CVDlsSetDenseJacFn,libsundials_cvode),Cint,(CVODEMemPtr,CVDlsDenseJacFn),cvode_mem,jac)
end

function CVDlsSetDenseJacFn(cvode_mem,jac)
    __CVDlsSetDenseJacFn(convert(CVODEMemPtr,cvode_mem),jac)
end

function __CVDlsSetBandJacFn(cvode_mem::CVODEMemPtr,jac::CVDlsBandJacFn)
    ccall((:CVDlsSetBandJacFn,libsundials_cvode),Cint,(CVODEMemPtr,CVDlsBandJacFn),cvode_mem,jac)
end

function CVDlsSetBandJacFn(cvode_mem,jac)
    __CVDlsSetBandJacFn(convert(CVODEMemPtr,cvode_mem),jac)
end

function __CVDlsGetWorkSpace(cvode_mem::CVODEMemPtr,lenrwLS::Ptr{Clong},leniwLS::Ptr{Clong})
    ccall((:CVDlsGetWorkSpace,libsundials_cvode),Cint,(CVODEMemPtr,Ptr{Clong},Ptr{Clong}),cvode_mem,lenrwLS,leniwLS)
end

function CVDlsGetWorkSpace(cvode_mem,lenrwLS,leniwLS)
    __CVDlsGetWorkSpace(convert(CVODEMemPtr,cvode_mem),pointer(lenrwLS),pointer(leniwLS))
end

function __CVDlsGetNumJacEvals(cvode_mem::CVODEMemPtr,njevals::Ptr{Clong})
    ccall((:CVDlsGetNumJacEvals,libsundials_cvode),Cint,(CVODEMemPtr,Ptr{Clong}),cvode_mem,njevals)
end

function CVDlsGetNumJacEvals(cvode_mem,njevals)
    __CVDlsGetNumJacEvals(convert(CVODEMemPtr,cvode_mem),pointer(njevals))
end

function __CVDlsGetNumRhsEvals(cvode_mem::CVODEMemPtr,nfevalsLS::Ptr{Clong})
    ccall((:CVDlsGetNumRhsEvals,libsundials_cvode),Cint,(CVODEMemPtr,Ptr{Clong}),cvode_mem,nfevalsLS)
end

function CVDlsGetNumRhsEvals(cvode_mem,nfevalsLS)
    __CVDlsGetNumRhsEvals(convert(CVODEMemPtr,cvode_mem),pointer(nfevalsLS))
end

function __CVDlsGetLastFlag(cvode_mem::CVODEMemPtr,flag::Ptr{Clong})
    ccall((:CVDlsGetLastFlag,libsundials_cvode),Cint,(CVODEMemPtr,Ptr{Clong}),cvode_mem,flag)
end

function CVDlsGetLastFlag(cvode_mem,flag)
    __CVDlsGetLastFlag(convert(CVODEMemPtr,cvode_mem),pointer(flag))
end

function __CVDlsGetReturnFlagName(flag::Clong)
    ccall((:CVDlsGetReturnFlagName,libsundials_cvode),Ptr{UInt8},(Clong,),flag)
end

function CVDlsGetReturnFlagName(flag)
    __CVDlsGetReturnFlagName(convert(Clong,flag))
end
# Julia wrapper for header: /home/astukalov/.julia/v0.5/Sundials/deps/usr/include/cvode/cvode_spils.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function __CVSpilsSetPrecType(cvode_mem::CVODEMemPtr,pretype::Cint)
    ccall((:CVSpilsSetPrecType,libsundials_cvode),Cint,(CVODEMemPtr,Cint),cvode_mem,pretype)
end

function CVSpilsSetPrecType(cvode_mem,pretype)
    __CVSpilsSetPrecType(convert(CVODEMemPtr,cvode_mem),convert(Cint,pretype))
end

function __CVSpilsSetGSType(cvode_mem::CVODEMemPtr,gstype::Cint)
    ccall((:CVSpilsSetGSType,libsundials_cvode),Cint,(CVODEMemPtr,Cint),cvode_mem,gstype)
end

function CVSpilsSetGSType(cvode_mem,gstype)
    __CVSpilsSetGSType(convert(CVODEMemPtr,cvode_mem),convert(Cint,gstype))
end

function __CVSpilsSetMaxl(cvode_mem::CVODEMemPtr,maxl::Cint)
    ccall((:CVSpilsSetMaxl,libsundials_cvode),Cint,(CVODEMemPtr,Cint),cvode_mem,maxl)
end

function CVSpilsSetMaxl(cvode_mem,maxl)
    __CVSpilsSetMaxl(convert(CVODEMemPtr,cvode_mem),convert(Cint,maxl))
end

function __CVSpilsSetEpsLin(cvode_mem::CVODEMemPtr,eplifac::realtype)
    ccall((:CVSpilsSetEpsLin,libsundials_cvode),Cint,(CVODEMemPtr,realtype),cvode_mem,eplifac)
end

function CVSpilsSetEpsLin(cvode_mem,eplifac)
    __CVSpilsSetEpsLin(convert(CVODEMemPtr,cvode_mem),eplifac)
end

function __CVSpilsSetPreconditioner(cvode_mem::CVODEMemPtr,pset::CVSpilsPrecSetupFn,psolve::CVSpilsPrecSolveFn)
    ccall((:CVSpilsSetPreconditioner,libsundials_cvode),Cint,(CVODEMemPtr,CVSpilsPrecSetupFn,CVSpilsPrecSolveFn),cvode_mem,pset,psolve)
end

function CVSpilsSetPreconditioner(cvode_mem,pset,psolve)
    __CVSpilsSetPreconditioner(convert(CVODEMemPtr,cvode_mem),pset,psolve)
end

function __CVSpilsSetJacTimesVecFn(cvode_mem::CVODEMemPtr,jtv::CVSpilsJacTimesVecFn)
    ccall((:CVSpilsSetJacTimesVecFn,libsundials_cvode),Cint,(CVODEMemPtr,CVSpilsJacTimesVecFn),cvode_mem,jtv)
end

function CVSpilsSetJacTimesVecFn(cvode_mem,jtv)
    __CVSpilsSetJacTimesVecFn(convert(CVODEMemPtr,cvode_mem),jtv)
end

function __CVSpilsGetWorkSpace(cvode_mem::CVODEMemPtr,lenrwLS::Ptr{Clong},leniwLS::Ptr{Clong})
    ccall((:CVSpilsGetWorkSpace,libsundials_cvode),Cint,(CVODEMemPtr,Ptr{Clong},Ptr{Clong}),cvode_mem,lenrwLS,leniwLS)
end

function CVSpilsGetWorkSpace(cvode_mem,lenrwLS,leniwLS)
    __CVSpilsGetWorkSpace(convert(CVODEMemPtr,cvode_mem),pointer(lenrwLS),pointer(leniwLS))
end

function __CVSpilsGetNumPrecEvals(cvode_mem::CVODEMemPtr,npevals::Ptr{Clong})
    ccall((:CVSpilsGetNumPrecEvals,libsundials_cvode),Cint,(CVODEMemPtr,Ptr{Clong}),cvode_mem,npevals)
end

function CVSpilsGetNumPrecEvals(cvode_mem,npevals)
    __CVSpilsGetNumPrecEvals(convert(CVODEMemPtr,cvode_mem),pointer(npevals))
end

function __CVSpilsGetNumPrecSolves(cvode_mem::CVODEMemPtr,npsolves::Ptr{Clong})
    ccall((:CVSpilsGetNumPrecSolves,libsundials_cvode),Cint,(CVODEMemPtr,Ptr{Clong}),cvode_mem,npsolves)
end

function CVSpilsGetNumPrecSolves(cvode_mem,npsolves)
    __CVSpilsGetNumPrecSolves(convert(CVODEMemPtr,cvode_mem),pointer(npsolves))
end

function __CVSpilsGetNumLinIters(cvode_mem::CVODEMemPtr,nliters::Ptr{Clong})
    ccall((:CVSpilsGetNumLinIters,libsundials_cvode),Cint,(CVODEMemPtr,Ptr{Clong}),cvode_mem,nliters)
end

function CVSpilsGetNumLinIters(cvode_mem,nliters)
    __CVSpilsGetNumLinIters(convert(CVODEMemPtr,cvode_mem),pointer(nliters))
end

function __CVSpilsGetNumConvFails(cvode_mem::CVODEMemPtr,nlcfails::Ptr{Clong})
    ccall((:CVSpilsGetNumConvFails,libsundials_cvode),Cint,(CVODEMemPtr,Ptr{Clong}),cvode_mem,nlcfails)
end

function CVSpilsGetNumConvFails(cvode_mem,nlcfails)
    __CVSpilsGetNumConvFails(convert(CVODEMemPtr,cvode_mem),pointer(nlcfails))
end

function __CVSpilsGetNumJtimesEvals(cvode_mem::CVODEMemPtr,njvevals::Ptr{Clong})
    ccall((:CVSpilsGetNumJtimesEvals,libsundials_cvode),Cint,(CVODEMemPtr,Ptr{Clong}),cvode_mem,njvevals)
end

function CVSpilsGetNumJtimesEvals(cvode_mem,njvevals)
    __CVSpilsGetNumJtimesEvals(convert(CVODEMemPtr,cvode_mem),pointer(njvevals))
end

function __CVSpilsGetNumRhsEvals(cvode_mem::CVODEMemPtr,nfevalsLS::Ptr{Clong})
    ccall((:CVSpilsGetNumRhsEvals,libsundials_cvode),Cint,(CVODEMemPtr,Ptr{Clong}),cvode_mem,nfevalsLS)
end

function CVSpilsGetNumRhsEvals(cvode_mem,nfevalsLS)
    __CVSpilsGetNumRhsEvals(convert(CVODEMemPtr,cvode_mem),pointer(nfevalsLS))
end

function __CVSpilsGetLastFlag(cvode_mem::CVODEMemPtr,flag::Ptr{Clong})
    ccall((:CVSpilsGetLastFlag,libsundials_cvode),Cint,(CVODEMemPtr,Ptr{Clong}),cvode_mem,flag)
end

function CVSpilsGetLastFlag(cvode_mem,flag)
    __CVSpilsGetLastFlag(convert(CVODEMemPtr,cvode_mem),pointer(flag))
end

function __CVSpilsGetReturnFlagName(flag::Clong)
    ccall((:CVSpilsGetReturnFlagName,libsundials_cvode),Ptr{UInt8},(Clong,),flag)
end

function CVSpilsGetReturnFlagName(flag)
    __CVSpilsGetReturnFlagName(convert(Clong,flag))
end
# Julia wrapper for header: /home/astukalov/.julia/v0.5/Sundials/deps/usr/include/cvode/cvode_band.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function __CVBand(cvode_mem::CVODEMemPtr,N::Clong,mupper::Clong,mlower::Clong)
    ccall((:CVBand,libsundials_cvode),Cint,(CVODEMemPtr,Clong,Clong,Clong),cvode_mem,N,mupper,mlower)
end

function CVBand(cvode_mem,N,mupper,mlower)
    __CVBand(convert(CVODEMemPtr,cvode_mem),convert(Clong,N),convert(Clong,mupper),convert(Clong,mlower))
end
# Julia wrapper for header: /home/astukalov/.julia/v0.5/Sundials/deps/usr/include/cvode/cvode_bandpre.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function __CVBandPrecInit(cvode_mem::CVODEMemPtr,N::Clong,mu::Clong,ml::Clong)
    ccall((:CVBandPrecInit,libsundials_cvode),Cint,(CVODEMemPtr,Clong,Clong,Clong),cvode_mem,N,mu,ml)
end

function CVBandPrecInit(cvode_mem,N,mu,ml)
    __CVBandPrecInit(convert(CVODEMemPtr,cvode_mem),convert(Clong,N),convert(Clong,mu),convert(Clong,ml))
end

function __CVBandPrecGetWorkSpace(cvode_mem::CVODEMemPtr,lenrwLS::Ptr{Clong},leniwLS::Ptr{Clong})
    ccall((:CVBandPrecGetWorkSpace,libsundials_cvode),Cint,(CVODEMemPtr,Ptr{Clong},Ptr{Clong}),cvode_mem,lenrwLS,leniwLS)
end

function CVBandPrecGetWorkSpace(cvode_mem,lenrwLS,leniwLS)
    __CVBandPrecGetWorkSpace(convert(CVODEMemPtr,cvode_mem),pointer(lenrwLS),pointer(leniwLS))
end

function __CVBandPrecGetNumRhsEvals(cvode_mem::CVODEMemPtr,nfevalsBP::Ptr{Clong})
    ccall((:CVBandPrecGetNumRhsEvals,libsundials_cvode),Cint,(CVODEMemPtr,Ptr{Clong}),cvode_mem,nfevalsBP)
end

function CVBandPrecGetNumRhsEvals(cvode_mem,nfevalsBP)
    __CVBandPrecGetNumRhsEvals(convert(CVODEMemPtr,cvode_mem),pointer(nfevalsBP))
end
# Julia wrapper for header: /home/astukalov/.julia/v0.5/Sundials/deps/usr/include/cvode/cvode_bbdpre.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function __CVBBDPrecInit(cvode_mem::CVODEMemPtr,Nlocal::Clong,mudq::Clong,mldq::Clong,mukeep::Clong,mlkeep::Clong,dqrely::realtype,gloc::CVLocalFn,cfn::CVCommFn)
    ccall((:CVBBDPrecInit,libsundials_cvode),Cint,(CVODEMemPtr,Clong,Clong,Clong,Clong,Clong,realtype,CVLocalFn,CVCommFn),cvode_mem,Nlocal,mudq,mldq,mukeep,mlkeep,dqrely,gloc,cfn)
end

function CVBBDPrecInit(cvode_mem,Nlocal,mudq,mldq,mukeep,mlkeep,dqrely,gloc,cfn)
    __CVBBDPrecInit(convert(CVODEMemPtr,cvode_mem),convert(Clong,Nlocal),convert(Clong,mudq),convert(Clong,mldq),convert(Clong,mukeep),convert(Clong,mlkeep),dqrely,gloc,cfn)
end

function __CVBBDPrecReInit(cvode_mem::CVODEMemPtr,mudq::Clong,mldq::Clong,dqrely::realtype)
    ccall((:CVBBDPrecReInit,libsundials_cvode),Cint,(CVODEMemPtr,Clong,Clong,realtype),cvode_mem,mudq,mldq,dqrely)
end

function CVBBDPrecReInit(cvode_mem,mudq,mldq,dqrely)
    __CVBBDPrecReInit(convert(CVODEMemPtr,cvode_mem),convert(Clong,mudq),convert(Clong,mldq),dqrely)
end

function __CVBBDPrecGetWorkSpace(cvode_mem::CVODEMemPtr,lenrwLS::Ptr{Clong},leniwLS::Ptr{Clong})
    ccall((:CVBBDPrecGetWorkSpace,libsundials_cvode),Cint,(CVODEMemPtr,Ptr{Clong},Ptr{Clong}),cvode_mem,lenrwLS,leniwLS)
end

function CVBBDPrecGetWorkSpace(cvode_mem,lenrwLS,leniwLS)
    __CVBBDPrecGetWorkSpace(convert(CVODEMemPtr,cvode_mem),pointer(lenrwLS),pointer(leniwLS))
end

function __CVBBDPrecGetNumGfnEvals(cvode_mem::CVODEMemPtr,ngevalsBBDP::Ptr{Clong})
    ccall((:CVBBDPrecGetNumGfnEvals,libsundials_cvode),Cint,(CVODEMemPtr,Ptr{Clong}),cvode_mem,ngevalsBBDP)
end

function CVBBDPrecGetNumGfnEvals(cvode_mem,ngevalsBBDP)
    __CVBBDPrecGetNumGfnEvals(convert(CVODEMemPtr,cvode_mem),pointer(ngevalsBBDP))
end
# Julia wrapper for header: /home/astukalov/.julia/v0.5/Sundials/deps/usr/include/cvode/cvode_dense.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function __CVDense(cvode_mem::CVODEMemPtr,N::Clong)
    ccall((:CVDense,libsundials_cvode),Cint,(CVODEMemPtr,Clong),cvode_mem,N)
end

function CVDense(cvode_mem,N)
    __CVDense(convert(CVODEMemPtr,cvode_mem),convert(Clong,N))
end
# Julia wrapper for header: /home/astukalov/.julia/v0.5/Sundials/deps/usr/include/cvode/cvode_diag.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function __CVDiag(cvode_mem::CVODEMemPtr)
    ccall((:CVDiag,libsundials_cvode),Cint,(CVODEMemPtr,),cvode_mem)
end

function CVDiag(cvode_mem)
    __CVDiag(convert(CVODEMemPtr,cvode_mem))
end

function __CVDiagGetWorkSpace(cvode_mem::CVODEMemPtr,lenrwLS::Ptr{Clong},leniwLS::Ptr{Clong})
    ccall((:CVDiagGetWorkSpace,libsundials_cvode),Cint,(CVODEMemPtr,Ptr{Clong},Ptr{Clong}),cvode_mem,lenrwLS,leniwLS)
end

function CVDiagGetWorkSpace(cvode_mem,lenrwLS,leniwLS)
    __CVDiagGetWorkSpace(convert(CVODEMemPtr,cvode_mem),pointer(lenrwLS),pointer(leniwLS))
end

function __CVDiagGetNumRhsEvals(cvode_mem::CVODEMemPtr,nfevalsLS::Ptr{Clong})
    ccall((:CVDiagGetNumRhsEvals,libsundials_cvode),Cint,(CVODEMemPtr,Ptr{Clong}),cvode_mem,nfevalsLS)
end

function CVDiagGetNumRhsEvals(cvode_mem,nfevalsLS)
    __CVDiagGetNumRhsEvals(convert(CVODEMemPtr,cvode_mem),pointer(nfevalsLS))
end

function __CVDiagGetLastFlag(cvode_mem::CVODEMemPtr,flag::Ptr{Clong})
    ccall((:CVDiagGetLastFlag,libsundials_cvode),Cint,(CVODEMemPtr,Ptr{Clong}),cvode_mem,flag)
end

function CVDiagGetLastFlag(cvode_mem,flag)
    __CVDiagGetLastFlag(convert(CVODEMemPtr,cvode_mem),pointer(flag))
end

function __CVDiagGetReturnFlagName(flag::Clong)
    ccall((:CVDiagGetReturnFlagName,libsundials_cvode),Ptr{UInt8},(Clong,),flag)
end

function CVDiagGetReturnFlagName(flag)
    __CVDiagGetReturnFlagName(convert(Clong,flag))
end
# Julia wrapper for header: /home/astukalov/.julia/v0.5/Sundials/deps/usr/include/cvode/cvode_impl.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

# Julia wrapper for header: /home/astukalov/.julia/v0.5/Sundials/deps/usr/include/cvode/cvode_spbcgs.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function __CVSpbcg(cvode_mem::CVODEMemPtr,pretype::Cint,maxl::Cint)
    ccall((:CVSpbcg,libsundials_cvode),Cint,(CVODEMemPtr,Cint,Cint),cvode_mem,pretype,maxl)
end

function CVSpbcg(cvode_mem,pretype,maxl)
    __CVSpbcg(convert(CVODEMemPtr,cvode_mem),convert(Cint,pretype),convert(Cint,maxl))
end
# Julia wrapper for header: /home/astukalov/.julia/v0.5/Sundials/deps/usr/include/cvode/cvode_spgmr.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function __CVSpgmr(cvode_mem::CVODEMemPtr,pretype::Cint,maxl::Cint)
    ccall((:CVSpgmr,libsundials_cvode),Cint,(CVODEMemPtr,Cint,Cint),cvode_mem,pretype,maxl)
end

function CVSpgmr(cvode_mem,pretype,maxl)
    __CVSpgmr(convert(CVODEMemPtr,cvode_mem),convert(Cint,pretype),convert(Cint,maxl))
end
# Julia wrapper for header: /home/astukalov/.julia/v0.5/Sundials/deps/usr/include/cvode/cvode_sptfqmr.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function __CVSptfqmr(cvode_mem::CVODEMemPtr,pretype::Cint,maxl::Cint)
    ccall((:CVSptfqmr,libsundials_cvode),Cint,(CVODEMemPtr,Cint,Cint),cvode_mem,pretype,maxl)
end

function CVSptfqmr(cvode_mem,pretype,maxl)
    __CVSptfqmr(convert(CVODEMemPtr,cvode_mem),convert(Cint,pretype),convert(Cint,maxl))
end


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

mutable struct VirtualMachineScaleSetReimageParameters <: SwaggerModel
    tempDisk::Any # spec type: Union{ Nothing, Bool } # spec name: tempDisk

    function VirtualMachineScaleSetReimageParameters(;tempDisk=nothing)
        o = new()
        validate_property(VirtualMachineScaleSetReimageParameters, Symbol("tempDisk"), tempDisk)
        setfield!(o, Symbol("tempDisk"), tempDisk)
        o
    end
end # type VirtualMachineScaleSetReimageParameters

const _property_map_VirtualMachineScaleSetReimageParameters = Dict{Symbol,Symbol}(Symbol("tempDisk")=>Symbol("tempDisk"))
const _property_types_VirtualMachineScaleSetReimageParameters = Dict{Symbol,String}(Symbol("tempDisk")=>"Bool")
Base.propertynames(::Type{ VirtualMachineScaleSetReimageParameters }) = collect(keys(_property_map_VirtualMachineScaleSetReimageParameters))
Swagger.property_type(::Type{ VirtualMachineScaleSetReimageParameters }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_VirtualMachineScaleSetReimageParameters[name]))}
Swagger.field_name(::Type{ VirtualMachineScaleSetReimageParameters }, property_name::Symbol) =  _property_map_VirtualMachineScaleSetReimageParameters[property_name]

function check_required(o::VirtualMachineScaleSetReimageParameters)
    true
end

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


text:
"""
control points of Bézier curve from given function and range
"""
function BézPts(f,t0,t1)
    f(t0),
    3*(f(2*t0/3+t1/3)-(8*f(t0)+f(t1))/27)-3*(f(t0/3+2*t1/3)-(f(t0)+8*f(t1))/27)/2,
    -3*(f(2*t0/3+t1/3)-(8*f(t0)+f(t1))/27)/2+3*(f(t0/3+2*t1/3)-(f(t0)+8*f(t1))/27),
    f(t1)
end

function LxrPt(p::AbstractVector{<:Real},unitlength)
    Point(unitlength*[1,-1].*p...)
end

"""
export svg file
"""
function save_svg(name::String, M::AbstractBSplineManifold; up=5, down=-5, right=5, left=-5, mesh=(10,10), unitlength=100, points=true, thickness=1, backgroundcolor=RGB(1,1,1), maincolor=RGB(1,0,0))
    if split(name,'.')[end] ≠ "svg"
        name = name * ".svg"
    end
    if dim(M) == 1
        _save_luxor_1d2d(name, M, up=up, down=down, right=right, left=left, mesh=mesh, unitlength=unitlength, points=points, thickness=thickness, backgroundcolor=backgroundcolor, maincolor=maincolor)
    elseif dim(M) == 2
        _save_luxor_2d2d(name, M, up=up, down=down, right=right, left=left, mesh=mesh, unitlength=unitlength, points=points, thickness=thickness, backgroundcolor=backgroundcolor, maincolor=maincolor)
    else
        error("the dimension of B-spline manifold must be 2 or less")
    end
end

"""
export png file
"""
function save_png(name::String, M::AbstractBSplineManifold; up=5, down=-5, right=5, left=-5, mesh=(10,10), unitlength=100, points=true, thickness=1, backgroundcolor=RGB(1,1,1), maincolor=RGB(1,0,0))
    if split(name,'.')[end] ≠ "png"
        name = name * ".png"
    end
    if dim(M) == 1
        _save_luxor_1d2d(name, M, up=up, down=down, right=right, left=left, mesh=mesh, unitlength=unitlength, points=points, thickness=thickness, backgroundcolor=backgroundcolor, maincolor=maincolor)
    elseif dim(M) == 2
        _save_luxor_2d2d(name, M, up=up, down=down, right=right, left=left, mesh=mesh, unitlength=unitlength, points=points, thickness=thickness, backgroundcolor=backgroundcolor, maincolor=maincolor)
    else
        error("the dimension of B-spline manifold must be 2 or less")
    end
end

"""
export png file
"""
function save_png(name::String, M::AbstractBSplineManifold, colors::AbstractArray{<:Colorant,2}; up=5, down=-5, right=5, left=-5, unitlength=100)
    if split(name,'.')[end] ≠ "png"
        name = name * ".png"
    end

    if dim(M) == 2
        _save_luxor_2d2d_color(name, M, colors, up=up, down=down, right=right, left=left, unitlength=unitlength)
    else
        error("the dimension of B-spline manifold must be 2")
    end
end

"""
export png file
"""
function save_png(name::String, M::AbstractBSplineManifold, colorfunc::Function; up=5, down=-5, right=5, left=-5, unitlength=100)
    if split(name,'.')[end] ≠ "png"
        name = name * ".png"
    end

    if dim(M) == 2
        _save_luxor_2d2d_color(name, M, colorfunc, up=up, down=down, right=right, left=left, unitlength=unitlength)
    else
        error("the dimension of B-spline manifold must be 2")
    end
end

function _save_luxor_2d2d(name::String, M::AbstractBSplineManifold; up=5, down=-5, right=5, left=-5, mesh=(10,10), unitlength=100, points=true, thickness=1, backgroundcolor=RGB(1,1,1), maincolor=RGB(1,0,0), subcolor=RGB(0.5,0.5,0.5))
    linecolor = maincolor
    fillcolor = weighted_color_mean(0.5, maincolor, colorant"white")
    segmentcolor = subcolor
    pointcolor = weighted_color_mean(0.5, subcolor, colorant"black")

    P1, P2 = P = collect(bsplinespaces(M))
    p¹, p² = p = degree.(P)
    k¹, k² = k = knots.(P)
    a = controlpoints(M)
    n¹, n² = n = dim.(P)

    K¹, K² = K = [unique(k[i][1+p[i]:end-p[i]]) for i in 1:2]
    N¹, N² = length.(K).-1
    m¹, m² = mesh

    Drawing(unitlength*(right-left),unitlength*(up-down),name)
    Luxor.origin(-unitlength*left,unitlength*up)
    setline(thickness)
    background(backgroundcolor)

    setcolor(fillcolor)
    drawbezierpath(BezierPath(vcat(
        [BezierPathSegment(map(p->LxrPt(p,unitlength),BézPts(u¹->M([u¹,K²[1]]),K¹[i],K¹[i+1]))...) for i in 1:N¹],
        [BezierPathSegment(map(p->LxrPt(p,unitlength),BézPts(u²->M([K¹[end],u²]),K²[i],K²[i+1]))...) for i in 1:N²],
        [BezierPathSegment(map(p->LxrPt(p,unitlength),BézPts(u¹->M([u¹,K²[end]]),K¹[end-i+1],K¹[end-i]))...) for i in 1:N¹],
        [BezierPathSegment(map(p->LxrPt(p,unitlength),BézPts(u²->M([K¹[1],u²]),K²[end-i+1],K²[end-i]))...) for i in 1:N²]
    )),:fill,close=true)

    setcolor(linecolor)
    for u¹ in range(K¹[1],stop=K¹[end],length=m¹+1)
        drawbezierpath(BezierPath([BezierPathSegment(map(p->LxrPt(p,unitlength),BézPts(u²->M([u¹,u²]),K²[i],K²[i+1]))...) for i in 1:N²]),:stroke)
    end
    for u² in range(K²[1],stop=K²[end],length=m²+1)
        drawbezierpath(BezierPath([BezierPathSegment(map(p->LxrPt(p,unitlength),BézPts(u¹->M([u¹,u²]),K¹[i],K¹[i+1]))...) for i in 1:N¹]),:stroke)
    end

    if points
        CtrlPts = [LxrPt(a[i,j,:],unitlength) for i in 1:size(a)[1], j in 1:size(a)[2]]

        setcolor(segmentcolor)
        setline(thickness)
        for i in 1:n¹
            poly(CtrlPts[i,:], :stroke)
        end
        for j in 1:n²
            poly(CtrlPts[:,j], :stroke)
        end

        setcolor(pointcolor)
        map(p->circle(p,3*thickness,:fill), CtrlPts)
    end
    finish()
    return nothing
end

function _save_luxor_1d2d(name::String, M::AbstractBSplineManifold; up=5, down=-5, right=5, left=-5, mesh=10, unitlength=100, points=true, thickness=1, backgroundcolor=RGB(1,1,1), maincolor=RGB(1,0,0), subcolor=RGB(.5,.5,.5))
    linecolor = maincolor
    segmentcolor = subcolor
    pointcolor = weighted_color_mean(0.5, subcolor, colorant"black")

    P1, = P = collect(bsplinespaces(M))
    p¹, = p = degree.(P)
    k¹, = k = knots.(P)
    a = controlpoints(M)
    n¹, = n = dim.(P)

    K¹, = K = [unique(k[i][1+p[i]:end-p[i]]) for i in 1:1]
    N¹, = length.(K).-1
    m¹, = mesh

    Drawing(unitlength*(right-left),unitlength*(up-down),name)
    Luxor.origin(-unitlength*left,unitlength*up)
    setline(2*thickness)
    background(backgroundcolor)

    setcolor(linecolor)
    drawbezierpath(BezierPath([BezierPathSegment(map(p->LxrPt(p,unitlength),BézPts(u¹->M([u¹]),K¹[i],K¹[i+1]))...) for i in 1:N¹]),:stroke)

    if points
        CtrlPts = [LxrPt(a[i,:],unitlength) for i in 1:size(a)[1]]

        setcolor(segmentcolor)
        setline(thickness)
        poly(CtrlPts[:], :stroke)

        setcolor(pointcolor)
        map(p->circle(p,3*thickness,:fill), CtrlPts)
    end
    finish()
    return nothing
end

function _save_luxor_2d2d_color(name::String, M::AbstractBSplineManifold, colors::AbstractArray{<:Colorant,2}; up=5, down=-5, right=5, left=-5, unitlength=100)
    P = collect(bsplinespaces(M))
    colorfunc(u) = sum(bsplinebasis(P,u).*colors)
    _save_luxor_2d2d_color(name, M, colorfunc; up=up, down=down, right=right, left=left, unitlength=unitlength)
end

function _save_luxor_2d2d_color(name::String, M::AbstractBSplineManifold, colorfunc::Function; up=5, down=-5, right=5, left=-5, unitlength=100)
    mesh = 10

    P = collect(bsplinespaces(M))
    p¹, p² = p = degree.(P)
    k¹, k² = k = knots.(P)
    a = controlpoints(M)
    n¹, n² = n = dim.(P)

    D = [k[i][1+p[i]]..k[i][end-p[i]] for i in 1:2]

    K¹ = unique(vcat([collect(range(k¹[i], k¹[i+1], length=mesh+1)) for i in 1+p¹:length(k¹)-p¹-1]...))
    K² = unique(vcat([collect(range(k²[i], k²[i+1], length=mesh+1)) for i in 1+p²:length(k²)-p²-1]...))

    Drawing(unitlength*(right-left),unitlength*(up-down),name)
    Luxor.origin(-unitlength*left,unitlength*up)
    background(RGBA(0,0,0,0))

    for I₁ in 1:length(K¹)-1, I₂ in 1:length(K²)-1
        BézPth=BezierPath([
                BezierPathSegment(map(p->LxrPt(p,unitlength),BézPts(t->M([t,K²[I₂]]),K¹[I₁],K¹[I₁+1]))...),
                BezierPathSegment(map(p->LxrPt(p,unitlength),BézPts(t->M([K¹[I₁+1],t]),K²[I₂],K²[I₂+1]))...),
                BezierPathSegment(map(p->LxrPt(p,unitlength),BézPts(t->M([t,K²[I₂+1]]),K¹[I₁+1],K¹[I₁]))...),
                BezierPathSegment(map(p->LxrPt(p,unitlength),BézPts(t->M([K¹[I₁],t]),K²[I₂+1],K²[I₂]))...)])
        mesh1 = Luxor.mesh(BézPth, [
            colorfunc([K¹[I₁], K²[I₂]]),
            colorfunc([K¹[I₁+1], K²[I₂]]),
            colorfunc([K¹[I₁+1], K²[I₂+1]]),
            colorfunc([K¹[I₁], K²[I₂+1]])
            ])
        setmesh(mesh1)
        box(LxrPt([right+left,up+down]/2,unitlength), (right-left)*unitlength,(up-down)*unitlength,:fill)
    end
    finish()
    return nothing
end


text:
using Test, SpecialFunctions
import TruncatedNormal as TN

@test TN.tnmean(-Inf, Inf) == 0
@test TN.tnmean(0, +Inf) ≈ +√(2/π)
@test TN.tnmean(-Inf, 0) ≈ -√(2/π)
@test TN.tnmean(-Inf, -Inf) == -Inf
@test TN.tnmean(+Inf, +Inf) == +Inf

for x = -10:10
    @test TN.tnmean(x, +Inf) ≈ +√(2/π) / erfcx(+x/√2)
    @test TN.tnmean(-Inf, x) ≈ -√(2/π) / erfcx(-x/√2)
end

for x = exp.(-100:100)
    @test TN.tnmean(-x, x) ≈ 0
    @test TN.tnmean(0, x) ≈ -√(2/π) * expm1(-x^2 / 2) / erf(x / √2)
end

@test TN.tnmean(1e-44, 1e-43) ≈ 5.4999999999999999999999999999999999999999e-44

@test TN.tnmean(-Inf, Inf, 0, 1) == TN.tnmean(-Inf, Inf)
for x = -100:100
    @test TN.tnmean(x, +Inf, 0, 1) ≈ TN.tnmean(x, +Inf)
    @test TN.tnmean(-Inf, x, 0, 1) ≈ TN.tnmean(-Inf, x)
end

@test TN.tnmean(-1e4, 1e4, 0, 1) == TN.tnmean(-1e4, 1e4)

for x = -100:100, y = x + 1 : 100
    @test TN.tnmean(x, y, 0, 1) == TN.tnmean(x, y)
end

@test TN.tnmean(100, 115) ≈ 100.00999800099926070518490239457545847490332879043
@test TN.tnmean(-1e6, -999000) ≈ -999000.00000100100100099899498898098
@test TN.tnmean(+1e6, +Inf) ≈ +1.00000000000099999999999800000e6
@test TN.tnmean(-Inf, -1e6) ≈ -1.00000000000099999999999800000e6

@test TN.tnmean(-1e200, 1e200) ≈ 0
@test TN.tnmean(0, +1e200) ≈ +0.797884560802865355879892119869
@test TN.tnmean(-1e200, 0) ≈ -0.797884560802865355879892119869

@test TN.tnmean(50, 70, -2, 3) ≈ 50.171943499898757645751683644632860837133138152489
@test TN.tnmean(-100.0, 0.0, 0.0, 2.0986317998643735) ≈ -1.6744659119217125058885983754999713622460154892645
@test TN.tnmean(0.0, 0.9, 0.0, 0.07132755843183151) ≈ 0.056911157632522598806524588414964004271754161737065
@test TN.tnmean(-100.0, 100.0, 0.0, 17.185261847875548) ≈ 0
@test TN.tnmean(-100.0, 0.5, 0.0, 0.47383322897860064) ≈ -0.1267981330521791493635176736743283314399
@test TN.tnmean(-100.0, 100.0, 0.0, 17.185261847875548) ≈ 0

for x = 1e1 .^ (-400:100:400)
    @test TN.tnmean(-x, x) == TN.tnmean(-x, x) || isnan(TN.tnmean(-x, x)) && isnan(TN.tnmean(-x, x))
end

for a = exp.(-10:10), b = exp.(-10:10)
    a ≤ b || continue
    @test a ≤ TN.tnmean(a, b) ≤ b
end

# https://github.com/JuliaStats/Distributions.jl/issues/827
@test TN.tnmean(0, 1000, 1000000, 1) ≈ 999.99999899899899900100501101901899090472046236710608108591983


text:
"""
    update_fista!(B, B_prev, A, resid, resid_B, grad, X, Y, Z, norms, lambda,
                  reg, iter, stepsize)

Updates coefficient estimates in place for each FISTA iteration based on the Elastic-net
silution, when `X` and `Z` are both standardized.

# Arguments

- B = 2d array of floats consisting of coefficient estimates
- B_prev = 2d array of floats consisting of coefficient estimates saved from
  the previous iteration
- A = 2d array of floats consisting of extrapolated coefficients
- resid = 2d array of floats consisting of the residuals calculated from the
  extrapolated coefficients
- resid_B = 2d array of floats consisting of the residuals calculated from
  the coefficient estimates
- grad = 2d array of floats consisting of the gradient
- X = 2d array of floats consisting of the row covariates, with all
  categorical variables coded in appropriate contrasts
- Y = 2d array of floats consisting of the multivariate response observations
- Z = 2d array of floats consisting of the column covariates, with all
  categorical variables coded in appropriate contrasts
- norms = `nothing`
- lambdaL1 = l1 penalty, a floating scalar
- lambdaL2 = l2 penalty, a floating scalar
- reg = 2d array of bits, indicating whether or not to regularize each of the
  coefficients
- iter = 1d array consisting of a single integer keeping track of how many
  iterations have been computed
- stepsize = 1d array consisting of a float; step size of updates

# Value

None; updates coefficients in place

"""
function update_fista!(B::AbstractArray{Float64,2},
                       B_prev::AbstractArray{Float64,2},
                       A::AbstractArray{Float64,2},
                       resid::AbstractArray{Float64,2},
                       resid_B::AbstractArray{Float64,2},
                       grad::AbstractArray{Float64,2},
                       X::AbstractArray{Float64,2},
                       Y::AbstractArray{Float64,2},
                       Z::AbstractArray{Float64,2},
                       norms::Nothing, lambdaL1::Float64, lambdaL2::Float64,
                       reg::BitArray{2},
                       iter::AbstractArray{Int64,1},
                       stepsize::AbstractArray{Float64,1})

    # Update residuals based on extrapolated coefficients
    calc_resid!(resid, X, Y, Z, A)
    # Update gradient
    calc_grad!(grad, X, Z, resid)

    # Cycle through all the coefficients to perform assignments
    for j = 1:size(B,2), i = 1:size(B,1)
        # Update coefficient estimate in the previous copy
        if B_prev[i,j] != B[i,j]
            B_prev[i,j] = B[i,j]
        end

        B[i,j] = A[i,j] - stepsize[1] * grad[i,j] # RSS updates
        # Apply shrinkage to regularized coefficients
        if reg[i,j]
            B[i,j] = prox(A[i,j], grad[i,j], sign(B[i,j]), lambdaL1, norms,
                          stepsize[1])/(1+2*lambdaL2*stepsize[1])
        end

        # Update extrapolated coefficient
        A[i,j] = B[i,j] + ((iter[1]-1.0)/(iter[1]+2.0))*(B[i,j] - B_prev[i,j])
    end

    # Update residuals based on coefficient estimates
    calc_resid!(resid_B, X, Y, Z, B)
end

"""
    update_fista!(B, B_prev, A, resid, resid_B, grad, X, Y, Z, norms, lambda,
                  reg, iter, stepsize)

Updates coefficient estimates in place for each FISTA iteration when `X` and
`Z` are not standardized.

# Arguments

- B = 2d array of floats consisting of coefficient estimates
- B_prev = 2d array of floats consisting of coefficient estimates saved from
  the previous iteration
- A = 2d array of floats consisting of extrapolated coefficients
- resid = 2d array of floats consisting of the residuals calculated from the
  extrapolated coefficients
- resid_B = 2d array of floats consisting of the residuals calculated from
  the coefficient estimates
- grad = 2d array of floats consisting of the gradient
- X = 2d array of floats consisting of the row covariates, with all
  categorical variables coded in appropriate contrasts
- Y = 2d array of floats consisting of the multivariate response observations
- Z = 2d array of floats consisting of the column covariates, with all
  categorical variables coded in appropriate contrasts
- norms = 2d array of floats consisting of the norms corresponding to each
  coefficient
- lambdaL1 = l1 penalty, a floating scalar
- lambdaL2 = l2 penalty, a floating scalar
- reg = 2d array of bits, indicating whether or not to regularize each of the
  coefficients
- iter = 1d array consisting of a single integer keeping track of how many
  iterations have been computed
- stepsize = 1d array consisting of a float; step size of updates

# Value

None; updates coefficients in place

"""
function update_fista!(B::AbstractArray{Float64,2},
                       B_prev::AbstractArray{Float64,2},
                       A::AbstractArray{Float64,2},
                       resid::AbstractArray{Float64,2},
                       resid_B::AbstractArray{Float64,2},
                       grad::AbstractArray{Float64,2},
                       X::AbstractArray{Float64,2},
                       Y::AbstractArray{Float64,2},
                       Z::AbstractArray{Float64,2},
                       norms::AbstractArray{Float64,2},
                       lambdaL1::Float64, lambdaL2::Float64,
                       reg::BitArray{2},
                       iter::AbstractArray{Int64,1},
                       stepsize::AbstractArray{Float64,1})

    # Update residuals based on extrapolated coefficients
    calc_resid!(resid, X, Y, Z, A)
    # Update gradient
    calc_grad!(grad, X, Z, resid)

    # Cycle through all the coefficients to perform assignments
    for j = 1:size(B,2), i = 1:size(B,1)
        # Update coefficient estimate in the previous copy
        if B_prev[i,j] != B[i,j]
            B_prev[i,j] = B[i,j]
        end

        B[i,j] = A[i,j] - stepsize[1] * grad[i,j]./norms[i,j] # L2 updates
        # Apply shrinkage to regularized coefficients
        if reg[i,j]
            B[i,j] = prox(A[i,j], grad[i,j], sign(B[i,j]), lambdaL1, norms[i,j],
                          stepsize[1])/(1+2*lambdaL2*stepsize[1])
        end

        # Update extrapolated coefficient
        A[i,j] = B[i,j] + ((iter[1]-1.0)/(iter[1]+2.0))*(B[i,j] - B_prev[i,j])
    end

    # Update residuals based on coefficient estimates
    calc_resid!(resid_B, X, Y, Z, B)
end

"""
    fista!(X, Y, Z, lambda, B, regXidx, regZidx, reg, norms;
           isVerbose, stepsize, thresh, maxiter)

Performs the Elastic-net version FISTA with fixed step size.

# Arguments

- X = 2d array of floats consisting of the row covariates, with all
  categorical variables coded in appropriate contrasts
- Y = 2d array of floats consisting of the multivariate response
- Z = 2d array of floats consisting of the column covariates, with all
  categorical variables coded in appropriate contrasts
- lambda = penalty parameter, a floating scalar
- alpha = parameter (ϵ[0, 1]) determining the mix of penalties between L1 and L2
- B = 2d array of floats consisting of starting coefficient estimates
- regXidx = 1d array of indices corresponding to regularized X covariates
- regZidx = 1d array of indices corresponding to regularized Z covariates
- reg = 2d array of bits, indicating whether or not to regularize each of the
  coefficients
- norms = 2d array of floats consisting of the norms corresponding to each
  coefficient or `nothing`

# Keyword arguments

- isVerbose = boolean flag indicating whether or not to print messages.
  Defaults to `true`.
- stepsize = float; step size for updates. Defaults to `0.01`.
- thresh = threshold at which the coefficients are considered to have
  converged, a floating scalar. Defaults to `10^(-7)`.
- maxiter = maximum number of update iterations. Defaults to `10^10`.

# Value

None; updates coefficients in place

# Some notes

Convergence is determined as when the log ratio of the current and previous
criteria is less than the threshold `thresh`.

The default method for choosing the fixed step size for `fista!` is to use the
reciprocal of the product of the maximum eigenvalues of `X*transpose(X)`
and `Z*transpose(Z)`. This is computed by the `mlmnet` function when `fista!`
is passed into the `fun` argument and `setStepSize` is set to `true`. If
`setStepSize` is set to `false`, the value of the `stepsize` argument will be
used as the fixed step size. Note that obtaining the eigenvalues when `X`
and/or `Z` are very large may exceed computational limitations.

"""
function fista!(X::AbstractArray{Float64,2}, Y::AbstractArray{Float64,2},
                Z::AbstractArray{Float64,2},
                lambda::Float64, alpha::Float64,
                B::AbstractArray{Float64,2},
                regXidx::AbstractArray{Int64,1},
                regZidx::AbstractArray{Int64,1}, reg::BitArray{2}, norms;
                isVerbose::Bool=true, stepsize::Float64=0.01,
                thresh::Float64=10.0^(-7), maxiter::Int=10^10)

    # Re-parametrize the Elastic-net tuning parameters
    lambdaL1 = lambda*alpha
    lambdaL2 = lambda*(1-alpha)

    # Pre-allocating arrays to store residuals for extrapolated coefficients
    # and gradient.
    resid = Array{Float64}(undef, size(Y))
    grad = Array{Float64}(undef, size(B))

    # Calculate residuals
    calc_resid!(resid, X, Y, Z, B)
    # Initialize residuals using coefficient estimates
    resid_B = copy(resid)
    # Store coefficients from previous iteration
    B_prev = copy(B)
    # Store extrapolated coefficients
    A = copy(B)

    # Denominators of criterion
    crit_denom = [size(X,1)*size(Z,1), size(X,2)*size(Z,2)]
    # Placeholder to store the old criterion
    oldcrit = 1.0
    # Calculate the current criterion
    crit = criterion(B[regXidx, regZidx], resid_B, lambdaL1, lambdaL2, crit_denom)

    stepsize = [stepsize]
    iter = [0]
    # Iterate until coefficients converge or maximum iterations have been
    # reached.
    while (abs(log(crit/oldcrit)) > thresh) && (iter[1] < maxiter)
        # Store the current criterion
        oldcrit = crit
        # Update the coefficient estimates
        update_fista!(B, B_prev, A, resid, resid_B, grad, X, Y, Z, norms,
                      lambdaL1, lambdaL2, reg, iter, stepsize)
        # Calculate the criterion after updating
        crit = criterion(B[regXidx, regZidx], resid_B, lambdaL1, lambdaL2, crit_denom)

        iter[:] = iter .+ 1 # Increment the number of iterations
        # Warning message if coefficient estimates do not converge.
        if iter==(maxiter)
            println("Estimates did not converge.")
        end
    end

    println_verbose(string("Criterion: ", crit), isVerbose)
    println_verbose(string("Number of iterations: ", iter[1]), isVerbose)
end