text:
module SimplePlotPerf

using SimplePlot

export prplot,precisionrecall,kaplanmeier

function precisionrecall(sortedLabels; samplePoints=nothing)
    p = Float64[]
    r = Float64[]
    totalRight = sum(sortedLabels)
    foundRight = 0
    for i in 1:length(sortedLabels)
        foundRight += sortedLabels[i]
        push!(p, foundRight/i)
        push!(r, foundRight/totalRight)
    end

    if samplePoints == nothing
        return r,p
    else

        # re-sample the curve at the given recall points
        pout = zeros(length(samplePoints))
        ind = 1
        lastp = p[1]
        lastr = r[1]
        for i in 1:length(samplePoints)
            while samplePoints[i] > r[ind] && ind < length(r)
                lastp = p[ind]
                lastr = r[ind]
                ind += 1
            end
            if samplePoints[i] <= r[ind]
                # pick the nearest recall point
                pout[i] = r[ind]-samplePoints[i] < lastr-samplePoints[i] ? p[ind] : lastp
            else
                pout[i] = lastp # we are past any recorded recall points
            end
        end
        return samplePoints,pout
    end
end
function precisionrecall(pred, labels; samplePoints=nothing)
    @assert length(pred) == length(labels)
    precisionrecall(labels[sortperm(pred, rev=true)], samplePoints=samplePoints)
end

function prplot(methods::AbstractArray, truth; samplePoints=linspace(0,1,200), auc=true,
        colors=SimplePlot.colors,
        randomMeanLayer=line(),
        randomSampleLayer=line(), kwargs...)
    rmean = zeros(length(samplePoints))
    pmean = zeros(length(samplePoints))
    layers = Any[]
    numRandom = 10
    for i in 1:numRandom
        r,p = precisionrecall(rand(length(truth)), truth, samplePoints=samplePoints)
        rmean .+= r
        pmean .+= p

        layerParams = merge(Dict(
            :color => "#cccccc"
        ), randomSampleLayer.params, randomSampleLayer.unusedParams)
        push!(layers, line(r, p; layerParams...))
    end
    rmean ./= numRandom
    pmean ./= numRandom
    aucRand = @sprintf("%0.03f", area_under_curve(rmean, pmean))

    methodLayers = Any[]
    for (ind,method) in enumerate(methods)
        r,p = nothing,nothing
        if typeof(method[2][1]) <: Array
            r = zeros(length(samplePoints))
            p = zeros(length(samplePoints))

            for member in method[2]
                rtmp,ptmp = precisionrecall(member, truth, samplePoints=samplePoints)
                r .+= rtmp
                p .+= ptmp
                layerParams = merge(Dict(
                    :color => colors[ind],
                    :alpha => 0.5,
                ))
                push!(layers, line(rtmp, ptmp; layerParams...))
            end
            r ./= length(method[2])
            p ./= length(method[2])
        else
            r,p = precisionrecall(method[2], truth, samplePoints=samplePoints)
        end

        aucStr = @sprintf("%0.03f", area_under_curve(r, p))
        layerParams = Dict{Any,Any}(
            :linewidth => 2,
            :color => colors[ind],
            :label => method[1]*(auc ? " (AUC $aucStr)" : "")
        )
        if length(method) > 2
            @assert typeof(method[3]) <: SimplePlot.LineLayer "Third slot for a method must be a template line layer"
            merge!(layerParams, method[3].params, method[3].unusedParams)
        end
        push!(methodLayers, line(r, p; layerParams...))
    end

    aucStr = @sprintf("%0.03f", area_under_curve(rmean, pmean))
    layerParams = merge(Dict(
        :color => "#666666",
        :linewidth => 2,
        :label => "Random"*(auc ? " (AUC $aucStr)" : "")
    ), randomMeanLayer.params, randomMeanLayer.unusedParams)
    axis(
        methodLayers...,
        line(rmean, pmean; layerParams...),
        layers...;
        xlabel="recall",
        ylabel="precision",
        kwargs...
    )
end

function kaplanmeier(times, events)
    xs = Float64[]
    ys = Float64[]

    probSurv = 1
    N = length(times)
    numEventsNow = 0
    lastTime = -1.0

    count = 0
    for i in sortperm(times)
        count += 1

        # if we have a new unique time point
        if times[i] != lastTime && lastTime >= 0.0
            numLeft = N-count+1
            probSurv *= (numLeft - numEventsNow)/numLeft
            numEventsNow = 0
            push!(xs, times[i])
            push!(ys, probSurv)
        end

        lastTime = times[i]
        numEventsNow += 1
    end
    xs,ys
end

# must be sorted by increasing x
function area_under_curve(x, y)
    area = 0.0
    for i in 2:length(x)
        area += (y[i-1]+y[i])/2 * (x[i]-x[i-1])
    end
    area
end

end # module


text:
struct ConceptNet{L<:Language, K<:AbstractString, E<:Real}
    embeddings::Dict{L, Dict{K,Vector{E}}}
    width::Int
    fuzzy_words::Dict{L, Vector{K}}
end

ConceptNet(embeddings::Dict{K,Vector{E}}, width::Int) where
        {K<:AbstractString, E<:Real} =
    ConceptNet{Languages.English(), K, E}(embeddings, width, Dict(Languages.English()=>K[]))

# Aliases
const ConceptNetMulti{L,E<:Real} = ConceptNet{L, String, E}
const ConceptNetMultiCompressed{L} = ConceptNet{L, String, Int8}
const ConceptNetEnglish{E} = ConceptNet{Languages.English, String, E}

# Show methods
show(io::IO, conceptnet::ConceptNetMultiCompressed{L}) where L = begin
    nlanguages = length(conceptnet.embeddings)
    print(io, "ConceptNet{$L} (compressed): $nlanguages language(s)",
          ", $(length(conceptnet)) embeddings")
end

show(io::IO, conceptnet::ConceptNetMulti{L,E}) where {L,E} = begin
    nlanguages = length(conceptnet.embeddings)
    print(io, "ConceptNet{$L,$E}: $nlanguages language(s)",
          ", $(length(conceptnet)) embeddings")
end

show(io::IO, conceptnet::ConceptNetEnglish{E}) where E =
    print(io, "ConceptNet{English,$E}: $(length(conceptnet)) embeddings")

# Overloaded `get` method for a ConceptNet language dictionary
# Example: the embedding corresponding to "###_something" is returned for any search query
#          of two words where the first word in made out out 3 letters followed by
#          the word 'something'
function get(embeddings::Dict{K,Vector{E}},
             keyword::K,
             default::Vector{E},
             fuzzy_words::Vector{K};
             wildcard_matching::Bool=true) where
        {K<:AbstractString, E<:Real}
    if haskey(embeddings, keyword)
        # The keyword exists in the dictionary
        return embeddings[keyword]
    else
        if wildcard_matching
            # The keyword is not found; try fuzzy matching
            ω = 0.4 # weight assinged to matching a #, 1-w weight assigned to a matching letter
            L = length(keyword)
            matches = (word for word in fuzzy_words
                       if length(word) == L &&
                          occursin(Regex(replace(word,"#"=>".")), keyword))
            if isempty(matches)
                return default
            else
                best_match = ""
                max_score = 0
                for match in matches
                    l = length(replace(match,"#"=>"")) # number of letters matched
                    score = ω*(L-l)/L + (1-ω)*l/L
                    if score > max_score
                        best_match = match
                        max_score = score
                    end
                end
                return embeddings[best_match]
            end
        else
            # The keyword is not found; no fuzzy matching
            return default
        end
    end
end

function get(embeddings::Dict{K,Vector{E}},
             keywords::AbstractVector{K},
             default::Vector{E},
             fuzzy_words::Vector{K};
             wildcard_matching::Bool=true,
             n::Int=0) where
        {K<:AbstractString, E<:Real}
    p = length(keywords)
    keywords_embedded = Matrix{E}(undef, n, p)
    for i in 1:p
        keywords_embedded[:,i] = get(embeddings,
                                     keywords[i],
                                     default,
                                     fuzzy_words,
                                     wildcard_matching=wildcard_matching)
    end
    return keywords_embedded
end

# Indexing
# Generic indexing, multiple words
# Example: julia> conceptnet[Languages.English(), ["another", "word"])
getindex(conceptnet::ConceptNet{L,K,E}, language::L, words::S) where
        {L<:Language, K, E<:Real, S<:AbstractVector{<:AbstractString}} =
    get(conceptnet.embeddings[language],
        words,
        zeros(E, conceptnet.width),
        conceptnet.fuzzy_words[language],
        wildcard_matching=true,
        n=conceptnet.width)

# Generic indexing, multiple words
# Example: julia> conceptnet[:en, ["another", "word"]]
getindex(conceptnet::ConceptNet{L,K,E}, language::Symbol, words::S) where
        {L<:Language, K, E<:Real, S<:AbstractVector{<:AbstractString}} =
    conceptnet[LANGUAGES[language], words]

# Generic indexing, single word
# Example: julia> conceptnet[Languages.English(), "word"]
getindex(conceptnet::ConceptNet{L,K,E}, language::L, word::S) where
        {L<:Language, K, E<:Real, S<:AbstractString} =
    get(conceptnet.embeddings[language],
        word,
        zeros(E, conceptnet.width),
        conceptnet.fuzzy_words[language],
        wildcard_matching=true)

# Generic indexing, single word
# Example: julia> conceptnet[:en, "word"]
getindex(conceptnet::ConceptNet{L,K,E}, language::Symbol, word::S) where
        {L<:Language, K, E<:Real, S<:AbstractString} =
    conceptnet[LANGUAGES[language], word]

# Single-language indexing: conceptnet[["another", "word"]], if language==Languages.English()
getindex(conceptnet::ConceptNet{L,K,E}, words::S) where
        {L<:Languages.Language, K, E<:Real, S<:AbstractVector{<:AbstractString}} =
    conceptnet[L(), words]

# Single-language indexing: conceptnet["word"], if language==Languages.English()
getindex(conceptnet::ConceptNet{L,K,E}, word::S) where
        {L<:Languages.Language, K, E<:Real, S<:AbstractString} =
    conceptnet[L(), word]

# Index by language (returns a Dict{word=>embedding})
getindex(conceptnet::ConceptNet, language::L) where {L<:Languages.Language} =
    conceptnet.embeddings[language]

# Index by language (returns a Dict{word=>embedding})
getindex(conceptnet::ConceptNet, language::Symbol) =
    conceptnet.embeddings[LANGUAGES[language]]

# length methods
length(conceptnet::ConceptNet) = begin
    if !isempty(conceptnet.embeddings)
        return mapreduce(x->length(x[2]), +, conceptnet.embeddings)
    else
        return 0
    end
end

# size methods
size(conceptnet::ConceptNet) = (conceptnet.width, length(conceptnet))

size(conceptnet::ConceptNet, inds...) = (conceptnet.width, length(conceptnet))[inds...]

# in
function in(key::S, conceptnet::ConceptNet) where S<:AbstractString
    found = false
    for lang in keys(conceptnet.embeddings)
        if haskey(conceptnet.embeddings[lang], key)
            found = true
            break
        end
    end
    return found
end

function in(lang::L, conceptnet::ConceptNet) where L<:Languages.Language
    return haskey(conceptnet.embeddings, lang)
end

# Keys
keys(conceptnet::ConceptNet) =
    Iterators.flatten(keys(conceptnet.embeddings[lang]) for lang in keys(conceptnet.embeddings))

ext: jl


text:
function XC_c_vwn_spin( Rhoe, zeta )

    third = 1.0/3.0
    pi34 = 0.6203504908994  # pi34=(3/4pi)^(1/3)
    rs = pi34/Rhoe^third

    A      = ( 0.0310907, 0.01554535, -0.01688686394039 )
    x0     = ( -0.10498, -0.32500, -0.0047584 )
    b      = ( 3.72744, 7.06042, 1.13107 )
    c      = ( 12.9352, 18.0578, 13.0045 )
    Q      = ( 6.15199081975908, 4.73092690956011, 7.12310891781812 )
    tbQ    = ( 1.21178334272806, 2.98479352354082, 0.31757762321188 )
    fx0    = ( 12.5549141492, 15.8687885, 12.99914055888256 )
    bx0fx0 = ( -0.03116760867894, -0.14460061018521, -0.00041403379428 )

    # N.B.: A is expressed in Hartree
    # Q = sqrt(4*c - b^2)
    # tbQ = 2*b/Q
    # fx0 = X(x_0) = x_0^2 + b*x_0 + c
    # bx0fx0 = b*x_0/X(x_0)

    cfz = 2.0^(4.0/3.0) - 2.0
    cfz1 = 1.0 / cfz
    cfz2 = 4.0/3.0 * cfz1
    iddfz0 = 9.0 / 8.0 *cfz
    sqrtrs = sqrt(rs)
    zeta3 = zeta^3
    zeta4 = zeta3*zeta
    trup = 1.0 + zeta
    trdw = 1.0 - zeta
    trup13 = trup^(1.0/3.0)
    trdw13 = trdw^(1.0/3.0)
    fz = cfz1 * (trup13*trup + trdw13*trdw - 2.0)         # f(zeta)
    dfz = cfz2 * (trup13 - trdw13)     # d f / d zeta

    ecP, vcP = padefit( sqrtrs, 1, x0, Q, b, c, A, tbQ, bx0fx0 )    # ecF = e_c Paramagnetic
    ecF, vcF = padefit( sqrtrs, 2, x0, Q, b, c, A, tbQ, bx0fx0 )    # ecP = e_c Ferromagnetic
    ac, dac  = padefit( sqrtrs, 3, x0, Q, b, c, A, tbQ, bx0fx0 )    # ac = "spin stiffness"

    ac = ac * iddfz0
    dac = dac * iddfz0
    De = ecF - ecP - ac # e_c[F] - e_c[P] - alpha_c/(ddf/ddz(z=0))
    fzz4 = fz * zeta4
    ec = ecP + ac * fz  + De * fzz4

    dec1 = vcP + dac*fz + (vcF - vcP - dac) * fzz4       # e_c - (r_s/3)*(de_c/dr_s)
    dec2 = ac*dfz + De*(4.0*zeta3*fz + zeta4*dfz)      # de_c/dzeta

    # v_c[s] = e_c - (r_s/3)*(de_c/dr_s) + [sign(s)-zeta]*(de_c/dzeta)
    vcup = dec1 + (1.0 - zeta)*dec2
    vcdw = dec1 - (1.0 + zeta)*dec2

    return ec, vcup, vcdw
end

function padefit(x, i, x0, Q, b, c, A, tbQ, bx0fx0)
    # implements formula [4.4] in:
    # S.H. Vosko, L. Wilk, and M. Nusair, Can. J. Phys. 58, 1200 (1980)

    # Pade fit calculated in x and its derivative w.r.t. rho
    # rs = inv((rho*)^(1/3)) = x^2
    # fit  [eq. 4.4]
    # dfit/drho = fit - (rs/3)*dfit/drs = ec - (x/6)*dfit/dx

    sqx = x * x                          # x^2 = r_s
    xx0 = x - x0[i]                      # x - x_0
    Qtxb = Q[i] / (2.0*x + b[i])      # Q / (2x+b)
    atg = atan(Qtxb)                     # tan^-1(Q/(2x+b))
    fx = sqx + b[i]*x + c[i]             # X(x) = x^2 + b*x + c

    fit = A[i] * (  log(sqx/fx) + tbQ[i]*atg -
            bx0fx0[i] * ( log(xx0*xx0/fx) + (tbQ[i] + 4.0*x0[i]/Q[i]) * atg )  )

    txb = 2.0*x + b[i]
    txbfx = txb / fx
    itxbQ = 1.0 / (txb*txb + Q[i]*Q[i])

    dfit = fit - A[i] / 3.0 + A[i]*x/6.0 * (  txbfx + 4.0*b[i]*itxbQ +
              bx0fx0[i] * ( 2.0/xx0 - txbfx - 4.0*( b[i] + 2.0*x0[i])*itxbQ )  )

    return fit, dfit
end


text:
mutable struct CodeGen
    context::LLVM.Context
    builder::LLVM.Builder
    jit::LLVM.ExecutionEngine
    current_scope::CurrentScope
    function_protos::Dict{String, PrototypeAST}
    mod::LLVM.Module
    pass_manager::LLVM.FunctionPassManager
    function CodeGen(context::LLVM.Context, builder::LLVM.Builder, jit::LLVM.ExecutionEngine,
                     current_scope::CurrentScope, function_protos::Dict{String, PrototypeAST})
        new(context, builder, jit, current_scope, function_protos)
    end
end

function CodeGen()
    ctx = LLVM.Context()
    builder = LLVM.Builder(ctx)
    jit = LLVM.JIT(LLVM.Module(""))
    return CodeGen(
        ctx,
        builder,
        jit,
        CurrentScope(),
        Dict{String, PrototypeAST}(),
    )
end

current_scope(cg::CodeGen) = cg.current_scope

function new_scope(f, cg::CodeGen)
    open_scope!(current_scope(cg))
    f()
    pop!(current_scope(cg))
end
Base.show(io::IO, cg::CodeGen) = print(io, "CodeGen")

function initialize_module_and_pass_manager!(cg::CodeGen)
    cg.mod = LLVM.Module("KaleidoscopeModule")
    cg.pass_manager = LLVM.FunctionPassManager(cg.mod)
    # data layout
    LLVM.instruction_combining!(cg.pass_manager)
    LLVM.reassociate!(cg.pass_manager)
    LLVM.gvn!(cg.pass_manager)
    LLVM.cfgsimplification!(cg.pass_manager)
    LLVM.promote_memory_to_register!(cg.pass_manager)
    LLVM.initialize!(cg.pass_manager)
    return cg
end

function get_function(cg::CodeGen, name::String)
    if haskey(LLVM.functions(cg.mod), name)
        return LLVM.functions(cg.mod)[name]
    end

    if haskey(cg.function_protos, name)
        return codegen(cg, cg.function_protos[name])
    end

    error("encountered undeclared function $name")
end

function create_entry_block_allocation(fn::LLVM.Function, varname::String)
    local alloc
    LLVM.Builder() do builder
        entry_block = LLVM.entry(fn)
        position_before_terminator!(builder, entry_block)
        alloc = LLVM.alloca!(builder, LLVM.DoubleType(), varname)
    end
    return alloc
end

function codegen(cg::CodeGen, expr::NumberExprAST)::LLVM.Value
    return LLVM.ConstantFP(LLVM.DoubleType(), expr.val)
end

function codegen(cg::CodeGen, expr::VariableExprAST)::LLVM.Value
    # TODO: Error handling if argument is not found
    V = get(current_scope(cg), expr.name, nothing)
    V == nothing && error("did not find variable $(expr.name)")
    return LLVM.load!(cg.builder, V, expr.name)
end

function codegen(cg::CodeGen, expr::BinaryExprAST)::LLVM.Value
    if expr.op == Kinds.EQUAL
        var = expr.lhs
        if !(var isa VariableExprAST)
            error("destination of '=' must be a variable")
        end
        R = codegen(cg, expr.rhs)
        V = get(current_scope(cg), var.name, nothing)
        V == nothing && error("unknown variable name $(var.name)")
        LLVM.store!(cg.builder, R, V)
        return R
    end
    L = codegen(cg, expr.lhs)
    R = codegen(cg, expr.rhs)

    if expr.op == Kinds.PLUS
        return LLVM.fadd!(cg.builder, L, R, "addtmp")
    elseif expr.op == Kinds.MINUS
        return LLVM.fsub!(cg.builder, L, R, "subtmp")
    elseif expr.op == Kinds.STAR
        return LLVM.fmul!(cg.builder, L, R, "multmp")
    elseif expr.op == Kinds.SLASH
        return LLVM.fdiv!(cg.builder, L, R, "divtmp")
    elseif expr.op == Kinds.LESS
        L = LLVM.fcmp!(cg.builder, LLVM.API.LLVMRealOLT, L, R, "cmptmp")
        return LLVM.uitofp!(cg.builder, L, LLVM.DoubleType(), "booltmp")
    elseif expr.op == Kinds.GREATER
        L = LLVM.fcmp!(cg.builder, LLVM.API.LLVMRealOGT, L, R, "cmptmp")
        return LLVM.uitofp!(cg.builder, L, LLVM.DoubleType(), "booltmp")
    else
        error("Unhandled binary operator $(expr.op)")
    end
end

function codegen(cg::CodeGen, expr::CallExprAST)::LLVM.Value
    func = get_function(cg, expr.callee)

    if length(LLVM.parameters(func)) != length(expr.args)
        error("number of parameters mismatch")
    end

    args = LLVM.Value[]
    for v in expr.args
        push!(args, codegen(cg, v))
    end

    return LLVM.call!(cg.builder, func, args, "calltmp")
end

function codegen(cg::CodeGen, expr::PrototypeAST)::LLVM.Value
    if haskey(LLVM.functions(cg.mod), expr.name)
        func = LLVM.functions(cg.mod)[expr.name]

        if length(LLVM.parameters(func)) != length(expr.args)
            error("number of parameters mismatch")
        end

        if length(LLVM.blocks(func)) != 0
            error("Existing function exists with a body")
        end
    else
        args = [LLVM.DoubleType() for i in 1:length(expr.args)]
        func_type = LLVM.FunctionType(LLVM.DoubleType(), args)
        func = LLVM.Function(cg.mod, expr.name, func_type)
        LLVM.linkage!(func, LLVM.API.LLVMExternalLinkage)

        for (i, param) in enumerate(LLVM.parameters(func))
            LLVM.name!(param, expr.args[i])
        end
    end
    return func
end

function codegen(cg::CodeGen, expr::FunctionAST)::LLVM.Value
    cg.function_protos[expr.proto.name] = expr.proto
    the_function = get_function(cg, expr.proto.name)

    entry = LLVM.BasicBlock(the_function, "entry")
    LLVM.position!(cg.builder, entry)

    for (i, param) in enumerate(LLVM.parameters(the_function))
        argname = expr.proto.args[i]
        alloc = create_entry_block_allocation(the_function, argname)
        LLVM.store!(cg.builder, param, alloc)
        current_scope(cg)[argname] = alloc
    end

    new_scope(cg) do
        body = codegen(cg, expr.body)
        # TODO, delete function on error
        LLVM.ret!(cg.builder, body)
        status = convert(Core.Bool, LLVM.API.LLVMVerifyFunction(LLVM.ref(the_function), LLVM.API.LLVMPrintMessageAction))
        if status
            throw(LLVM.LLVMException("broken function: $the_function"))
        end
    end
    LLVM.run!(cg.pass_manager, the_function)
    return the_function
end

function codegen(cg::CodeGen, expr::IfExprAST)
    func = LLVM.parent(LLVM.position(cg.builder))
    then = LLVM.BasicBlock(func, "then")
    elsee = LLVM.BasicBlock(func, "else")
    merge = LLVM.BasicBlock(func, "ifcont")

    local phi
    new_scope(cg) do
        # if
        cond = codegen(cg, expr.cond)
        zero = LLVM.ConstantFP(LLVM.DoubleType(), 0.0)
        condv = LLVM.fcmp!(cg.builder, LLVM.API.LLVMRealONE, cond, zero, "ifcond")
        LLVM.br!(cg.builder, condv, then, elsee)

        # then
        LLVM.position!(cg.builder, then)
        thencg = codegen(cg, expr.then)
        LLVM.br!(cg.builder, merge)
        then_block = position(cg.builder)

        # else
        LLVM.position!(cg.builder, elsee)
        elsecg = codegen(cg, expr.elsee)
        LLVM.br!(cg.builder, merge)
        else_block = position(cg.builder)

        # merge
        LLVM.position!(cg.builder, merge)
        phi = LLVM.phi!(cg.builder, LLVM.DoubleType(), "iftmp")
        append!(LLVM.incoming(phi), [(thencg, then_block), (elsecg, else_block)])
    end

    return phi
end

function codegen(cg::CodeGen, expr::ForExprAST)

    new_scope(cg) do
        startblock = position(cg.builder)
        func = LLVM.parent(startblock)
        alloc = create_entry_block_allocation(func, expr.varname)
        start = codegen(cg, expr.start)
        LLVM.store!(cg.builder, start, alloc)
        loopblock = LLVM.BasicBlock(func, "loop")
        LLVM.br!(cg.builder, loopblock)
        LLVM.position!(cg.builder, loopblock)
        current_scope(cg)[expr.varname] = alloc
        codegen(cg, expr.body)
        step = codegen(cg, expr.step)
        endd = codegen(cg, expr.endd)

        curvar = LLVM.load!(cg.builder, alloc, expr.varname)
        nextvar = LLVM.fadd!(cg.builder, curvar, step, "nextvar")
        LLVM.store!(cg.builder, nextvar, alloc)

        endd = LLVM.fcmp!(cg.builder, LLVM.API.LLVMRealONE, endd,
            LLVM.ConstantFP(LLVM.DoubleType(), 1.0))

        loopendblock = position(cg.builder)
        afterblock = LLVM.BasicBlock(func, "afterloop")

        LLVM.br!(cg.builder, endd, loopblock, afterblock)
        LLVM.position!(cg.builder, afterblock)
    end

    return LLVM.ConstantFP(LLVM.DoubleType(), 0.0)
end

function codegen(cg::CodeGen, expr::VarExprAST)
    func = LLVM.parent(LLVM.position(cg.builder))
    local initval
    for (varname, init) in expr.varnames
        initval = codegen(cg, init)
        alloca = create_entry_block_allocation(func, varname)
        LLVM.store!(cg.builder, initval, alloca)
        current_scope(cg)[varname] = alloca
    end
    return initval
end

function codegen(cg::CodeGen, expr::BlockExprAST)
    local v
    # what about empty blocks?
    new_scope(cg) do
        for expr in expr.exprs
            v = codegen(cg, expr)
        end
    end
    return v
end


text:
import Pfam

msa = Pfam.MSA("PF00011")
df = Pfam.load(msa)

@test size(df, 1) > 10000
@test size(df, 2) == 4
@test "H0TRS9_9BRAD" ∈ df.uniprot_id


text:
#=
# Rfam.jl usage example
=#

import Rfam

# Use `Rfam.fetch` to get the set of sequences in a family, as a
# Vector{FASTA.Record}:

fasta = Rfam.fetch("RF00162")


text:
"""    approximate!(dvars, svars, Z, Δτ, n)

Assigns to the i-th columns of `F`, for `i = 1:n` the short-time limit of the
intermediate scattering function, and calculates the memory function for the
same time-points.
"""
function approximate!(dvars, svars, Z, Δτ, n)
    @unpack F, Fˢ, ζ = dvars
    @unpack S, Sˢ, B, Bˢ, w, υ = svars

    @inbounds for i = 1:n
        τ = i * Δτ

        for j in eachindex(S)
            F[i, j]  = short_times_limit!(S[j], B[j], τ)
            Fˢ[i, j] = short_times_limit!(Sˢ[j], Bˢ[j], τ)
            Z[j] = product(w[j], F[i, j], Fˢ[i, j])
        end

        ζ[i] = υ * sum(Z)
    end

    return dvars
end

short_times_limit!(Sⱼ, β, τ) = exp(-τ * reduce_dof(β)) * Sⱼ


text:
const XC_LDA_X = 1   # Exchange
const XC_LDA_C_VWN = 7  # Vosko, Wilk, & Nusair (5)
const XC_GGA_X_PBE = 101
const XC_GGA_C_PBE = 130

function test_main(id_x, id_c)

  ρ = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6]
  σ = [0.2, 0.3, 0.4, 0.5, 0.6, 0.7]
  N = size(ρ)[1]

  ϵ_x = zeros(N)
  ϵ_c = zeros(N)

  ccall( (:ffrlibxc, "ffrlibxc"), Void,
         (Int, Int, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}),
         id_x, N, ρ, σ, ϵ_x )

  ccall( (:ffrlibxc, "ffrlibxc"), Void,
         (Int, Int, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}),
          id_c, N, ρ, σ, ϵ_c )

  @printf("\n")
  for i = 1:N
    @printf("%d %f %f %f\n", i, ρ[i], σ[i], ϵ_x[i] + ϵ_c[i])
  end

end

include("LDA_VWN.jl")
function test_my_VWN()
  ρ = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6]
  ϵ = excVWN(ρ)
  N = size(ρ)[1]
  for i = 1:N
    @printf("%d %f %f\n", i, ρ[i], ϵ[i])
  end
end

test_main(XC_LDA_X, XC_LDA_C_VWN)

test_my_VWN()

#test_main(XC_GGA_X_PBE, XC_GGA_C_PBE)


text:
module AboutWindow

using Gtk, Gtk.ShortNames, Gtk.GConstants, Pkg
using ..Ahorn

window = nothing

function spawnWindowIfAbsent!()
    if window === nothing
        icon = Pixbuf(filename = Ahorn.logoFile, width = 256, height = 256, preserve_aspect_ratio = true)
        version = ""
        try
            local ctx = Pkg.Types.Context()
            local pkg = PackageSpec(name="Ahorn", uuid=ctx.env.project.deps["Ahorn"])
            Pkg.Operations.resolve_versions!(ctx, [pkg])
            version *= "Hash " * string(pkg.repo.tree_sha)[1:7]
        catch e
            version *= "Don't expect a version number"
        end
        global window = AboutDialog(title = "About $(Ahorn.baseTitle)",
        icon = Ahorn.windowIcon, logo = icon,
        program_name = Ahorn.baseTitle,
        comments="Level Editor for the Game Celeste\n\nPowered by Julia and Maple",
        version = version,
        website = "https://github.com/CelestialCartographers/Ahorn",
        website_label = "GitHub  Repository")

        # Hide window instead of destroying it
        signal_connect(hideAboutWindow, window, "delete_event")
        signal_connect(hideAboutWindow, window, "response")

        GAccessor.transient_for(window, Ahorn.window)
        GAccessor.authors(window, ["Cruor", "Vexatos"])
        GAccessor.artists(window, ["Vexatos"])
    end
end

function hideAboutWindow(widget=nothing, event=nothing)
    if event == GtkResponseType.GTK_RESPONSE_DELETE_EVENT || isa(event, Gtk.GdkEvent)
        spawnWindowIfAbsent!()
        visible(window, false)
        return true
    end
end

function showAboutWindow(widget)
    spawnWindowIfAbsent!()
    visible(window, true)
    return true
end

end


text:
maxdims["trg"] = 10:10:50
descriptions["trg"] = "TRG, 2D classical Ising model\nN → ∞, 20 iterations\nβ = 1.001 βc"

function runbenchmark(::Val{:trg};
                      maxdim::Int, nsweeps::Int = 20,
                      outputlevel::Int = 0, cutoff::Float64 = 0.0,
                      β::Float64 = 1.001 * βc, splitblocks::Bool = false)
  if splitblocks
    println("Benchmark trg doesn't support splitblocks $splitblocks.")
    return nothing
  end
  # Make Ising model MPO
  s = Index(2)
  sₕ = addtags(s, "horiz")
  sᵥ = addtags(s, "vert")
  T = ising_mpo(sₕ, sᵥ, β)
  κ, T = trg(T; χmax = maxdim, cutoff = cutoff, nsteps = nsweeps)
  if outputlevel > 0
    @show nsweeps
    @show β
    @show cutoff
    @show κ, exp(-β * ising_free_energy(β))
  end
  return ITensors.maxdim(T)
end


text:
# Automatically generated using Clang.jl

const GSW_INVALID_VALUE = 9.0e15
const GSW_ERROR_LIMIT = 1.0e10
const INTERP_METHOD_LINEAR = 1
const INTERP_METHOD_PCHIP = 2


text:
include("../evo/data.jl")

function test_classification(cont_constructor::Function;
                             cfg::Dict=Config())
    c = cont_constructor(cfg)
    m = Model(cfg, c)
    X, Y = get_iris(0)
    nin = size(X, 1)
    nout = length(unique(Y))
    layered_init!(m, nin, nout; nhidden=nin, nreward=1)
    fitness = classify(m, X, Y)
    println(fitness)
    @test fitness[1] >= -1.0
end

# @testset "Random controller" begin
#     test_classification(rand_controller)
# end

@testset "Static controller" begin
    test_classification(static_controller; cfg = Config(Config()))
end

# @testset "SNN controller" begin
#     test_classification(snn_controller;
#                         cfg = Config(Config(), "cfg/snn.yaml"))
# end

# @testset "CGP controller" begin
#     include("../evo/cgp.jl")
#     test_classification(cgp_controller)
# end


text:
import Base.print, Base.show

export print, print_pointsummary

function print(io::IO, pb::SDCOContext)
    println(io, " * Nb of scalar variables:     $(pb.nscalvar)")
    println(io, " * Nb of matrix variables:     $(length(pb.nsdpvar))")
    println(io, "         individual sizes:     $(pb.nsdpvar)")
    println(io, " * Nb of constraints:          $(pb.m)")
    println(io)
    println(io, " * Objective:")
    print(io, " ", pb.c)
    println(io)
    println(io, " * Constraints:")
    for (i, ai) in enumerate(pb.A)
        println(io, " **** ctr $i: rhs = $(pb.b[i])")
        print(io, ai)
    end
    println(io)
    for (opt, val) in pb.options
        println(io, opt, "   => ", val)
    end
end

function print(io::IO, pt::PointE)
    for (i, mat) in enumerate(pt.mats)
        print(io, " - mat $i: ")
        display(mat)
    end
    print(io, " - vec: ")
    display(pt.vec)
end

function show(io::IO, pt::PointE)
    for (i, mat) in enumerate(pt.mats)
        print(io, " - mat $i - ")
        display(mat)
    end
    print(io, " - vec - ")
    display(pt.vec)
end

function show(io::IO, z::PointPrimalDual)
    println(io, "   * Primal space - x:")
    display(z.x)
    println(io, "   * Dual space - y:")
    display(z.y)
    println(io, "   * Dual space - s:")
    display(z.s)
end

function print_pointsummary(pb::SDCOContext, z::PointPrimalDual)
    println("Primal / dual obj      :    ", get_primobj(pb, z), " / ", get_dualobj(pb, z))
    println("Primal feasability err :    ", get_primfeaserr(pb, z))
    println("Dual feasability err   :    ", get_dualfeaserr(pb, z))
    println("Min_K(x)               :    ", min(pb, z.x))
    println("Min_K(s)               :    ", min(pb, z.s))
    println("x . s                  :    ", hadamard(z.x, z.s))
end


text:
using Test, Random, DataFrames
using Nanocsv

function rwtest(df; sep=',', na="")
    fname = "nanocsv_test.txt"
    try
        write_csv("test.csv", df, sep=sep, na=na)
        df2 = read_csv("test.csv", sep=sep, na=na)
        @test isequal(df, df2)
    finally
        isfile(fname) && rm(fname)
    end
end

function rwthrows(df; err=ArgumentError, sep=',', na="")
    fname = "nanocsv_test.txt"
    try
        write_csv("test.csv", df, sep=sep, na=na)
        @test_throws err read_csv("test.csv", sep=sep, na=na)
    finally
        isfile(fname) && rm(fname)
    end
end

@testset "basic read and write" begin
    Random.seed!(1)
    rwtest(DataFrame(rand(Int, 5, 5)))
    rwtest(DataFrame(rand(Int, 5, 5)), sep=';', na="NA")
end

@testset "some more complex test" begin
    df = DataFrame(i = [1,2,missing], f=[missing,1.0, 2.0], m=missing,
                   s1=["s1",",",";"], s2=["s2",missing,missing])
    rwtest(df)
    rwtest(df, sep=';', na="NA")
end

@testset "minimal DataFrames" begin
    rwthrows(DataFrame())
    rwthrows(DataFrame(), sep=';', na="NA")
    rwtest(DataFrame(:a=>[], :b=>[], Symbol("")=>[]))
    rwtest(DataFrame(:a=>[], :b=>[], Symbol("")=>[]), sep=';', na="NA")
    rwtest(DataFrame(Symbol("")=>missing))
    rwtest(DataFrame(Symbol("")=>missing), sep=';', na="NA")
    rwtest(DataFrame(Symbol("")=>["",missing,""]))
    rwtest(DataFrame(Symbol("")=>["",missing,""]), sep=';', na="NA")
end


text:
module VarBCFiles

using Dates, Printf, DataStructures, DataFrames
import Base: read, write, show, size, getindex, ==, merge!, length, filter!
import DataFrames: DataFrame
export VarBC, VarBCRecord, increment, label

include("VarBCRecord.jl")

struct VarBC  # <: AbstractDict{String, VarBCRecord}
    datetime::DateTime
    version::String
    header1::String
    header2::String
    records::OrderedDict{String,VarBCRecord}
end

length(a::VarBC) = length(a.records)
show(io::IO, ::MIME"text/plain", a::VarBC) = println(io, "VarBC with $(length(a)) records for $(a.datetime)")

"""
    filter!(f,a)

Update VarBCRecords in `a`, removing elements for which `f` is `false`. The function `f` is passed `key=>value` pairs

# Example
```jldoctest

julia> filter!(p-> p.second.label == "sfcobs", a)

```
"""
filter!(f,a::VarBC) = filter!(f,a.records)

"""
    merge!(a,b)

Merges the VarBC Records from `b` into `a`. If a key is present in both `a` and `b` the value from `b` is used.
Header information in `a` is not taken from `b`
"""
merge!(a::VarBC,b::VarBC) = merge!(a.records,b.records)

"""
    a == b

Returns `true` if `a` and `b` are equal
"""
function ==(a::VarBC,b::VarBC)
    a.datetime == b.datetime &&
    a.header1 == b.header1 &&
    b.header2 == b.header2 &&
    a.version == a.version &&
    a.records == a.records
end

"""
    read(filename,VarBC)

returns a object of type `VarBC`
"""
function read(fname::String,::Type{VarBC})
    io = open(fname)
    version = readline(io); @assert version == "VARBC_cycle.version006"
    header1 = readline(io)
    header2 = readline(io)
    str, date, time = split(header1)
    numrecord, dummy1 = split(header2)
    dt = DateTime("$date$(lpad(time,6,"0"))","yyyymmddHHMMSS")
    out = VarBC(dt,version,header1, header2,OrderedDict{String,VarBCRecord}())
    while !eof(io)
       push!(out.records, read(io,VarBCRecord))
    end
    close(io)
    return out
end

"""
    write(fname,a)

Write VarBC struct to file
"""
function write(fname::String,a::VarBC)
    io = open(fname,"w")
    write(io,a.version,"\n")
    write(io,a.header1,"\n")
    write(io,a.header2,"\n")
    i=1
    for (key,rec) in a.records
        write(io,rec,key,i)
        i=i+1
    end
    close(io)
end

"""
    DataFrame(a)

Constructs a DataFrame with param0 values from a VarBC struct.

Columnnames are:
datetime, key, p0, p1, ..... , p18
"""
function DataFrame(a::VarBC)
    colnames = ["datetime", "key" , "p".*string.(0:18)...]
    coltypes = [DateTime,  String, fill(Float64,19)...]
    named_tuple = (; zip(Symbol.(colnames), type[] for type in coltypes )...)
    df = DataFrames.DataFrame(named_tuple)
    for (key,val) in a.records
        par = fill(NaN,19)
        par[val.predcs.+1] = val.param0
        push!(df,[a.datetime,key,par...])
    end
    return df
end

end # module


text:
@enum MODE PRE=1 CURR=2 POST=3
struct TestContState
    x::Float64
    term::Bool
end
const TestStateType = CMSSPState{MODE,TestContState}

@enum MODESWITCH STAY=1 SWITCH=2

const TestActionType = CMSSPAction{MODESWITCH, Float64}

const TestCMSSPType = CMSSP{MODE,TestContState,MODESWITCH,Float64}

# Conversion functions
function POMDPs.convert_s(::Type{V}, s::TestStateType,
                          mdp::TestCMSSPType) where {V <: AbstractVector{Float64}}
    v = SVector{3,Float64}(convert(Float64, Int(s.mode)), s.continuous.x, convert(Float64, s.continuous.term))
    return v
end

function POMDPs.convert_s(::Type{TestStateType}, v::AbstractVector{Float64}, mdp::TestCMSSPType)
    mode = MODE(convert(Int64,v[1]))
    continuous = TestContState(v[2], convert(Bool,v[3]))
end


text:
#=
Copyright (c) 2015, Intel Corporation

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

    * Redistributions of source code must retain the above copyright notice,
      this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name of Intel Corporation nor the names of its contributors
      may be used to endorse or promote products derived from this software
      without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=#

export MaxPoolingLayer

@neuron type MaxNeuron
    maxidx :: Batch{Int}
end

MaxNeuron() = MaxNeuron(Batch(0))

@neuron forward(neuron::MaxNeuron) do
    begin
        maxval = -Inf32
        max_idx = 1
        for i in 1:length(neuron.inputs[1])
            curr = neuron.inputs[1][i]
            if curr > maxval
                maxval = curr
                max_idx = i
            end
        end
        neuron.value = maxval
        neuron.maxidx = max_idx
    end
end

@neuron backward(neuron::MaxNeuron) do
    begin
        idx = neuron.maxidx
        neuron.∇inputs[1][idx] += neuron.∇
    end
end

function MaxPoolingLayer(name::Symbol, net::Net, input_ensemble::AbstractEnsemble,
                         kernel::Int, stride::Int, pad::Int; ceil_mode=false)
    width, height, channels = size(input_ensemble)
    if ceil_mode
        # Round up
        div_op = cld
    else
        div_op = fld
    end
    width_out = div_op((width - kernel + 2 * pad), stride) + 1
    height_out = div_op((height - kernel + 2 * pad), stride) + 1
    neurons =
        [MaxNeuron() for _ = 1:width_out, _ = 1:height_out, _ = 1:channels]

    pool = Ensemble(net, name, neurons)
    add_connections(net, input_ensemble, pool, function (i, j, c)
        in_x = (i-1)*stride - pad
        in_y = (j-1)*stride - pad
        return (in_x+1:in_x+kernel, in_y+1:in_y+kernel, c:c)
    end; padding=pad)
    pool
end

export MeanPoolingLayer

@neuron type MeanNeuron
end

@neuron forward(neuron::MeanNeuron) do
    begin
        num_inputs = length(neuron.inputs[1])
        the_sum = 0.0f0
        for i in 1:num_inputs
            the_sum += neuron.inputs[1][i]
        end
        neuron.value = the_sum / num_inputs
    end
end

@neuron backward(neuron::MeanNeuron) do
    begin
        num_inputs = length(neuron.∇inputs[1])
        val = neuron.∇ / num_inputs
        for i in 1:num_inputs
            neuron.∇inputs[1][i] += val
        end
    end
end

function MeanPoolingLayer(name::Symbol, net::Net, input_ensemble::AbstractEnsemble,
                          kernel::Int, stride::Int, pad::Int; ceil_mode=false)
    if ceil_mode
        # Round up
        div_op = cld
    else
        div_op = fld
    end
    width, height, channels = size(input_ensemble)
    width_out = div_op((width - kernel + 2 * pad), stride) + 1
    height_out = div_op((height - kernel + 2 * pad), stride) + 1
    neurons =
        [MeanNeuron() for _ = 1:width_out, _ = 1:height_out, _ = 1:channels]

    pool = Ensemble(net, name, neurons)
    add_connections(net, input_ensemble, pool, function (i, j, c)
        in_x = (i-1)*stride - pad
        in_y = (j-1)*stride - pad
        return (in_x+1:in_x+kernel, in_y+1:in_y+kernel, c:c)
    end;padding=pad)
    pool
end


text:
seq = hofstQseq(1:100_000)
cnt = count(diff(seq) .< 0)
println("$cnt elements are less than the preceding one.")

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
@testset "reshape" begin
    A = rand(4, 5)
    B, pullback = rrule(reshape, A, (5, 4))
    @test B == reshape(A, (5, 4))
    Ȳ = randn(4, 5)

    (s̄, Ā, d̄) = pullback(Ȳ)
    @test s̄ == NO_FIELDS
    @test d̄ isa DoesNotExist
    @test extern(Ā) == reshape(Ȳ, (5, 4))

    B, pullback = rrule(reshape, A, 5, 4)
    @test B == reshape(A, 5, 4)

    Ȳ = randn(4, 5)
    (s̄, Ā, d̄1, d̄2) = pullback(Ȳ)
    @test s̄ == NO_FIELDS
    @test d̄1 isa DoesNotExist
    @test d̄2 isa DoesNotExist
    @test extern(Ā) == reshape(Ȳ, 5, 4)
end

@testset "hcat" begin
    A = randn(3, 2)
    B = randn(3)
    C = randn(3, 3)
    H, pullback = rrule(hcat, A, B, C)
    @test H == hcat(A, B, C)
    H̄ = randn(3, 6)
    (ds, dA, dB, dC) = pullback(H̄)
    @test ds == NO_FIELDS
    @test dA ≈ view(H̄, :, 1:2)
    @test dB ≈ view(H̄, :, 3)
    @test dC ≈ view(H̄, :, 4:6)
end

@testset "reduce hcat" begin
    A = randn(3, 2)
    B = randn(3, 1)
    C = randn(3, 3)
    x = [A, B, C]
    H, pullback = rrule(reduce, hcat, x)
    @test H == reduce(hcat, x)
    H̄ = randn(3, 6)
    x̄ = randn.(size.(x))
    rrule_test(reduce, H̄, (hcat, nothing), (x, x̄))
end

@testset "vcat" begin
    A = randn(2, 4)
    B = randn(1, 4)
    C = randn(3, 4)
    V, pullback = rrule(vcat, A, B, C)
    @test V == vcat(A, B, C)
    V̄ = randn(6, 4)
    (ds, dA, dB, dC) = pullback(V̄)
    @test ds == NO_FIELDS
    @test dA ≈ view(V̄, 1:2, :)
    @test dB ≈ view(V̄, 3:3, :)
    @test dC ≈ view(V̄, 4:6, :)
end

@testset "reduce vcat" begin
    A = randn(2, 4)
    B = randn(1, 4)
    C = randn(3, 4)
    x = [A, B, C]
    V, pullback = rrule(reduce, vcat, x)
    @test V == reduce(vcat, x)
    V̄ = randn(6, 4)
    x̄ = randn.(size.(x))
    rrule_test(reduce, V̄, (vcat, nothing), (x, x̄))
end

@testset "fill" begin
    y, pullback = rrule(fill, 44, 4)
    @test y == [44, 44, 44, 44]
    (ds, dv, dd) = pullback(ones(4))
    @test ds === NO_FIELDS
    @test dd isa DoesNotExist
    @test extern(dv) == 4

    y, pullback = rrule(fill, 2.0, (3, 3, 3))
    @test y == fill(2.0, (3, 3, 3))
    (ds, dv, dd) = pullback(ones(3, 3, 3))
    @test ds === NO_FIELDS
    @test dd isa DoesNotExist
    @test dv ≈ 27.0
end


text:
# not testable
# # this must be the test
# @testset "Eager evaluation" begin
#     reset_default_graph()
#     sess = Session()
#     enable_eager_execution()
#     g = rand(10,10)
#     o = Variable(g)
#     @test value(o) ≈ g
# end

@testset "control_dependency" begin
    a = Variable(ones(10))
    b = Variable(ones(10))
    for i = 1:10
        control_dependencies(a) do
            a = scatter_add(a, i, b[i])
        end
    end

    a_ = constant(0.0)
    b_ = Variable(0.0)
    op = assign(b_, 2.0)
    a_ = a_ + 1.0
    a_ = bind(a_, op)

    init(sess)
    run(sess, a_)
    @test run(sess, a)≈ones(10)*2
    @test run(sess, b_)≈2.0

end

@testset "while loop" begin
    a = constant(rand(10))
    condition = (i, var)->tf.less(i,11)
    function body(i, var)
        var = tf.cond(tf.equal(i, 1),
            ()->write(var, i, a),
            ()->write(var, i, 2*read(var, i-1))
        )
        i+1, var
    end
    i = constant(1, dtype=Int32)
    ta = TensorArray(10)
    out_i, out_ta = tf.while_loop(condition, body, [i,ta])
    ts = read(out_ta, 10)
    sum_ts = sum(ts)
    grd = gradients(sum_ts, a)
    @test run(sess, grd)≈512*ones(10)

    ta = TensorArray(10)
    a = constant(1.0)
    ta = write(ta, 1, a)
    ta = write(ta, 2, read(ta,1)+a) # 2a
    ta = write(ta, 3, read(ta,2)+a) # 3a
    ta = write(ta, 4, read(ta,2)+read(ta,3)+a) # 6a
    g = gradients(read(ta, 4), a)
    @test run(sess, g)≈6

    ta = TensorArray(9)

    function condition(i, ta)
        tf.less(i, 10)
    end
    ta = write(ta, 1, constant(1.0))
    ta = write(ta, 2, constant(1.0))
    function body(i, ta)
        ta = write(ta, i, 2*read(ta,i-1)+read(ta,i-2)+a)
        i+1, ta
    end

    i = constant(3, dtype=Int32)
    _, ta_out = while_loop(condition, body, [i,ta])
    g = gradients(read(ta_out,9),a)
    @test run(sess, g)≈288

end

@testset "if_clause" begin
    a = constant(2.0)
    b = constant(1.0)
    c = constant(0.0)
    bb = b*2
    cc = c + 3.0
    res = if_else(a>1.0, bb, cc)
    @test run(sess, res)≈2.0
end

@testset "if_else: tf.where" begin
    condition = [true true false false
                true true true true]
    a = [1 1 1 1
         1 1 1 1]
    b = [2 2 2 2
         2 2 2 2]
    res = if_else(condition, a, b)
    @test run(sess, res)≈[1 1 2 2
                          1 1 1 1]

    res1 = if_else(condition[1,:], a[1,:], b[1,:])
    @test run(sess, res1)≈[1;1;2;2]
end

@testset "get and add collection" begin
    a = Variable(ones(10), name="my_collect1")
    b = 2a
    @test get_collection("my_collect")==[a]
end

@testset "has_gpu" begin
    @test has_gpu()==false
end

@testset "timeline" begin
    a = normal(2000, 5000)
    b = normal(5000, 1000)
    res = a*b
    run_profile(sess, res)
    save_profile("test.json")
    rm("test.json")
end

@testset "independent" begin
    x = constant(rand(10))
    y = 2*x
    z = sum(independent(y))
    @test isnothing(gradients(z, x))
end


text:
using Base.Test
using Reactive

facts("Flatten") do

    a = Signal(0)
    b = Signal(1)

    c = Signal(a)

    d = flatten(c)
    cnt = foldp((x, y) -> x+1, 0, d)

    context("Signal{Signal} -> flat Signal") do
        # Flatten implies:
        @fact value(c) --> a
        @fact value(d) --> value(a)
    end

    context("Initial update count") do

        @fact value(cnt) --> 0
    end

    context("Current signal updates") do
        push!(a, 2)
        step()

        @fact value(cnt) --> 1
        @fact value(d) --> value(a)
    end

    context("Signal swap") do
        push!(c, b)
        step()
        @fact value(cnt) --> 2
        @fact value(d) --> value(b)

        push!(a, 3)
        step()
        @fact value(cnt) --> 2
        @fact value(d) --> value(b)

        push!(b, 3)
        step()

        @fact value(cnt) --> 3
        @fact value(d) --> value(b)
    end
end


text:
module GrandGameOfLife

using Random

export LifeCells, makestep, create_initial_state

struct Cell
    x::Int32
    y::Int32
end

Base.:(+)(c::Cell, d::Tuple) = Cell(c.x + d[1], c.y + d[2])

struct LifeCells
    bag::Set{Cell}
    bag2::Set{Cell}
end

LifeCells(v) = LifeCells(Set(v), Set(v))

has(L::LifeCells, c::Cell) = c in L.bag

# These are only needed for debugging purposes
# function to_xy(c::UInt64)
#     (c >> 32) % Int, (c & 0x00000000ffffffff) % Int
# end

# function to_cell(x::Int, y::Int)
#     return (x << 32 | y) % UInt64
# end

# function to_cell(x::Tuple{Int, Int})
#     return to_cell(x[1], x[2])
# end

struct Neighbourhood
    c::Cell
end

function Base.iterate(iter::Neighbourhood, state = 1)
    if state == 1
        return (iter.c + (-1, -1), 2)
    elseif state == 2
        return (iter.c + (-1, 0), 3)
    elseif state == 3
        return (iter.c + (-1, 1), 4)
    elseif state == 4
        return (iter.c + (0, -1), 5)
    elseif state == 5
        return (iter.c + (0, 1), 6)
    elseif state == 6
        return (iter.c + (1, -1), 7)
    elseif state == 7
        return (iter.c + (1, 0), 8)
    elseif state == 8
        return (iter.c + (1, 1), 9)
    else
        return nothing
    end
end

# function neighbours!(neighbours_container, big_neighbourhood_container, L::LifeCells)
#     idx = 0
#     for coord in big_neighbourhood_container
#         if coord in L.bag
#             idx += 1
#             neighbours_container.v[idx] = coord
#         end
#     end
#     neighbours_container.l = idx
#     neighbours_container
# end
#
# function is_in_neighbours(neighbours_container, cell)
#     for i in 1:neighbours_container.l
#         neighbours_container.v[i] == cell && return 1
#     end
#     return 0
# end

function makestep(L::LifeCells)
    # initialize bag with the original size
    bag = L.bag2
    empty!(bag)
    sizehint!(bag, length(L.bag))

    for cell in L.bag
        cnt = 0
        for c in Neighbourhood(cell)
            k = has(L, c)
            cnt += k
            k && continue
            c in bag && continue

            cnt2 = 0
            for c2 in Neighbourhood(c)
                cnt2 += has(L, c2)
            end
            if cnt2 == 3
                push!(bag, c)
            end
        end
        if (cnt == 2) | (cnt == 3)
            push!(bag, cell)
        end
    end

    return LifeCells(bag, L.bag)
end

function create_initial_state(n, m, density = 0.7, rng = Random.GLOBAL_RNG)
    bag = Set{Cell}()
    sizehint!(bag, Int(ceil(n*m*density)))
    for i in 1:n, j in 1:m
        if rand(rng) <= density
            push!(bag, Cell(i, j))
        end
    end

    LifeCells(bag)
end

end # module


text:
# ===============================
# First-order approximation model
# ===============================

"""
    FirstOrder{EM} <: AbstractApproximationModel

First-order approximation model.

### Fields

- `exp` -- (optional, default: `BaseExp`) exponentiation method

### Algorithm

The transformations are:

- ``Φ ← \\exp(Aδ)``,
- ``Ω_0 ← CH(\\mathcal{X}_0, Φ\\mathcal{X}_0 ⊕ δU ⊕ B_ε)``

where ``B_ε`` is the input ball of radius ``ε`` centered in the origin.

### Reference

C. Le Guernic and A. Girard: Reachability analysis of linear systems using
support functions. NAHS 2010.
"""
struct FirstOrder{EM} <: AbstractApproximationModel
    exp::EM
end

# convenience constructor
function FirstOrder(; exp=BaseExp)
    return FirstOrder(exp)
end

function Base.show(io::IO, alg::FirstOrder)
    print(io, "`FirstOrder` approximation model with: \n")
    print(io, "    - exponentiation method: $(alg.exp) \n")
end

Base.show(io::IO, m::MIME"text/plain", alg::FirstOrder) = print(io, alg)

# -----------------------------------------------
# FirstOrder approximation: Homogeneous case
# -----------------------------------------------

function discretize(ivp::IVP{<:CLCS, <:AbstractZonotope}, δ, alg::FirstOrder)
    A = state_matrix(ivp)
    X0 = initial_state(ivp)
    n = dim(X0)
    norm_A = opnorm(A, Inf)

    # matrix exponential
    Φ = _exp(A, δ, alg.exp)

    # compute Ω0
    Xδ = Φ * X0
    α = (exp(norm_A * δ) - 1 - norm_A * δ) * norm(X0, Inf)
    Ω0 = ConvexHull(X0, Xδ + BallInf(zeros(n), α))

    # create result
    X = stateset(ivp)
    Sdis = ConstrainedLinearDiscreteSystem(Φ, X)
    return InitialValueProblem(Sdis, Ω0)
end

# -------------------------------------------------
# FirstOrder approximation: Inhomogeneous case
# -------------------------------------------------

function discretize(ivp::IVP{<:CLCCS, <:AbstractZonotope}, δ, alg::FirstOrder)
    A = state_matrix(ivp)
    X0 = initial_state(ivp)
    n = dim(X0)
    norm_A = opnorm(A, Inf)

    # matrix exponential
    Φ = _exp(A, δ, alg.exp)

    # compute Ω0
    Xδ = Φ * X0
    U = next_set(inputset(ivp), 1)
    factor = (exp(norm_A * δ) - 1 - norm_A * δ)
    norm_U_over_A = norm(U, Inf) / norm_A
    α = factor * (norm(X0, Inf) + norm_U_over_A)
    Ω0 = ConvexHull(X0, Xδ + δ * U + BallInf(zeros(n), α))

    # discretize inputs
    β = factor * norm_U_over_A
    V = δ * U + BallInf(zeros(n), β)

    # create result
    B = IdentityMultiple(one(eltype(A)), size(A, 1))
    X = stateset(ivp)
    Sdis = ConstrainedLinearControlDiscreteSystem(Φ, B, X, V)
    return InitialValueProblem(Sdis, Ω0)
end


text:
@testset "JONSWAP" begin
    θ₀ = [0.7, 0.9, 3.3, 5.0]
    ω = 1.0

    @testset "Gradient" begin
        @test approx_gradient_uni(θ -> sdf(JONSWAP{1}(θ), ω), θ₀)      ≈ grad_sdf(JONSWAP{1}(θ₀), ω)
    end
    @testset "Hessian" begin
        @test approx_hessian_uni( θ -> grad_sdf(JONSWAP{1}(θ), ω), θ₀) ≈ hess_sdf(JONSWAP{1}(θ₀), ω)
    end

    @testset "Error handling" begin
        @test_throws ArgumentError JONSWAP{1}(-0.7, 0.9, 3.3, 5.0)
        @test_throws ArgumentError JONSWAP{1}(0.7, -0.9, 3.3, 5.0)
        @test_throws ArgumentError JONSWAP{1}(0.7, 0.9, 0.8, 5.0)
        @test_throws ArgumentError JONSWAP{1}(0.7, 0.9, 3.3, 0.9)
        @test_throws ArgumentError JONSWAP{1}(ones(3))
        @test_throws ArgumentError JONSWAP{1}(ones(5))
        @test_throws ErrorException JONSWAP(ones(4))
        @test_throws ErrorException JONSWAP(0.7, 0.9, 3.3, 5.0)
    end
end


text:
using ArgParse
using Crayons.Box
using DataStructures
import Statistics, StatsBase
using NumericalDataManipulation.CommonMath
using NumericalDataManipulation.Data.Manipulation
using NumericalDataManipulation.Data.Manager

# ---- system ----
Base.eval(:(have_color = true))
function Base.show(io::IO, intervals::Vector{Tuple{Float64, Float64, Symbol}})
    print(io, "{$(join(map(ivl -> "$ivl", intervals), ", "))}")
end
function Base.show(io::IO, ::MIME"text/plain", intervals::Vector{Tuple{Float64, Float64, Symbol}})
    show(io, intervals)
end
# ----------------

function main(args::Vector{String})
    try
        arguments = parse_commandline(args)
        work(arguments)
    catch error
        error_type = BOLD(UNDERLINE(string(typeof(error))))
        if isa(error, ArgParseError) || isa(error, ArgumentError)
            if isa(error, ArgParseError)
            @error "$error_type ⭃ $(error.msg)"
            elseif isa(error, ArgumentError)
                @error "$error_type ⭃ $(error.msg):\n$(sprint(io->Base.show_backtrace(io, catch_backtrace())))"
            end
        else
            #@error "$error_type ⭃ $(error)\n$(BOLD("𝐓race")):\n$(stacktrace(catch_backtrace()))"
            rethrow(error)
        end
    end
end
function parse_commandline(args::Vector{String})
    return arguments = parse_args(args, ARG_SETTINGS)
end
function work(arguments::Dict{String, Any})
    @assert haskey(arguments, "%COMMAND%")
    command = arguments["%COMMAND%"]
    options = arguments[command]
    if "merge-tables" == command
        file_a = options["file-a"]; file_a_skip_title = options["file-a-skip-title"]
        file_b = options["file-b"]; file_b_skip_title = options["file-b-skip-title"]
        master = options["master-file"]
        columns = options["columns"]
        intervals = options["intervals"]
        merge_function_strategy = options["function-merge-strategy"]
        grid_merge_parameters = options["grid-merge-parameters"]
        @info """Data table selective merge:
                 ✔ File A: '$file_a';
                 ✔ File B: '$file_b';
                 ✔ Columns: $columns;
                 ✔ Intervals: {$(join(map(ivl -> "[$(ivl[1]), $(ivl[2])]:$(ivl[3])", intervals), ", "))};
                 ✔ Function merge strategy: $merge_function_strategy;
                 ✔ Grid merge parameters: $grid_merge_parameters"""
        master_interval, slave_intervals = normalize_intervals(intervals, master)
        @info "Result: master - $master_interval; slaves - $slave_intervals"
        taks = MergeTwoTablesTask(file_a, file_b,
            master, columns, master_interval, slave_intervals,
            merge_function_strategy, grid_merge_parameters,
            file_a_skip_title, file_b_skip_title)
        Manager.merge_two_tables(taks)
    elseif "dummy-command" == command
        @info "Dummy command with options: $options"
    end
end

# --- utilities ---
function normalize_intervals(intervals::Vector{Tuple{Float64, Float64, Symbol}}, master::Symbol)
    # TODO: maybe a smarter approach is possible
    @assert !isempty(intervals)
    check_conflicts(intervals, master)
    normalized_intervals = Vector{Tuple{Float64, Float64, Symbol}}()

    master_interval = ifempty(filter(ivl -> master == ivl[3], intervals), [(-Inf, Inf, :a)])[1]
    slave = master == :a ? :b : :a
    slave_intervals = sort(filter(ivl -> slave == ivl[3], intervals); by = ivl -> ivl[1])

    return (master_interval, slave_intervals)
end
function check_conflicts(intervals::Vector{Tuple{Float64, Float64, Symbol}}, master::Symbol)
    # ---- duplicates
    conflicts = find_interval_file_duplicates(intervals)
    if !isempty(conflicts)
        sb = IOBuffer()
        println(sb, "Conflicting intervals are found:")
        for conflict ∈ conflicts
            civl = conflict.first
            file_ivls = conflict.second
            println(sb, "  - Conflict interval $civl: $file_ivls")
        end
        throw(ArgumentError(String(take!(sb))))
    end

    # ---- conflicts of master and slave intervals
    slave = master == :a ? :b : :a
    master_intervals = filter(ivl -> master == ivl[3], intervals)
    slave_intervals = filter(ivl -> slave == ivl[3], intervals)
    if length(master_intervals) > 1
        throw(ArgumentError("Only one master file \"$master\" interval can be specified. Given: $master_intervals"))
    end
    if isempty(slave_intervals)
        throw(ArgumentError("No intervals for the slave file \"$slave\" are specified."))
    end
    if length(master_intervals) == 1
        m_interval = master_intervals[1]
        conflicting_slaves = filter(s_ivl -> s_ivl[1] < m_interval[1] || s_ivl[2] > m_interval[2], slave_intervals)
        if !isempty(conflicting_slaves)
            throw(ArgumentError("There are \"$slave\" file intervals which are out of the master interval $m_interval: $conflicting_slaves"))
        end
    end

    # ---- overlap conflicts between slave intervals
    conflicts = find_overlapping_conflicts(slave_intervals)
    if !isempty(conflicts)
        sb = IOBuffer()
        println(sb, "Conflicting slave intervals are found:")
        for conflict ∈ conflicts
            civl = conflict.first
            file_ivls = conflict.second
            println(sb, "  - Conflict for the slave interval $civl: $file_ivls")
        end
        throw(ArgumentError(String(take!(sb))))
    end
end
function find_overlapping_conflicts(slave_intervals::Vector{Tuple{Float64, Float64, Symbol}})
    dictionary = OrderedDict{Tuple{Float64, Float64}, Vector{Tuple{Float64, Float64, Symbol}}}()
    for (i, ivl_i) ∈ enumerate(slave_intervals)
        ivl_i_conflicts =
            filter(it -> begin j, ivl_j = it
                    i ≠ j && (ivl_i[1] ≤ ivl_j[1] ≤ ivl_i[2] || ivl_i[1] ≤ ivl_j[2] ≤ ivl_i[2])
                end,
                collect(enumerate(slave_intervals)))
        if !isempty(ivl_i_conflicts)
            dictionary[(ivl_i[1], ivl_i[2])] = map(it -> it[2], ivl_i_conflicts)
        end
    end
    return dictionary
end
function find_interval_file_duplicates(intervals::Vector{Tuple{Float64, Float64, Symbol}})
    dictionary = OrderedDict{Tuple{Float64, Float64}, Vector{Tuple{Float64, Float64, Symbol}}}()
    for ivl ∈ intervals
        if !haskey(dictionary, (ivl[1], ivl[2]))
            dictionary[(ivl[1], ivl[2])] = [ivl]
        else
            push!(dictionary[(ivl[1], ivl[2])], ivl)
        end
    end
    duplicates = OrderedDict{Tuple{Float64, Float64}, Vector{Tuple{Float64, Float64, Symbol}}}()
    for entry ∈ dictionary
        if length(entry.second) > 1
            duplicates[entry.first] = entry.second
        end
    end
    return duplicates
end
function get_sorted_intervals(intervals::Vector{Tuple{Float64, Float64, Symbol}}, by::Symbol)
    @assert by ∈ [:left, :right]
    left_sorted = sort(intervals; by = ivl -> ivl[1])
    dictionary = OrderedDict{Float64, Vector{Tuple{Float64, Float64, Symbol}}}()
    keyfunction = by == :left ? ivl -> ivl[1] : ivl -> ivl[2]
    for ivl ∈ intervals
        if !haskey(dictionary, keyfunction(ivl))
            dictionary[keyfunction(ivl)] = [ivl]
        else
            push!(dictionary[keyfunction(ivl)], ivl)
        end
    end
    return dictionary
end
# -----------------
# --- common ---
function ifempty(object::Vector, default::Vector)
    return isempty(object) ? default : object
end
# --------------

# --------------
function ArgParse.parse_item(::Type{MergeFunctionStrategyName}, some_strategy::AbstractString)
    if !isempty(some_strategy) && some_strategy ∈ ["naive-join", "sigmoid-join"]
        throw(ArgParseError("Wrong merge function strategy name '$some_strategy'; Supportred values are 'naive-join' and 'sigmoid-join'."))
    end
    if "naive-join" == some_strategy
        return NAIVE_JOIN::MergeFunctionStrategyName
    elseif "sigmoid-join" == some_strategy
        return SIGMOID_JOIN::MergeFunctionStrategyName
    else
        throw(ArgParseError("Unexpected merge function strategy name '$some_strategy'"))
    end
end

const REGEX_GRID_PARS = r"^\d+\s*,\s*\d+(\.\d+)?\s*,\s*(mean_log10_golden_lower|mean_log10|geomean|harmmean|mean)$"
const REGEX_SPLIT_GRID_PARS = r"\s*,\s*"
function ArgParse.parse_item(::Type{GridMergeParameters}, parameters::AbstractString)
    if !occursin(REGEX_GRID_PARS, parameters)
        throw(ArgParseError("Wrong grid merge parameter string: $parameters"))
    end
    tail = parse(Int, match(r"\d+(?=\s*,\s*\d+(\.\d+)?)", parameters).match)
    ϵ = parse(Float64, match(r"\d+(\.\d+)?(?=\s*,\s*(mean_log10_golden_lower|mean_log10|geomean|harmmean|mean))", parameters).match)

    mean_function_string = match(r"mean_log10_golden_lower|mean_log10|mean", parameters).match
    mean_function = if "mean_log10" == mean_function_string
        mean_log10
    elseif "mean_log10_golden_lower" == mean_function_string
        mean_log10_golden_lower
    elseif "geomean" == mean_function_string
        StatsBase.geomean
    elseif "harmmean" == mean_function_string
        StatsBase.harmmean
    elseif "mean" == mean_function_string
        Statistics.mean
    end

    return GridMergeParameters(tail, ϵ, mean_function)
end

function ArgParse.parse_item(::Type{Symbol}, str::AbstractString)
    master_file = Symbol(str)
    if master_file ∉ [:a, :b]
        throw(ArgParseError("Choose the value \"a\" or \"b\" to specify the master file. Given value: $master_file"))
    end
    return master_file
end

const REGEX_COLUMNS = r"^([1-9][0-9]*\s*,\s*)*([1-9][0-9]*)$"
const REGEX_SPLIT_COLUMNS = r"\s*,\s*"
function ArgParse.parse_item(::Type{Vector{Int}}, columns_string::AbstractString)
    columns = strip(columns_string)
    if !occursin(REGEX_COLUMNS, columns)
        throw(ArgParseError("The column list is not a comma separated list of integers."))
    end
    column_numbers = map(c -> parse(Int, c), split(columns, REGEX_SPLIT_COLUMNS))
    if !allunique(column_numbers)
        throw(ArgParseError("The column list $column_numbers contains duplicates."))
    end
    return sort(column_numbers)
end

function ArgParse.parse_item(::Type{Symbol}, str::AbstractString)
    master_file = Symbol(str)
    if master_file ∉ [:a, :b]
        throw(ArgParseError("Choose the value \"a\" or \"b\" to specify the master file. Given value: $master_file"))
    end
    return master_file
end

const REGEX_INTERVALS = r"^(\[([0-9]+(.[0-9]+)?|Inf|-Inf)\s*,\s*([0-9]+(.[0-9]+)?|Inf|-Inf)\](:a|:b)\s*,\s*)*\[([0-9]+(.[0-9]+)?|Inf|-Inf)\s*,\s*([0-9]+(.[0-9]+)?|Inf|-Inf)\](:a|:b)$"
const REGEX_SPLIT_INTERVALS = r"(?<=\](:a|:b))\s*,\s*(?=\[)"
const REGEX_INTERVAL_LEFT = r"(?<=\[)([0-9]+(.[0-9]+)?|Inf|-Inf)(?=\s*,)"
const REGEX_INTERVAL_RIGHT = r"([0-9]+(.[0-9]+)?|Inf|-Inf)(?=\])"
const REGEX_INTERVAL_ORIGIN = r"(?<=\]:)(a|b)"
function ArgParse.parse_item(::Type{Vector{Tuple{Float64, Float64, Symbol}}}, intervals_string::AbstractString)
    intervals = strip(intervals_string)
    if !occursin(REGEX_INTERVALS, intervals)
        throw(ArgParseError("The given interval list is not a comma separated list of intervals. Example: \"[1.0, 2.0]:a, [29.6, 29.7]:b\"."))
    end
    interval_list = map(interval -> begin
            left = parse(Float64, match(REGEX_INTERVAL_LEFT, interval).match)
            right = parse(Float64, match(REGEX_INTERVAL_RIGHT, interval).match)
            origin = Symbol(match(REGEX_INTERVAL_ORIGIN, interval).match)
            if left == Inf
                throw(ArgParseError("The left border of the interval [$left, $right]:$origin cannot be Inf."))
            end
            if right == Inf
                throw(ArgParseError("The right border of the interval [$left, $right]:$origin cannot be Inf."))
            end
            return (left, right, origin)
        end,
        split(intervals, REGEX_SPLIT_INTERVALS))
    return interval_list
end

# ---- CLI CONFIGURATION ----
# --- LOGO ---
const LOGO = strip(
"""
╔╦╗┌─┐┌┐┌┬┌─┐┬ ┬┬  ┌─┐┌┬┐┬┌─┐┌┐┌
║║║├─┤││││├─┘│ ││  ├─┤ │ ││ ││││
╩ ╩┴ ┴┘└┘┴┴  └─┘┴─┘┴ ┴ ┴ ┴└─┘┘└┘
""")
# --- CLI settings and main() call ---
const ARG_SETTINGS = ArgParseSettings(;prog = "manipulation")
@add_arg_table ARG_SETTINGS begin
    "merge-tables"
        help = "Merge two tables from two given files."
        action = :command
    "dummy-command"
        help = "Dummy command"
        action = :command
end
# ---- @add_arg_table ARG_SETTINGS["merge-tables"]
@add_arg_table ARG_SETTINGS["merge-tables"] begin
    "--file-a", "-a"
        help = "Data file A. Supported file extension: '.csv'"
        arg_type = String
        required = true
    "--file-b", "-b"
        help = "Data file B. Supported file extension: '.csv'"
        arg_type = String
        required = true
    "--file-a-skip-title"
        help = "Skip the first line (usually a title) in the file A."
        arg_type = Bool
        default = true
        required = false
    "--file-b-skip-title"
        help = "Skip the first line (usually a title) in the file B."
        arg_type = Bool
        default = true
        required = false
    "--master-file", "-m"
        help = "Choose file \"a\" or \"b\" as a \"master\" data file to merge all data to."
        arg_type = Symbol
        default = :a
        required = false
    "--columns", "-c"
        help = """Select the column numbers to merge.
                  A comma separated list of integers.
                  The first column in files is always an argument grid (X).
                  $(BOLD("Always")) start your column numbering (1, 2, ...)
                  from the second column.
                  TODO: Support tables with different columns numbers."""
        arg_type = Vector{Int}
        required = true
    "--intervals", "-i"
        help = """Comma separated list of intervals in square brackets in quotes.
                  Example: \"[1.0, 2.0]:a, [29.6, 29.7]:b\".
                  It it allowed to use symbols Inf and -Inf in the interval values.
                  The intervals will be simplified automatically if possible.
                  The strings :a and :b denote which file must be used to
                  import data in the given interval."""
        arg_type = Vector{Tuple{Float64, Float64, Symbol}}
        required = true
    "--function-merge-strategy", "-f"
        help = """The function merge strategy name denotes a way of joining
                  the pieces of two similar functions from different tables.
                  Supported strategies:
                  $(UNDERLINE("✔ naive-join")):
                  Simple interval join of the table functions;
                  $(UNDERLINE("✔ sigmoid-join")):
                  The intervals are joined
                  with the sigmoid function with a rule
                  $(ITALICS("'f(x) = (1 - σ(x; x₀, α))⋅f⁽¹⁾(x) + σ(x; x₀, α)⋅f⁽²⁾(x)'")),
                  where σ is the sigmoid function,
                  parameter x₀ defines a boundary point between
                  interval 'a' and interval 'b',
                  parameter α defines the sigmoid function area width,
                  the functions f⁽¹⁾ and f⁽²⁾ originate from
                  $(BOLD("file-a")) and $(BOLD("file-b"))."""
        arg_type = MergeFunctionStrategyName
        default = SIGMOID_JOIN::MergeFunctionStrategyName
        required = false
    "--grid-merge-parameters", "-g"
        help = """Grid merging parameters (comma separated values):
                  1. Number of knots to merge the small intervals;
                  Default value - 5;
                  2. Relative smallness ϵʳᵉˡ of a grid interval in comparisson
                  to a mean grid value (usually logarithmic mean);
                  Default value - 0.4;
                  3. Function name to calculate grid step mean value;
                  Supported values:
                  mean_log10,
                  mean_log10_golden_lower,
                  geomean,
                  harmmean,
                  mean;
                  Default value - mean_log10
                  (mean_log10(Δx₁, Δx₂) = ((υ₁ + υ₂) / 2)^((o₁ + o₂) / 2),
                  where Δx = υ × 10^o)
                  """
        arg_type = GridMergeParameters
        default = GridMergeParameters(5, 0.4, mean_log10)
        required = false
end

@add_arg_table ARG_SETTINGS["dummy-command"] begin
    "--dummy-option"
        help = "Data file A"
        arg_type = String
        required = false
end

# ------------------------------------
if !isinteractive()
    println(RED_FG(LOGO))
    main(ARGS)
end
# ------------------------------------


text:
@testset "Attributes" begin
    @testset "is_set_by_optimize" begin
        @test MOI.is_set_by_optimize(MOI.TerminationStatus())
        @test !MOI.is_set_by_optimize(MOI.ConstraintSet())
        @test !MOI.is_set_by_optimize(MOI.ObjectiveSense())
    end
    @testset "is_copyable" begin
        @test !MOI.is_copyable(MOI.TerminationStatus())
        @test !MOI.is_copyable(MOI.ConstraintSet())
        @test MOI.is_copyable(MOI.ObjectiveSense())
    end
    @testset "supports" begin
        model = DummyModel()
        @test_throws ArgumentError MOI.supports(model, MOI.TerminationStatus())
        @test_throws ArgumentError begin
            MOI.supports(model, MOI.ConstraintSet(),
                         MOI.ConstraintIndex{MOI.SingleVariable,
                                             MOI.EqualTo{Float64}})
        end
        @test MOI.supports(model, MOI.ObjectiveSense())
    end
end


text:
include("header.jl")

# http://docs.julialang.org/en/latest/manual/arrays.html#man-supported-index-types-1

# Test KnetArray operations: cat, convert, copy, display, eachindex,
# eltype, endof, fill!, first, getindex, hcat, isempty, length,
# linearindexing, ndims, ones, pointer, rand!, reshape, setindex!,
# similar, size, stride, strides, summary, vcat, vec, zeros

if gpu() >= 0
    @testset "karray" begin
        a = rand(3,4)
        k = KnetArray(a)

        # getindex, setindex!
        # Index types: Integer, CartesianIndex, Vector{Int}, Array{Int}, EmptyArray, a:c, a:b:c, Colon, Bool
        # See http://docs.julialang.org/en/latest/manual/arrays.html#man-supported-index-types-1
        # check out http://docs.julialang.org/en/latest/manual/arrays.html#Cartesian-indices-1
        @testset "indexing" begin
            @test a == k                     		# Supported index types:
            for i in ((:,), (:,:),                      # Colon, Tuple{Colon}
                      (3,), (2,3),              	# Int, Tuple{Int}
                      (3:5,), (1:2,3:4),                # UnitRange, Tuple{UnitRange}
                      (2,:), (:,2),                     # Int, Colon
                      (1:2,:), (:,1:2),                 # UnitRange,Colon
                      (1:2,2), (2,1:2),                 # Int, UnitRange
                      (1:2:3,),                         # StepRange
                      (1:2:3,:), (:,1:2:3),             # StepRange,Colon
                      ([1,3],), ([2,2],),               # Vector{Int}
                      ([1,3],:), (:,[1,3]),             # Vector{Int},Colon
                      ([2,2],:), (:,[2,2]),             # Repeated index
                      ([],),                            # Empty Array
                      ((a.>0.5),),                      # BitArray
                      ([1 3; 2 4],),                    # Array{Int}
                      (CartesianIndex(3,),), (CartesianIndex(2,3),), # CartesianIndex
                      (:,a[1,:].>0.5),(a[:,1].>0.5,:),  # BitArray2 # FAIL for julia4
                      ([CartesianIndex(2,2), CartesianIndex(2,1)],), # Array{CartesianIndex} # FAIL for julia4
                      )
                #@show i
                @test a[i...] == k[i...]
                ai = a[i...]
                if isa(ai, Number)
                    a[i...] = 0
                    k[i...] = 0
                    @test a == k
                    a[i...] = ai
                    k[i...] = ai
                else
                    a[i...] .= 0
                    k[i...] .= 0
                    @test a == k
                    a[i...] .= ai
                    k[i...] .= ai
                end
                @test a == k
                @test gradcheck(getindex, a, i...; args=1)
                @test gradcheck(getindex, k, i...; args=1)
            end
            # make sure end works
            @test a[2:end] == k[2:end]
            @test a[2:end,2:end] == k[2:end,2:end]
            # k.>0.5 returns KnetArray{T}, no Knet BitArrays yet
            #TODO: @test a[a.>0.5] == k[k.>0.5]
        end

        # Unsupported indexing etc.:
        # @test_broken a[1:2:3,1:3:4] == Array(k[1:2:3,1:3:4]) # MethodError: no method matching getindex(::Knet.KnetArray{Float64,2}, ::StepRange{Int64,Int64}, ::StepRange{Int64,Int64})
        # @test_broken a[[3,1],[4,2]] == Array(k[[3,1],[4,2]]) # MethodError: no method matching getindex(::Knet.KnetArray{Float64,2}, ::Array{Int64,1}, ::Array{Int64,1})
        # @test_broken cat((1,2),a,a) == Array(cat((1,2),k,k)) # cat only impl for i=1,2

        a = rand(3,4)
        k = KnetArray(a)

        # AbstractArray interface
        @testset "abstractarray" begin

            for f in (copy, lastindex, first, isempty, length, ndims, vec, zero,
                      a->(eachindex(a);0), a->(eltype(a);0), # a->(Base.linearindexing(a);0),
                      a->collect(Float64,size(a)), a->collect(Float64,strides(a)),
                      a->cat(a,a;dims=1), a->cat(a,a;dims=2), a->hcat(a,a), a->vcat(a,a),
                      a->reshape(a,2,6), a->reshape(a,(2,6)),
                      a->size(a,1), a->size(a,2),
                      a->stride(a,1), a->stride(a,2), )

                #@show f
                @test f(a) == f(k)
                @test gradcheck(f, a)
                @test gradcheck(f, k)
            end

            @test convert(Array{Float32},a) == convert(KnetArray{Float32},k)
            @test fill!(similar(a),pi) == fill!(similar(k),pi)
            @test fill!(similar(a,(2,6)),pi) == fill!(similar(k,(2,6)),pi)
            @test fill!(similar(a,2,6),pi) == fill!(similar(k,2,6),pi)
            @test isa(pointer(k), Ptr{Float64})
            @test isa(pointer(k,3), Ptr{Float64})
            @test isempty(KnetArray{Float32}(undef,0))
            @test rand!(copy(a)) != rand!(copy(k))
            @test k == k
            @test a == k
            @test k == a
            @test isapprox(k,k)
            @test isapprox(a,k)
            @test isapprox(k,a)
            @test a == copyto!(similar(a),k)
            @test k == copyto!(similar(k),a)
            @test k == copyto!(similar(k),k)
            @test k == copy(k)
            @test pointer(k) != pointer(copy(k))
            @test k == deepcopy(k)
            @test pointer(k) != pointer(deepcopy(k))
        end

        a = rand(3,4)
        k = KnetArray(a)

        @testset "cpu2gpu" begin
            # cpu/gpu xfer with grad support
            @test gradcheck(x->Array(sin.(KnetArray(x))),a)
            @test gradcheck(x->KnetArray(sin.(Array(x))),k)
        end

        a = rand(3,4)
        k = KnetArray(a)

        @testset "reshape" begin
            a = KnetArray{Float32}(undef, 2, 2, 2)

            @test size(reshape(a, 4, :)) == size(reshape(a, (4, :))) == (4, 2)
            @test size(reshape(a, :, 4)) == size(reshape(a, (:, 4))) == (2, 4)
            @test size(reshape(a, :, 1, 4)) == (2, 1,  4)
        end
    end
end

nothing


text:
# Import packages
using  OrdinaryDiffEq, DiffEqDevTools, Test, Random

# Define test problems
# Note that the time span in DiffEqProblemLibrary.ODEProblemLibrary is given by
# Float64 numbers

linear = (u,p,t) -> (p*u)
linear_analytic = (u0,p,t) -> u0*exp(p*t)
prob_ode_bigfloatlinear = ODEProblem(
                          ODEFunction(linear,analytic=linear_analytic),
                          big"0.5",(big"0.0",big"1.0"),big"1.01")

f_2dlinear = (du,u,p,t) -> (@. du = p*u)
f_2dlinear_analytic = (u0,p,t) -> @. u0*exp(p*t)
prob_ode_bigfloat2Dlinear = ODEProblem(
                    ODEFunction(f_2dlinear,analytic=f_2dlinear_analytic),
                  rand(BigFloat,(4,2)),(big"0.0",big"1.0"),big"1.01")

# Prepare tests
Random.seed!(100)
problem_array = [prob_ode_bigfloatlinear,prob_ode_bigfloat2Dlinear]
dts = 1 .//2 .^(8:-1:1)

testTol = 0.2

@testset "Testing extrapolation methods" begin

    # Test AitkenNeville
    println("Testing AitkenNeville")
    @testset "Testing AitkenNeville" begin
        @testset "Testing sequential AitkenNeville" begin
            for prob in problem_array
            global dts

            #  Convergence test
            for j = 1:4
                sim = test_convergence(dts,prob,AitkenNeville(max_order = j,
                min_order = j, init_order = j, threading = false))
                @test sim.𝒪est[:final] ≈ j atol=testTol
            end

            # Regression test
            sol = solve(prob,AitkenNeville(max_order = 9, min_order = 1,
                init_order = 9, threading = false),reltol=1e-3)
            @test length(sol.u) < 15
            sol = solve(prob,AitkenNeville(max_order = 9, min_order = 1,
                init_order = 9, threading = false),reltol=1e-6)
            @test length(sol.u) < 18
            end
        end
    end # AitkenNeville

    # Define the subdividing sequences
    sequence_array =[:harmonic, :romberg, :bulirsch]

    println("Testing ImplicitEulerExtrapolation")
    @testset "Testing ImplicitEulerExtrapolation" begin
    for prob in problem_array, seq in sequence_array
        global dts

        newTol = 0.35
        #  Convergence test
        for j = 1:4
        alg = ImplicitEulerExtrapolation(min_order = j,
            init_order = j, max_order=j,
            sequence = seq, threading = false)
        sim = test_convergence(dts,prob,alg)
        @test sim.𝒪est[:final] ≈ alg.n_init+1.1 atol=newTol #Superconvergence
        end
        # Regression test
        sol = solve(prob,ImplicitEulerExtrapolation(max_order = 9, min_order = 1,
            init_order = 9,sequence = seq,threading = false),reltol=1e-3)
        @test length(sol.u) < 15
    end
    end

    println("Testing ImplicitEulerBarycentricExtrapolation")
    @testset "Testing ImplicitEulerBarycentricExtrapolation" begin
        for prob in problem_array, seq in sequence_array
            global dts

            newTol = 0.35
            #  Convergence test
            for j = 1:4
            alg = ImplicitEulerBarycentricExtrapolation(min_order = j,
                init_order = j, max_order=j,
                sequence = seq, threading = false)
            sim = test_convergence(dts,prob,alg)
            @test sim.𝒪est[:final] ≈ alg.n_init + 0.5 atol=newTol #Superconvergence
            end
            # Regression test
            sol = solve(prob,ImplicitEulerBarycentricExtrapolation(max_order = 9, min_order = 1,
                init_order = 9,sequence = seq, threading = false),reltol=1e-3)
            @test length(sol.u) < 15
        end
    end

    println("Testing ImplicitDeuflhardExtrapolation")
    @testset "Testing ImplicitDeuflhardExtrapolation" begin
    for prob in problem_array, seq in sequence_array
        global dts

        # Convergence test
        for j = 1:6
        alg = ImplicitDeuflhardExtrapolation(min_order = j,
            init_order = j, max_order=j,
            sequence = seq,threading=false)
        sim = test_convergence(dts,prob,alg)
        @test sim.𝒪est[:final] ≈ 2*(alg.n_init+1) atol=testTol
        end

        # Regression test
        alg = ImplicitDeuflhardExtrapolation(max_order=9, min_order=1,
        init_order=9, sequence=seq,threading=false)
        sol = solve(prob, alg, reltol=1e-3)
        @test length(sol.u) < 10
    end
    end

    println("Testing ImplicitHairerWannerExtrapolation")
    @testset "Testing ImplicitHairerWannerExtrapolation" begin
    for prob in problem_array, seq in sequence_array
        global dts

        # Convergence test
        for j = 1:6
        alg = ImplicitHairerWannerExtrapolation(min_order = j,
            init_order = j, max_order=j,
            sequence = seq, threading=false)
        sim = test_convergence(dts,prob,alg)
        @test sim.𝒪est[:final] ≈ 2*(alg.n_init+1) - 1 atol=testTol
        end

        alg = ImplicitHairerWannerExtrapolation(max_order=9, min_order=1,
        init_order=9, sequence=seq,threading=false)
        sol = solve(prob, alg, reltol=1e-3)
        @test length(sol.u) < 10
    end
    end

    # Test ExtrapolationMidpointDeuflhard

    println("Testing ExtrapolationMidpointDeuflhard")
    @testset "Testing ExtrapolationMidpointDeuflhard" begin
    @testset "Testing sequential ExtrapolationMidpointDeuflhard" begin
        for prob in problem_array, seq in sequence_array
        global dts

        # Convergence test
        for j = 1:6
            alg = ExtrapolationMidpointDeuflhard(min_order = j,
            init_order = j, max_order=j,
            sequence = seq, threading=false)
            sim = test_convergence(dts,prob,alg)
            @test sim.𝒪est[:final] ≈ 2*(alg.n_init+1) atol=testTol
        end

        # Regression test
        alg = ExtrapolationMidpointDeuflhard(max_order=9, min_order=1,
            init_order=9, sequence=seq, threading=false)
        sol = solve(prob, alg, reltol=1e-3)
        @test length(sol.u) < 10
        end
    end
    end # ExtrapolationMidpointDeuflhard

    # Test ExtrapolationMidpointHairerWanner
    println("Testing ExtrapolationMidpointHairerWanner")
    @testset "Testing ExtrapolationMidpointHairerWanner" begin
    @testset "Testing sequential ExtrapolationMidpointHairerWanner" begin
        for prob in problem_array, seq in sequence_array
        global dts

        # Convergence test
        for j = 1:6
            alg = ExtrapolationMidpointHairerWanner(min_order = j,
            init_order = j, max_order=j,
            sequence = seq, threading=false)
            sim = test_convergence(dts,prob,alg)
            @test sim.𝒪est[:final] ≈ 2*(alg.n_init+1) atol=testTol
        end

        # Regression test
        alg = ExtrapolationMidpointHairerWanner(max_order=9, min_order=2,
            init_order=9, sequence=seq, threading=false)
        sol = solve(prob, alg, reltol=1e-3)
        @test length(sol.u) < 10
        end
    end

    end # ExtrapolationMidpointHairerWanner

    println("Regression Test Float32 and Float64 Fallbacks")
    @testset "Regression Test Float32 and Float64 Fallbacks" begin
    prob_ode_2Dlinear = ODEProblem(
                        ODEFunction(f_2dlinear,analytic=f_2dlinear_analytic),
                        Float64.(prob_ode_bigfloat2Dlinear.u0),(0.0,1.0),1.01)
    s1 = solve(prob_ode_bigfloat2Dlinear,ExtrapolationMidpointDeuflhard())
    s2 = solve(prob_ode_2Dlinear,ExtrapolationMidpointDeuflhard())
    @test all(all(s1[i] - s2[i] .< 5e-14) for i in 1:length(s1))

    prob_ode_2Dlinear = ODEProblem(
                        ODEFunction(f_2dlinear,analytic=f_2dlinear_analytic),
                        Float32.(prob_ode_bigfloat2Dlinear.u0),(0.0f0,1.0f0),1.01f0)
    s1 = solve(prob_ode_bigfloat2Dlinear,ExtrapolationMidpointDeuflhard())
    s2 = solve(prob_ode_2Dlinear,ExtrapolationMidpointDeuflhard())
    @test all(all(s1[i] - s2[i] .< 5e-6) for i in 1:length(s1))
    end

end # Extrapolation methods


text:
function N(p) #Generates set of neighbors. To be replaced by solution to issue #6. For now, this will just be a matrix of MxK with K copies of p
      y = Array[]
      for i in 1:3
          push!(y,p)
      end
      return y
  end

max_stars_repo_name: ScPo-CompEcon/CourseMatch

max_issues_repo_name: ScPo-CompEcon/CourseMatch

max_forks_repo_name: ScPo-CompEcon/CourseMatch


text:
using Distances

export GibbsRandMove, GibbsScanMove

struct GibbsRandMove <: McmcMove
    ν::Int
    counts::Vector{Int}
    function GibbsRandMove(;ν::Int=1)
        new(ν, [0,0])
    end
end

struct GibbsScanMove <: McmcMove
    ν::Int
    counts::Vector{Int}
    function GibbsScanMove(;ν::Int=1)
        new(ν, [0,0])
    end
end

Base.show(io::IO, x::GibbsRandMove) = print(io, "GibbsRandMove(ν=$(x.ν))")
Base.show(io::IO, x::GibbsScanMove) = print(io, "GibbsScanMove(ν=$(x.ν))")

function accept_reject_entry!(
    x_curr::Vector{Int},
    x_prop::Vector{Int},
    i::Int,
    move::Union{GibbsRandMove,GibbsScanMove},
    model::VecMultigraphSNF
    )

    x_mode = model.mode
    # Proposal generation
    @inbounds w_curr = x_curr[i]
    d = 1-2rand(0:1)
    w_tmp = w_curr + d * rand(1:move.ν)
    w_prop = w_tmp ≥ 0 ? w_tmp : -w_tmp
    @inbounds x_prop[i] = w_prop

    # Eval acceptance probability
    γ, d = (model.γ, model.d)
    log_α = γ * (d(x_curr, x_mode) - d(x_prop, x_mode))

    if log(rand()) < log_α
        @inbounds x_curr[i] = w_prop
        move.counts[1] += 1
    else
        @inbounds x_prop[i] = w_curr
    end
end

# Optimised version for the hamming distances between mutlisets
function accept_reject_entry!(
    x_curr::Vector{Int},
    x_prop::Vector{Int},
    i::Int,
    move::Union{GibbsRandMove,GibbsScanMove},
    model::VecMultigraphSNF{Cityblock},
    )
    x_mode = model.mode
    # Proposal generation
    @inbounds w_curr = x_curr[i]
    d = 1-2rand(0:1)
    w_tmp = w_curr + d * rand(1:move.ν)
    w_prop = w_tmp ≥ 0 ? w_tmp : -w_tmp
    @inbounds x_prop[i] = w_prop

    # Eval acceptance probability
    γ, d = (model.γ, model.d)
    @inbounds w_mode = x_mode[i]
    log_α = γ * (abs(w_curr - w_mode) - abs(w_prop - w_mode))

    if log(rand()) < log_α
        @inbounds x_curr[i] = w_prop
        move.counts[1] += 1
    else
        @inbounds x_prop[i] = w_curr
    end
end

function accept_reject!(
    x_curr::Vector{Int},
    x_prop::Vector{Int},
    move::GibbsRandMove,
    model::VecMultigraphSNF
    )
    move.counts[2] += 1
    M = length(x_curr)
    i = rand(1:M)
    accept_reject_entry!(
        x_curr, x_prop,
        i,
        move, model
    )
end

function accept_reject!(
    x_curr::Vector{Int},
    x_prop::Vector{Int},
    move::GibbsScanMove,
    model::VecMultigraphSNF
    )
    move.counts[2] += length(x_curr)
    for i in eachindex(x_curr)
        accept_reject_entry!(
            x_curr, x_prop,
            i,
            move, model
        )
    end
end

# A stripped-back proposal generation function (used in posterio samplers)
function prop_sample!(
    x_curr::Vector{Int},
    x_prop::Vector{Int},
    move::GibbsRandMove
    )

    i = rand(1:length(x_curr))
    @inbounds w_curr = x_curr[i]
    d = 1-2rand(0:1)
    w_tmp = w_curr + d * rand(1:move.ν)
    w_prop = w_tmp ≥ 0 ? w_tmp : -w_tmp
    @inbounds x_prop[i] = w_prop

    return 0.0
end

function prop_sample_entry!(
    x_curr::Vector{Int},
    x_prop::Vector{Int},
    move::GibbsScanMove,
    i::Int
    )

    @inbounds w_curr = x_curr[i]
    d = 1-2rand(0:1)
    w_tmp = w_curr + d * rand(1:move.ν)
    w_prop = w_tmp ≥ 0 ? w_tmp : -w_tmp
    @inbounds x_prop[i] = w_prop

    return 0.0
end


text:
function gmm_matrix(b::AbstractParticleBelief, sigma::Float64)
    N = n_particles(b)
    P = zeros(N, N)
    gmm_matrix!(P, b, sigma)
end

function gmm_matrix!(P::Matrix{Float64}, b::AbstractParticleBelief, sigma::Float64)
    p = particles(b)
    N = length(p)
    P .= 1/(sqrt(2*pi)*sigma) # factor for all entries
    for (i,xi) in enumerate(p), j=i+1:N
        diff = xi - p[j]
        P[i,j] = P[j,i] *= exp(-diff^2 / (2*sigma^2))
    end
    P
end

function gmm_matrix(b::AbstractParticleBelief, sigma::Float64, tol_sigmas::Float64)
    p = particles(b)
    N = length(p)
    P = fill(1/(sqrt(2*pi)*sigma), N, N) # factor for all entries
    for (i,xi) in enumerate(p), j=i+1:N
        diff = xi - p[j]
        if abs(diff) < tol_sigmas
            P[i,j] = P[j,i] *= exp(-diff^2 / (2*sigma^2))
        else
            P[i,j] = P[j,i] = 0.0
        end
    end
    P
end


text:
# Memory copy 3D to return T_peak
using CUDA

function copy3D!(T2, T, Ci)
    ix = (blockIdx().x-1) * blockDim().x + threadIdx().x
    iy = (blockIdx().y-1) * blockDim().y + threadIdx().y
    iz = (blockIdx().z-1) * blockDim().z + threadIdx().z
    T2[ix,iy,iz] = T[ix,iy,iz] + Ci[ix,iy,iz]
    return
end

function memcopy3D()
    # Numerics
    nx, ny, nz = 1024, 512, 512                            # Number of grid points in dimensions x, y and z
    threads = (32, 8, 1)
    blocks  = (nx÷threads[1], ny÷threads[2], nz÷threads[3])

    nt  = 100                                              # Number of time steps

    # Array initializations
    T   = CUDA.zeros(Float64, nx, ny, nz)
    T2  = CUDA.zeros(Float64, nx, ny, nz)
    Ci  = CUDA.zeros(Float64, nx, ny, nz)

    # Initial conditions
    Ci .= 0.5
    T  .= 1.7
    T2 .= T

    t_tic = 0.0
    # Time loop
    for it = 1:nt
        if (it == 11) t_tic=time() end  # Start measuring time.
        @cuda blocks=blocks threads=threads copy3D!(T2, T, Ci)
        synchronize()
        T, T2 = T2, T
    end
    t_toc=time()-t_tic

    # Performance
    @show A_eff = (2*1+1)*1/1e9*nx*ny*nz*sizeof(Float64)      # Effective main memory access per iteration [GB] (Lower bound of required memory access: T has to be read and written: 2 whole-array memaccess; Ci has to be read: : 1 whole-array memaccess)
    t_it  = t_toc/(nt-10)                               # Execution time per iteration [s]
    T_eff = A_eff/t_it                                  # Effective memory throughput [GB/s]
    println("time_s=$t_toc T_eff=$T_eff")
end

memcopy3D()

max_stars_repo_name: luraess/amdgpu-test

max_issues_repo_name: luraess/amdgpu-test

max_forks_repo_name: luraess/amdgpu-test


text:
module TypeTraits

using ..Prerequisites: atol,rtol,rawtype

export efficientoperations
export MemoryOrder,FOrder,COrder
export forder,corder,indtosub,subtoind

struct EfficientOperations end
"""
    efficientoperations

Indicate that the efficient operations, i.e. "=="/"isequal", "<"/"isless" or "replace", will be used.
"""
const efficientoperations=EfficientOperations()

"""
    ==(::EfficientOperations,o1,o2) -> Bool
    isequal(::EfficientOperations,o1,o2) -> Bool

Compare two objects and judge whether they are eqaul to each other.
"""
@generated function Base.:(==)(::EfficientOperations,o1,o2)
    fcount=o1|>fieldcount
    if fcount==o2|>fieldcount
        if fcount==0
            return :(true)
        else
            expr=:(getfield(o1,1)==getfield(o2,1))
            for i=2:fcount
                expr=Expr(:&&,expr,:(getfield(o1,$i)==getfield(o2,$i)))
            end
            return expr
        end
    else
        return :(false)
    end
end
@generated function Base.isequal(::EfficientOperations,o1,o2)
    fcount=o1|>fieldcount
    if fcount==o2|>fieldcount
        if fcount==0
            return :(true)
        else
            expr=:(isequal(getfield(o1,1),getfield(o2,1)))
            for i=2:fcount
                expr=Expr(:&&,expr,:(isequal(getfield(o1,$i),getfield(o2,$i))))
            end
            return expr
        end
    else
        return :(false)
    end
end

"""
    <(::EfficientOperations,o1,o2) -> Bool
    isless(::EfficientOperations,o1,o2) -> Bool

Compare two objects and judge whether the first is less than the second.
"""
@generated function Base.:<(::EfficientOperations,o1,o2)
    n1,n2=o1|>fieldcount,o2|>fieldcount
    N=min(n1,n2)
    expr=n1<n2 ? Expr(:if,:(getfield(o1,$N)==getfield(o2,$N)),true,false) : false
    expr=Expr(:if,:(getfield(o1,$N)<getfield(o2,$N)),true,expr)
    for i in range(N-1,stop=1,step=-1)
        expr=Expr(:if,:(getfield(o1,$i)>getfield(o2,$i)),false,expr)
        expr=Expr(:if,:(getfield(o1,$i)<getfield(o2,$i)),true,expr)
    end
    return expr
end
@generated function Base.isless(::EfficientOperations,o1,o2)
    n1,n2=o1|>fieldcount,o2|>fieldcount
    N=min(n1,n2)
    expr=n1<n2 ? Expr(:if,:(getfield(o1,$N)==getfield(o2,$N)),true,false) : false
    expr=Expr(:if,:(isless(getfield(o1,$N),getfield(o2,$N))),true,expr)
    for i in range(N-1,stop=1,step=-1)
        expr=Expr(:if,:(isless(getfield(o2,$i),getfield(o1,$i))),false,expr)
        expr=Expr(:if,:(isless(getfield(o1,$i),getfield(o2,$i))),true,expr)
    end
    return expr
end

"""
    isapprox(::EfficientOperations,::Val{Names},o1,o2;atol::Real=atol,rtol::Real=rtol) where Names -> Bool

Compare two objects and judge whether they are inexactly equivalent to each other.
"""
@generated function Base.isapprox(::EfficientOperations,::Val{Names},o1,o2;atol::Real=atol,rtol::Real=rtol) where Names
    fcount=o1|>fieldcount
    if fcount==o2|>fieldcount
        if fcount==0
            return :(true)
        else
            exprs=[(name=QuoteNode(name);:(isapprox(getfield(o1,$name),getfield(o2,$name);atol=atol,rtol=rtol))) for name in Names]
            for name in fieldnames(o1)
                name ∉ Names && (name=QuoteNode(name);push!(exprs,:(isequal(getfield(o1,$name),getfield(o2,$name)))))
            end
            return Expr(:&&,exprs...)
        end
    else
        return :(false)
    end
end

"""
    replace(::EfficientOperations,o;kwargs...) -> typeof(o)

Return a copy of the input object with some of the field values replaced by the keyword arguments.
"""
@generated function Base.replace(::EfficientOperations,o;kwargs...)
    exprs=[:(get(kwargs,$name,getfield(o,$name))) for name in QuoteNode.(o|>fieldnames)]
    return :(rawtype(typeof(o))($(exprs...)))
end

abstract type MemoryOrder end
struct FOrder <: MemoryOrder end
struct COrder <: MemoryOrder end
"""
    forder

Indicate that the convertion between Cartesian index and linear index is using the Fortran order.
"""
const forder=FOrder()
"""
    corder

Indicate that the convertion between Cartesian index and linear index is using the C/C++ order.
"""
const corder=COrder()

@generated tail(ts::NTuple{N}) where N=Expr(:tuple,[:(ts[$i]) for i=2:N]...)
@generated head(ts::NTuple{N}) where N=Expr(:tuple,[:(ts[$i]) for i=1:N-1]...)

"""
    indtosub(dims::Tuple,ind::Int,order::FOrder) -> Tuple
    indtosub(dims::Tuple,ind::Int,order::COrder) -> Tuple

Convert an linear index to Cartesian index. Fortran-order or C-order can be assigned.
"""
function indtosub(dims::Tuple,ind::Int,order::FOrder)
    length(dims)==0 && return ()
    ((ind-1)%dims[1]+1,indtosub(tail(dims),(ind-1)÷dims[1]+1,order)...)
end
function indtosub(dims::Tuple,ind::Int,order::COrder)
    length(dims)==0 && return ()
    (indtosub(head(dims),(ind-1)÷dims[end]+1,order)...,(ind-1)%dims[end]+1)
end

"""
    subtoind(dims::NTuple{N,Int},inds::NTuple{N,Int},order::FOrder) where N -> Int
    subtoind(dims::NTuple{N,Int},inds::NTuple{N,Int},order::COrder) where N -> Int

Convert an Cartesian index to linear index. Fortran-order or C-order can be assigned.
"""
function subtoind(dims::NTuple{N,Int},inds::NTuple{N,Int},order::FOrder) where N
    length(dims)==0 && return 1
    (subtoind(tail(dims),tail(inds),order)-1)*dims[1]+inds[1]
end
function subtoind(dims::NTuple{N,Int},inds::NTuple{N,Int},order::COrder) where N
    length(dims)==0 && return 1
    (subtoind(head(dims),head(inds),order)-1)*dims[end]+inds[end]
end

end # module


text:
function read_obj(filepath::String)
    if splitext(filepath)[2] == ".obj"
        reader = vtk.vtkOBJReader()
        reader[:SetFileName](filepath)
        reader[:Update]()

        polydata = reader[:GetOutput]()
        return extract_simple_data(polydata)
    else
        throw("The file input is not an obj file.")
    end
end

function read_ply(filepath::String)
    if splitext(filepath)[2] == ".ply"
        reader = vtk.vtkPLYReader()
        reader[:SetFileName](filepath)
        reader[:Update]()

        polydata = reader[:GetOutput]()
        return extract_simple_data(polydata)
    else
        throw("The file input is not a ply file.")
    end
end


text:
using Documenter, ContinuousWavelets
ENV["GKSwstype"] = "100"
makedocs(sitename="ContinuousWavelets.jl",
         format=Documenter.HTML(),
         authors="David Weber",
         clean=true,
         pages=[
             "basic usage" => "index.md",
             "Install" => "installation.md",
             "CWT" => [
                 "Available Wavelet Families" => "coreType.md",
                 "CWT Construction" => "CWTConstruction.md",
                 "Wavelet Frequency Spacing" => "spacing.md",
                 "Boundary Conditions" => "bound.md",
                 "Inversion" => "inverse.md"
             ],
         ])

deploydocs(
    repo="github.com/dsweber2/ContinuousWavelets.jl.git",
)


text:
using NCDatasets

gsr = NCDataset("./silicon_TB09o_GSR.nc")
ref_energies = Dict{String, Float64}(
    "Entropy"        => gsr["e_entropy"][:],
    "Ewald"          => gsr["e_ewald"][:],
    "PspCorrection"  => gsr["e_corepsp"][:],
    "Xc"             => gsr["e_xc"][:],
    "Kinetic"        => gsr["e_kinetic"][:],
    "Hartree"        => gsr["e_hartree"][:],
    "AtomicLocal"    => gsr["e_localpsp"][:],
    # "AtomicNonlocal" => gsr["e_nonlocalpsp"][:],
)
ref_evals_ABINIT = gsr["eigenvalues"][:,:,1]
ref_ABINIT = gsr["etotal"][:]
ref_fermi = gsr["fermie"][:]

@show ref_energies
@show ref_evals_ABINIT
@show ref_fermi
@show ref_ABINIT