text:
using ..Registrator: decodeb64

function is_pfile_parseable(c::AbstractString)
    @debug("Checking whether (Julia)Project.toml is non-empty and parseable")
    if length(c) != 0
        try
            TOML.parse(c)
            return true, nothing
        catch ex
            if isa(ex, CompositeException) && isa(ex.exceptions[1], TOML.ParserError)
                err = "Error parsing project file"
                @debug(err)
                return false, err
            else
                rethrow(ex)
            end
        end
    else
        err = "Project file is empty"
        @debug(err)
        return false, err
    end
end

function pfile_hasfields(p::Pkg.Types.Project)
    @debug("Checking whether (Julia)Project.toml contains name, uuid and version")
    try
        if p.name === nothing || p.uuid === nothing || p.version === nothing
            err = "Project file should contain name, uuid and version"
            @debug(err)
            return false, err
        elseif !isempty(p.version.prerelease)
            err = "Pre-release version not allowed"
            @debug(err)
            return false, err
        elseif p.version == v"0"
            err = "Package version must be greater than 0.0.0"
            @debug(err)
            return false, err
        end
    catch ex
        err = "Error reading (Julia)Project.toml: $(ex.msg)"
        @debug(err)
        return false, err
    end

    return true, nothing
end

function verify_projectfile_from_sha(reponame, commit_sha; auth=GitHub.AnonymousAuth(), subdir = "")
    project = nothing
    projectfile_found = false
    projectfile_valid = false
    err = nothing
    @debug("Getting gitcommit object for sha")
    gcom = gitcommit(reponame, GitCommit(Dict("sha"=>commit_sha)); auth=auth)
    @debug("Getting tree object for sha")
    recurse = subdir != ""
    t = tree(reponame, Tree(gcom.tree); auth=auth, params = Dict(:recursive => recurse))
    tree_sha = t.sha
    project_files = joinpath.(subdir, Base.project_names)

    for tr in t.tree, file in project_files
        if tr["path"] == subdir
            tree_sha = tr["sha"]
        elseif tr["path"] == file
            projectfile_found = true
            @debug("(Julia)Project file found")

            @debug("Getting projectfile blob")
            if isa(auth, GitHub.AnonymousAuth)
                a = get_user_auth()
            else
                a = auth
            end
            b = blob(reponame, Blob(tr["sha"]); auth=a)

            @debug("Decoding base64 projectfile contents")
            projectfile_contents = decodeb64(b.content)

            @debug("Checking project file validity")
            projectfile_valid, err = is_pfile_parseable(projectfile_contents)

            if projectfile_valid
                project = Pkg.Types.read_project(copy(IOBuffer(projectfile_contents)))
                projectfile_valid, err = pfile_hasfields(project)
            end
            break
        end
    end

    return project, tree_sha, projectfile_found, projectfile_valid, err
end


text:
using GPUCompiler

if !@isdefined(TestRuntime)
    include("../util.jl")
end

# create a SPIRV-based test compiler, and generate reflection methods for it

function spirv_job(@nospecialize(func), @nospecialize(types); kernel::Bool=false, kwargs...)
    source = FunctionSpec(func, Base.to_tuple_type(types), kernel)
    target = SPIRVCompilerTarget()
    params = TestCompilerParams()
    CompilerJob(target, source, params), kwargs
end

function spirv_code_typed(@nospecialize(func), @nospecialize(types); kwargs...)
    job, kwargs = spirv_job(func, types; kwargs...)
    GPUCompiler.code_typed(job; kwargs...)
end

function spirv_code_warntype(io::IO, @nospecialize(func), @nospecialize(types); kwargs...)
    job, kwargs = spirv_job(func, types; kwargs...)
    GPUCompiler.code_warntype(io, job; kwargs...)
end

function spirv_code_llvm(io::IO, @nospecialize(func), @nospecialize(types); kwargs...)
    job, kwargs = spirv_job(func, types; kwargs...)
    GPUCompiler.code_llvm(io, job; kwargs...)
end

function spirv_code_native(io::IO, @nospecialize(func), @nospecialize(types); kwargs...)
    job, kwargs = spirv_job(func, types; kwargs...)
    GPUCompiler.code_native(io, job; kwargs...)
end

# aliases without ::IO argument
for method in (:code_warntype, :code_llvm, :code_native)
    spirv_method = Symbol("spirv_$(method)")
    @eval begin
        $spirv_method(@nospecialize(func), @nospecialize(types); kwargs...) =
            $spirv_method(stdout, func, types; kwargs...)
    end
end

# simulates codegen for a kernel function: validates by default
function spirv_code_execution(@nospecialize(func), @nospecialize(types); kwargs...)
    job, kwargs = spirv_job(func, types; kernel=true, kwargs...)
    GPUCompiler.compile(:asm, job; kwargs...)
end

ext: jl


text:
using Statistics, SparseArrays
using MAT, DelimitedFiles, ForwardDiff
"""
``
eqcond(m::TwoAssetHANK)
```

Expresses the equilibrium conditions in canonical form using Γ0, Γ1, C, Ψ, and Π matrices.

### Outputs

* `Γ0` (`n_states` x `n_states`) holds coefficients of current time states.
* `Γ1` (`n_states` x `n_states`) holds coefficients of lagged states.
* `C`  (`n_states` x `1`) is a vector of constants
* `Ψ`  (`n_states` x `n_shocks_exogenous`) holds coefficients of iid shocks.
* `Π`  (`n_states` x `n_states_expectational`) holds coefficients of expectational states.
"""
function eqcond_lite(m::TwoAssetHANK)
    # Read in parameters
    aalpha = m[:aalpha].value::Float64
    ddelta = m[:ddelta].value::Float64
    ddeath = m[:ddeath].value::Float64
    rrho   = m[:rrho].value::Float64
    chi0   = m[:chi0].value::Float64
    chi1   = m[:chi1].value::Float64
    chi2   = m[:chi2].value::Float64
    a_lb   = m[:a_lb].value::Float64
    kappa  = m[:kappa].value::Float64
    pam    = m[:pam].value::Float64
    xxi    = m[:xxi].value::Float64
    ggamma = m[:ggamma].value::Float64
    tau_I  = m[:tau_I].value::Float64
    trans  = m[:trans].value::Float64
    n_SS   = m[:n_SS].value::Float64

    nnu_aggZ    = m[:nnu_aggZ].value::Float64
    ssigma_aggZ = m[:ssigma_aggZ].value::Float64

    # Set liquid rates
    r_b_borr_SS  = m[:r_b_borr_SS].value::Float64
    borrwedge_SS = m[:borrwedge_SS].value::Float64

    lambda::Matrix{Float64} = get_setting(m, :lambda)

    K_liquid                   = get_setting(m, :K_liquid)::Bool
    aggregate_variables        = get_setting(m, :aggregate_variables)==1#::Bool
    distributional_variables   = get_setting(m, :distributional_variables)==1#::Int64
    distributional_variables_1 = get_setting(m, :distributional_variables_1)==1#::Int64
    permanent                  = get_setting(m, :permanent)::Bool

    I      = get_setting(m, :I)::Int64
    J      = get_setting(m, :J)::Int64
    I_g    = get_setting(m, :I_g)::Int64
    J_g    = get_setting(m, :J_g)::Int64
    N      = get_setting(m, :N)::Int64

    a      = get_setting(m, :a)::Vector{Float64}
    b      = get_setting(m, :b)::Vector{Float64}
    a_g    = get_setting(m, :a_g)::Vector{Float64}
    b_g    = get_setting(m, :b_g)::Vector{Float64}

    agrid_new = get_setting(m, :agrid_new)::Int64
    bgrid_new = get_setting(m, :bgrid_new)::Int64
    ygrid_new = get_setting(m, :ygrid_new)::Int64

    amin = get_setting(m, :amin)
    bmin = get_setting(m, :bmin)
    amax = get_setting(m, :amax)
    bmax = get_setting(m, :bmax)

    y      = vec(get_setting(m, :y))::Vector{Complex{Float64}}
    y_dist = get_setting(m, :y_dist)::Vector{Complex{Float64}}
    y_mean = get_setting(m, :y_mean)::Complex{Float64}
    KL      = get_setting(m, :KL_0)::Float64
    r_b_fix = get_setting(m, :r_b_fix) ? 1 : 0::Int64

    #--- Taken from what used to be inside
    #a, a_g, a_g_0pos          = create_a_grid(agrid_new, J, J_g, amin, amax)
    #    _, _, _, b, b_g, b_g_0pos = create_b_grid(bgrid_new, I, I_g)
    #lambda, y, y_mean, y_dist, _ = create_y_grid(N, ygrid_new)

    n_v = get_setting(m, :n_v)::Int64
    n_g = get_setting(m, :n_g)::Int64
    n_p = get_setting(m, :n_p)::Int64
    n_Z = get_setting(m, :n_Z)::Int64
    nVars    = get_setting(m, :nVars)::Int64
    nEErrors = get_setting(m, :nEErrors)::Int64

    #vars_SS     = vec(m[:vars_SS].value)::Vector{Float64}
    vars_SS = vec(load("data/vars_SS.jld2","vars_SS"))
    V_SS   = reshape(vars_SS[1:n_v], I*J, N)
    g_SS   = vars_SS[n_v + 1 : n_v + n_g]
    K_SS   = vars_SS[n_v + n_g + 1]
    r_b_SS = vars_SS[n_v + n_g + 2]

    # Aggregate output and aggregate consumption
    aggY_SS = aggregate_variables ? vars_SS[n_v+n_g+3] : 0.0 # aggregate output
    aggC_SS = aggregate_variables ? vars_SS[n_v+n_g+4] : 0.0 # aggregate consumption

    # Consumption and earnings inequality
    C_Var_SS    = distributional_variables ? vars_SS[n_v+n_g+3] : 0.0
    earn_Var_SS = distributional_variables ? vars_SS[n_v+n_g+4] : 0.0

    # Consumption of wealthy and poor hand-to-mouth
    C_WHTM_SS = distributional_variables_1 ? vars_SS[n_v+n_g+3] : 0.0
    C_PHTM_SS = distributional_variables_1 ? vars_SS[n_v+n_g+4] : 0.0

    aggZ_SS = vars_SS[n_v+n_g+n_p+1] # Aggregate Z

    # Construct problem functions
    util, deposit, cost = construct_problem_functions(ggamma, chi0, chi1, chi2, a_lb)

    @inline function get_residuals_lite(vars::Vector{T}) where {T<:Real}
        # ------- Unpack variables -------
        V::Matrix{T} = reshape(vars[1:n_v] .+ V_SS[:,indx], I, J)  # Value function
        g::Vector{T} = vars[n_v + 1 : n_v + n_g] .+ g_SS           # Distribution
        K::T         = vars[n_v + n_g + 1] + K_SS                  # Aggregate capital
        r_b::T       = vars[n_v + n_g + 2] + r_b_SS

        # Aggregate output, aggregate consumption
        aggY::T     = aggregate_variables      ? vars[n_v+n_g+3] + aggY_SS     : zero(T)
        aggC::T     = aggregate_variables      ? vars[n_v+n_g+4] + aggC_SS     : zero(T)
        C_Var::T    = distributional_variables ? vars[n_v+n_g+3] + C_Var_SS    : zero(T)
        earn_Var::T = distributional_variables ? vars[n_v+n_g+4] + earn_Var_SS : zero(T)
        C_WHTM::T = distributional_variables_1 ? vars[n_v+n_g+3] + C_WHTM_SS   : zero(T)
        C_PHTM::T = distributional_variables_1 ? vars[n_v+n_g+4] + C_PHTM_SS   : zero(T)

        # Aggregate Z
        aggZ::T             = vars[n_v+n_g+n_p+1] + aggZ_SS
        V_Dot::Vector{T}    = vars[    nVarsi +       1 : nVarsi + n_v]
        g_Dot::Vector{T}    = vars[    nVarsi + n_v + 1 : nVarsi + n_v + n_g]
        aggZ_Dot::T         = vars[    nVarsi + n_v + n_g + n_p + 1]
        VEErrors::Vector{T} = vars[2 * nVarsi + 1 : 2 * nVarsi + n_v]
        aggZ_Shock::T       = vars[2 * nVarsi + nEErrors + 1]

        # Prices, liquid rates
        w   = (1 - aalpha) * (K ^ aalpha) * n_SS ^ (-aalpha) *
            (!permanent ? exp(aggZ) ^ (1-aalpha) : 1.)
        r_a = aalpha * (K ^ (aalpha - 1)) * ((!permanent ? exp(aggZ) : 1.0) *
            n_SS) ^ (1 - aalpha) - ddelta
        r_b_borr = r_b .+ borrwedge_SS

        # Set grid
        dab_g_tilde = backward_difference(a_g, b_g)

        dab_g_o = reshape(repeat(reshape(dab_g_tilde, I_g, J_g), N, 1), I_g, J_g, N)
        g_end = (1 - sum(g .* vec(dab_g_o)[1:end-1])) / dab_g_tilde[I_g*J_g]
        @show size(dab_g_o)
        @show size(g .* vec(dab_g_o)[1:end-1])
        @show size(vec(broadcast(*, dab_g_tilde, reshape(vcat(g, 0.0), I_g*J_g, N)))[1:end-1])

        g      = vcat(g, g_end)
        g_inds = reshape(g, I_g*J_g, N)
        gg     = g_inds[:,indx]

        # Other necessary objects
        y_shock      = y .* exp.(kappa * aggZ * (y .- y_mean) ./ std(y))
        y_shock_mean = dot(y_shock, y_dist)
        y_shock      = real(y_shock ./ y_shock_mean .* y_mean)

        println("Timing: solve_hjb()")
        @time c, s, d  = solve_hjb_lite(V, a_lb, ggamma, permanent, ddeath, pam, aggZ,
                                        xxi, tau_I, w, trans, r_b, r_b_borr, y_shock[indx],
                                        a, b, cost, util, deposit)

        @time interp_decision = interp(b_g, a_g, b, a)
        d_g = reshape(interp_decision * vec(d), I_g, J_g)
        s_g = reshape(interp_decision * vec(s), I_g, J_g)
        c_g = reshape(interp_decision * vec(c), I_g, J_g)

        @show size(interp_decision)
        # Derive transition matrices
        println("Timing: transition_deriva()")
        @time A, AT = transition_deriva_lite(permanent, ddeath, pam, xxi, w, a_lb, aggZ,
                                        d, d_g, s, s_g, r_a, a, a_g, b, b_g, y_shock[indx],
                                        lambda[indx,indx], cost, util, deposit)

        #----------------------------------------------------------------
        # KFE
        #----------------------------------------------------------------
        @time gIntermediate, cc = solve_kfe(dab_g_tilde, b, g_inds, gg, ddeath,
                                            lambda, AT, indx, V_SS)

        #----------------------------------------------------------------
        # Compute equilibrium conditions
        #----------------------------------------------------------------
        # HJB equation
        perm_mult   = !permanent ? rrho + ddeath : rrho + ddeath - (1 - ggamma) * aggZ
        hjbResidual = vec(util.(c)) + A * vec(V) + cc + V_Dot + VEErrors - perm_mult * vec(V)

        # KFE
        a_gg = repeat(a_g, inner=I_g)
        b_gg = repeat(b_g, outer=J_g)

        @show length(g_Dot), length(gIntermediate)
        gResidual  = g_Dot[I_g*J_g*(indx-1)+1 : I_g*J_g*indx] - gIntermediate
        K_Residual = (K_liquid ? sum((a_gg .+ b_gg) .* gg .* dab_g_tilde) :
                                 sum( a_gg          .* gg .* dab_g_tilde)) - K
        r_b_Residual = 0.0
        if r_b_fix      == 1
            r_b_out      = r_b_SS
            r_b_Residual = r_b_out - r_b
        elseif r_b_phi  == 1
            r_b_out      = sum(b_gg .* gg .* dab_g_tilde)
            r_b_Residual = r_b_out - B_SS * exp(1 / pphi * (r_b - r_b_SS))
        elseif B_fix    == 1
            # Find death-corrected savings
            b_save       = dot(gIntermediate, dab_g_tilde .* b_gg)
            r_b_out      = 0.0
            r_b_Residual = r_b_out - b_save
        elseif K_liquid
            r_b_out      = r_a_out - illiquid_wedge
            r_b_Residual = r_b_out - r_b
        end

        Y_Residual        = zero(T)
        C_Residual        = zero(T)

        C_Var_Residual    = Array{T}(undef, 0)
        earn_Var_Residual = Array{T}(undef, 0)

        C_WHTM_Residual   = Array{T}(undef, 0)
        C_PHTM_Residual   = Array{T}(undef, 0)

        if aggregate_variables

            aggY_out = (K ^ aalpha) * (!permanent ? exp(aggZ) * n_SS : n_SS ) ^ (1 - aalpha)
            aggC_out = sum(vec(c_g) .* gg .* dab_g_tilde)
            Y_Residual = aggY_out - aggY
            C_Residual = aggC_out - aggC

        elseif distributional_variables
            r_b_g_vec = r_b .* (b_g .>= 0) + r_b_borr .* (b_g .< 0)
            C_Var_out = sum(log(vec(c_g)).^2 .* gg .* dab_g_tilde) -
                sum(log(vec(c_g)) .* gg .* dab_g_tilde) ^ 2
            earn = log.((1-tau_I) * w * repeat(repeat(vec(y), inner=I_g), inner=J_g) .+
                         b_gg .* (repeat(repeat(r_b_g_vec, outer=J_g), outer=N) .+ ddeath*pam) .+
                         trans .+ a_gg .* (r_a + ddeath*pam))
            earn_Var_out = sum(vec(earn).^2 .* gg .* dab_g_tilde) -
                sum(vec(earn) .* gg .* dab_g_tilde) ^ 2

            C_Var_Residual    = C_Var_out - C_Var
            earn_Var_Residual = earn_Var_out - earn_Var

        elseif distributional_variables_1

            WHTM_indicator      = zeros(I_g,J_g,N)
            WHTM_indicator[b_g_0pos:b_g_0pos+1,a_g_0pos+2:end,:] .= 1.
            WHTM_out            = sum(vec(WHTM_indicator) .* gg .* dab_g_tilde)
            C_WHTM_out          = sum(vec(WHTM_indicator) .* vec(c_g) .* gg .* dab_g_tilde)

            PHTM_indicator      = zeros(I_g,J_g,N)
            PHTM_indicator[b_g_0pos:b_g_0pos+1,a_g_0pos:a_g_0pos+2:end,:] .= 1.
            PHTM_out            = sum(vec(PHTM_indicator) .* gg .* dab_g_tilde)
            C_PHTM_out          = sum(vec(PHTM_indicator) .* vec(c_g) .* gg .* dab_g_tilde)

            C_WHTM_Residual = C_WHTM_out - C_WHTM
            C_PHTM_Residual = C_PHTM_out - C_PHTM
        end

        # Law of motion for aggregate tfp shock
        aggZ_Residual = aggZ_Dot - (-nnu_aggZ * aggZ + ssigma_aggZ * aggZ_Shock)

        # Return equilibrium conditions
        #if aggregate_variables == 1
            return [hjbResidual; gResidual; K_Residual; r_b_Residual; Y_Residual;
                    C_Residual; aggZ_Residual]
        #elseif distributional_variables == 1
        #    return [hjbResidual; gResidual; K_Residual; r_b_Residual; C_Var_Residual;
        #            earn_Var_Residual; aggZ_Residual]
        #elseif distributional_variables_1 == 1
        #    return [hjbResidual; gResidual; K_Residual; r_b_Residual; C_WHTM_Residual;
        #            C_PHTM_Residual; aggZ_Residual]
        #end
        #return [hjbResidual; gResidual; K_Residual; r_b_Residual; aggZ_Residual]
    end

    indx = 1
    n_v = Int(n_v/N)
    #n_g = Int((n_g + 1)/N)
    nVarsi = n_v + n_g + n_p + 1

    out = get_residuals_lite(zeros(Float64, 2 * nVarsi + nEErrors + 1))
    @time out = get_residuals_lite(zeros(Float64, 2 * nVarsi + nEErrors + 1))

    #JLD2.jldopen("reference/eqcond_after_.jld2", true, true, true, IOStream) do file
    #    file["residuals"] = out
    #end

    my_out = vec(DelimitedFiles.readdlm("reference/my_residuals.csv", ','))

    @assert maximum(abs.(my_out[1:2000] - vec(out)[1:2000])) < 1e-7
    @show maximum(abs.(my_out[1:2000] - vec(out)[1:2000]))
    @assert isapprox(my_out[60001:62250], vec(out)[2001:4250], rtol=1e-4)

    x = zeros(Float64, 2 * nVarsi + nEErrors + 1)
    @time derivs = ForwardDiff.sparse_jacobian(get_residuals_lite, x)

    nstates = nVars    # n_states(m)
    n_s_exp = nEErrors # n_shocks_expectational(m)
    n_s_exo = n_Z      # n_shocks_exogenous(m)

    # vars = zeros(Float64, 2 * nstates + n_s_exp + n_s_exo)

    Γ1 = spzeros(Float64, nstates, nstates)#-derivs[:, 1:nstates]
    Γ0 = spzeros(Float64, nstates, nstates)# derivs[:,   nstates +  1:2 * nstates]
    Π  = spzeros(Float64, nstates, n_s_exp)#-derivs[:, 2*nstates +  1:2 * nstates + n_s_exp]
    Ψ  = spzeros(Float64, nstates, n_s_exo)#-derivs[:, 2*nstates + n_s_exp + 1:2 * nstates + n_s_exp + n_s_exo]
    C  =  spzeros(Float64, nstates)

    for i=1:1
        indx = i
        x = zeros(Float64, 2 * nVarsi + nEErrors + 1)
        @time derivs = ForwardDiff.sparse_jacobian(get_residuals_lite, x)
        #Γ1[] = -
        #Γ0[] =
        #Π[]  =
        #Ψ[]  =
    end

    if typeof(Ψ) == Vector{Float64}
        Ψ = reshape(Ψ, length(Ψ), n_s_exo)
    end

    test_out = load("data/eqcond_output_matlab.jld2")
    @assert test_out["g0"] == Γ0
    @assert test_out["g1"] ≈ Γ1
    @assert test_out["psi"] == Ψ
    @assert test_out["pi"] == Π

    return Γ0, Γ1, Ψ, Π, C
end


text:
module CLI
using Formatting
include("action.jl")

function reallocate_actions(
    proposed_ipfs::Vector{T},
    existing_ipfs::Vector{T},
    proposed_allocations::Dict{T,<:Real},
    existing_allocations::Dict{T,T},
) where {T<:AbstractString}
    ipfses = reallocate_ipfs(existing_ipfs, proposed_ipfs)
    actions = map(
        ipfs ->
            "\e[0mgraph indexer rules stop $(ipfs)\n\e[1m\e[38;2;255;0;0;249mCheck before submitting: \e[0mgraph indexer rules set $(ipfs) decisionBasis always allocationAmount $(format(proposed_allocations[ipfs]))",
        ipfses,
    )
    return actions, ipfses
end

function allocate_actions(
    proposed_ipfs::Vector{T},
    reallocate_ipfs::Vector{T},
    proposed_allocations::Dict{T,<:Real},
) where {T<:AbstractString}
    ipfses = open_ipfs(proposed_ipfs, reallocate_ipfs)
    actions = map(
        ipfs ->
            "\e[0mgraph indexer rules set $(ipfs) decisionBasis always allocationAmount $(format(proposed_allocations[ipfs]))",
        ipfses,
    )
    return actions, ipfses
end

function unallocate_actions(
    existing_ipfs::Vector{T}, reallocate_ipfs::Vector{T}, frozenlist::Vector{T}
) where {T<:AbstractString}
    ipfses = setdiff(setdiff(existing_ipfs, reallocate_ipfs), frozenlist)
    actions = map(ipfs -> "\e[0mgraph indexer rules stop $(ipfs)", ipfses)
    return actions, ipfses
end
end


text:
using Documenter, Secp256k1

makedocs(sitename="Secp256k1",
         doctest=true)


text:
@testset "gpar" begin
    N = 10
    D_in = 2
    D_out = 3

    rng = MersenneTwister(123456)
    f = GPAR([GP(SEKernel()) for _ in 1:D_out])
    x = ColVecs(randn(D_in, N))
    Σs = [rand(rng) + 0.1 for _ in 1:D_out]

    # Check that the output of rand is of the correct dimensionality.
    y = rand(rng, f(x, Σs))
    @test length(y) == N
    @test size(y.X) == (D_out, N)

    # Check that logpdf returns something, and said this is real number.
    @test logpdf(f(x, Σs), y) isa Real

    # Check that the posterior is another GPAR.
    f_post = posterior(f(x, Σs), y)
    @test f_post isa GPAR

    # Check that quantities can be computed under the posterior.
    N_pred = 5
    x_pred = ColVecs(randn(rng, D_in, N_pred))
    y_post = rand(rng, f_post(x_pred, Σs))
    @test size(rand(rng, f_post(x_pred, Σs))) == size(y_post)
    @test logpdf(f_post(x_pred, Σs), y_post) isa Real
end


text:
@testset "lmm" begin
    rng = MersenneTwister(123)
    FDM = FiniteDifferences.central_fdm(5, 1)
    N = 6
    in_dim = 3
    out_dim = 3
    x1IO = KernelFunctions.MOInputIsotopicByOutputs(
        [rand(rng, in_dim) for _ in 1:N], out_dim
    )
    x2IO = KernelFunctions.MOInputIsotopicByOutputs(
        [rand(rng, in_dim) for _ in 1:N], out_dim
    )
    x3IO = KernelFunctions.MOInputIsotopicByOutputs(
        [rand(rng, in_dim) for _ in 1:div(N, 2)], out_dim
    )

    latentkernels = [Matern32Kernel(), SqExponentialKernel(), FBMKernel(), Matern32Kernel()]
    H = rand(length(latentkernels), out_dim)
    k = LinearMixingModelKernel(latentkernels, H)

    badH = rand(length(latentkernels) - 1, out_dim)
    @test_throws AssertionError LinearMixingModelKernel(latentkernels, badH)

    @test k isa LinearMixingModelKernel
    @test k isa MOKernel
    @test k isa Kernel
    @test k(x1IO[1], x2IO[1]) isa Real

    @test string(k) == "Linear Mixing Model Multi-Output Kernel"
    @test repr("text/plain", k) == (
        "Linear Mixing Model Multi-Output Kernel. Kernels:\n" *
        "\tMatern 3/2 Kernel (metric = Euclidean(0.0))\n" *
        "\tSquared Exponential Kernel (metric = Euclidean(0.0))\n" *
        "\tFractional Brownian Motion Kernel (h = 0.5)\n" *
        "\tMatern 3/2 Kernel (metric = Euclidean(0.0))"
    )

    TestUtils.test_interface(k, x1IO, x2IO, x3IO)

    x1IF = KernelFunctions.MOInputIsotopicByFeatures(x1IO.x, out_dim)
    x2IF = KernelFunctions.MOInputIsotopicByFeatures(x2IO.x, out_dim)
    x3IF = KernelFunctions.MOInputIsotopicByFeatures(x3IO.x, out_dim)

    TestUtils.test_interface(k, x1IF, x2IF, x3IF)

    k = LinearMixingModelKernel(SEKernel(), H)

    @test k isa LinearMixingModelKernel
    @test k isa MOKernel
    @test k isa Kernel
    @test length(k.K) == 4
    for kernel in k.K
        @test isa(kernel, SEKernel)
    end

    @test string(k) == "Linear Mixing Model Multi-Output Kernel"
    @test repr("text/plain", k) == (
        "Linear Mixing Model Multi-Output Kernel. Kernels:\n" *
        "\tSquared Exponential Kernel (metric = Euclidean(0.0))\n" *
        "\tSquared Exponential Kernel (metric = Euclidean(0.0))\n" *
        "\tSquared Exponential Kernel (metric = Euclidean(0.0))\n" *
        "\tSquared Exponential Kernel (metric = Euclidean(0.0))"
    )

    test_ADs(k)
end


text:
function compact_ata(;
    settings::InputSettings = InputSettings(),
    bank::DataFrames.DataFrame = DataFrames.DataFrame(),
    settings_file = "",
    bank_file = "",
    bank_delim = ";",
    add_constraints = true,
    constraints::DataFrames.DataFrame = DataFrames.DataFrame(),
    constraints_file = "",
    constraints_delim = ";",
    add_overlap = true,
    overlap::Matrix{Float64} = Matrix{Float64}(undef, 0, 0),
    overlap_file = "",
    overlap_delim = ";",
    add_obj_fun = true,
    solver = "siman", #"jump",
    return_results = true,
    results_folder = "results",
    return_plots = true,
    plots_folder = "plots",
    sim_pool = DataFrames.DataFrame(),
    psychometrics = false,
    kwargs...,
)
    ata_model = start_ata(;
        settings = settings,
        bank = bank,
        settings_file = settings_file,
        bank_file = bank_file,
        bank_delim = bank_delim,
    )

    if add_constraints
        add_constraints!(
            ata_model;
            constraints = constraints,
            constraints_file = constraints_file,
            constraints_delim = constraints_delim,
        )

    end
    if add_overlap
        add_overlap!(
            ata_model;
            overlap = overlap,
            overlap_file = overlap_file,
            overlap_delim = overlap_delim,
        )

    end
    if add_obj_fun
        add_obj_fun!(ata_model; psychometrics = psychometrics)

    end
    # 9. assemble
    assemble!(ata_model; solver = solver, kwargs...)
    if return_results
        print_results(ata_model; results_folder = results_folder, sim_pool = sim_pool)
    end
    if return_plots
        plot_results(ata_model; plots_folder = plots_folder)
    end
    print_infos(ata_model)
    return ata_model::AbstractModel
end


text:
include("../../neat_test.jl")

const DEFAULT_DATASET_PATH = "../../../../src/datasets"

println("--> Dataset 4")
config = "test_dataset4_config.txt"
#config = ARGS[1]
#generations = parse(Int, ARGS[2])
generations = 5000
st4 = run_NEAT_epochs(config, generations, false)
print_best_solution_perfomance(st4)
#print("* Calculating information plane points")
#calculate_information_plane_points(st4)
println()
print(st4.population.best_fitness[end])

plot_neat_fitness(st4, imagename="fitness-dt4.png")
plot_neat_IP(st4, imagename="IP-dt4.png", nbins=10)
plot_output_distribution(st4, imagename="distribution-dt4.png")
plot_complexity(st4, imagename="solution-complexity-dt4.png")
plot_total_accuracy(st4, imagename="total-accuracy-dt4.png")
plot_output_distribution_test(st4, imagename="output-test-distribution-best-NEAT-dt4.png")

#=
plot_neat_fitness(st4, imagename="fitness-dt4.png")
plot_neat_IP(st4, imagename="IP-dt4.png", nbins=10)
plot_neat_species(st4, imagename="species-dt4.png")
plot_output_distribution(st4, imagename="distribution-dt4.png")
plot_complexity(st4, imagename="solution-complexity-dt4.png")
plot_accuracy(st4, imagename="accuracy-dt4.png")
plot_total_accuracy(st4, imagename="total-accuracy-dt4.png")
plot_output_distribution_test(st4, imagename="output-test-distribution-best-NEAT-dt4.png")
plot_confusion_matrix(st4; imagename="confusion-matrix-testset-dt4.png", set=:test)
=#

# nohup julia neat_test_dt4.jl &> nohup-NEAT-dt4-relatorio.out &

max_stars_repo_name: Lucas-Muniz/MNEAT

max_issues_repo_name: Lucas-Muniz/MNEAT

max_forks_repo_name: Lucas-Muniz/MNEAT


text:
using Cascadia
using Gumbo
using Requests
using PlotlyJS
using DataStructures: OrderedDict

include("scraper/heise.jl")
include("scraper/golem.jl")
include("scraper/spiegel.jl")

ignore_year = 2019

data = Dict()
#               scrape_...(maxpages, ignore_year)
data["heise"] = scrape_heise(43, ignore_year)
data["golem"] = scrape_golem(6, ignore_year)
data["spiegel"] = scrape_spiegel(20, ignore_year)

# plotting the results
trace_heise = scatter(;x=collect(keys(data["heise"])), y=collect(values(data["heise"])), name= "heise.de")
trace_golem = scatter(;x=collect(keys(data["golem"])), y=collect(values(data["golem"])), name= "golem.de")
trace_spiegel = scatter(;x=collect(keys(data["spiegel"])), y=collect(values(data["spiegel"])), name="spiegel-online.de")

records = [trace_heise, trace_golem, trace_spiegel]

layout = Layout(;
    xaxis = Dict(
        :title => "Jahr der Veröffentlichung",
        :range => [2000, ignore_year-1],
        :dtick => 2
    ),
    margin = Dict(
        :l => 50,
        :r => 10,
        :b => 60,
        :t => 50
    ),
    yaxis = Dict(
        :title => "Anzahl der Artikel"
    ),
    title = "Artikel über künstliche Intelligenz pro Jahr (Stand: $(Dates.today()))"
)

println("plotting...")
plot(records, layout)


text:
### A Pluto.jl notebook ###
# v0.19.0

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ aa854f02-46bf-4a2b-b1e9-5fb52b99925c
begin
    import Pkg
    Pkg.activate(Base.current_project())
    Pkg.instantiate()
end

# ╔═╡ 8b4426b8-588b-49ad-83c5-a79ed698d704
using PlutoUI

# ╔═╡ d47c1cf4-4844-4c66-8a21-b62ad7740dcf
using JLD2

# ╔═╡ 3c07c33b-7f91-4e2d-a27a-09a8924f328c
using DataFrames

# ╔═╡ b59dd2c0-4c3e-4c3d-a4d4-f8da9dcbe286
using CSV

# ╔═╡ 76e2e659-c708-4a08-8263-def6c3dfa930
using Dates

# ╔═╡ 41b0f95e-d806-4167-8199-7392717782f9
using OffsetArrays

# ╔═╡ 95f2bd0b-7d8b-43b9-b2f9-cfa6628ec485
using LinearAlgebra

# ╔═╡ 356cc4c2-e565-4847-9286-cfc1f835654d
using StatsBase

# ╔═╡ 730340bb-8176-4282-8559-8d60b1d28d57
using Statistics

# ╔═╡ f07ab63f-dfb4-493c-a40c-6b3363a858c3
using Plots

# ╔═╡ 342cbcff-5a05-4199-b5b5-6cc0fa70c202
using TemplateMatching

# ╔═╡ 7ddcab3c-a100-44be-baa3-af28cc08b2d7
using Interpolations

# ╔═╡ 341b9066-dcb3-4d99-b985-82e24c11f28e
using Optim

# ╔═╡ c2f22858-6d3c-4f92-a05b-659deae566f6
md"## Reading Continuous Data"

# ╔═╡ a86b6c8a-9b90-4388-8208-d257ec6951d2
md"""Path of the JLD2 file containing continuous data: $(@bind datapath TextField(default="../data/2021-01-12_20-25-30/data.jld2"))"""

# ╔═╡ f3f17c65-dab6-46da-83a9-87aad1181524
data = load(datapath, "data")

# ╔═╡ 8271f7e6-d8ac-472a-9538-75d0afacf42e
samplefreq = 2 # MHz

# ╔═╡ 57b35604-727b-4eac-bcbd-ea302e4c79a2
md"## Reading events catalogue"

# ╔═╡ b62d7318-3a6d-4442-b6cc-3166b87dff8a
md"""Path of the CSV catalogue $(@bind cataloguepath TextField(default="../data/2021-01-12_20-25-30/cat.csv"))"""

# ╔═╡ f3139f19-e386-4bb7-8aca-fa54062d4093
begin
	import Base.+, Base.show

	+(x::Dates.UTInstant{Microsecond}, d::Microsecond) = Dates.UTInstant(x.periods + d)

	function show(io::IO, ::MIME"text/plain", x::Dates.UTInstant{Microsecond})
		date = DateTime(Dates.UTM(round(Int, 1e-3 * x.periods.value)))
		show(io, MIME"text/plain"(), date)
	end
end

# ╔═╡ 9de113e1-36ce-4251-9fa5-5baa27941a85
function DateTimeMicrosecond(y, m, d, h, mi, s, us)
	rata = us + 1_000_000 * (s + 60mi + 3600h + 86400 * Dates.totaldays(y, m, d))
	Dates.UTInstant{Microsecond}(Microsecond(rata))
end

# ╔═╡ 96cf6e3b-1c2b-49f6-979d-0fce28a55a06
starttime = DateTimeMicrosecond(2021, 01, 12, 20, 25, 30, 0)

# ╔═╡ 5d5f3ca7-260d-4dbb-bde1-59f2bd58c0c9
begin
	df = CSV.read(cataloguepath, DataFrame)
	sec = floor.(Int, df.Second)
	usec = round.(Int, 1e6 * (df.Second - sec))
	templates = DataFrame()
	templates.datetime = DateTimeMicrosecond.(df.Year, df.Month, df.Day, df.Hour, df.Minute, sec, usec)
	templates.sample = map(x -> round(Int, x.value * samplefreq), templates.datetime .- starttime)
	templates.north = df.North
	templates.east = df.East
	templates.up = df.Up
	templates
end

# ╔═╡ 6ac6cce3-09ae-421a-b502-b5c8807bf555
md"## Reading Sensors Positions"

# ╔═╡ be090485-7917-464a-8e8c-aec4a2c009d3
begin
	sensors_positions = CSV.read("../data/2021-01-12_20-25-30/passive-xyz.csv", DataFrame, header=[:north, :east, :up])
	sensors_positions .*= 100
	sensors_positions.sensor = axes(sensors_positions, 1) .- 1
	filter!(r -> r.sensor in keys(data), sensors_positions)
	sensors_positions
end

# ╔═╡ 6ab4a8da-4613-4004-807f-c36876285ebd
md" ## Preparing a single template"

# ╔═╡ fe42e7e6-2c56-4f11-aae6-ec06ad7a219f
md"""template number: $(@bind template_num Select(1:nrow(templates), default=116))"""

# ╔═╡ 298010b3-c34c-40ba-9e1b-b79f7378f1d4
begin
	displacement = DataFrame()
	for s in [:north, :east, :up]
		displacement[!, s] = sensors_positions[!, s] .- templates[template_num, s]
	end
	displacement
end

# ╔═╡ a1fc5e84-48df-45ea-b6a1-0073ed9f8340
distances = hypot.(displacement[!, :north], displacement[!, :east], displacement[!, :up]) # in cm

# ╔═╡ c7f0582a-8997-4641-b607-0b95a7f73a6c
v_p = 0.67345 # in cm/us

# ╔═╡ 074983c0-2747-4b09-a0c0-4d4adba24c0b
offsets = @. round(Int, samplefreq * (distances / v_p))

# ╔═╡ cd5f912a-fafb-4947-bd70-318151c2c49c
shifts = templates[template_num, :sample] .+ offsets

# ╔═╡ 033473d3-6958-4e7b-8c95-a9e15d882118
md"""
Template window pre (in samples): $(@bind t_pre Slider(0:10:1000, default=100, show_value=true))

Template window post (in samples): $(@bind t_post Slider(0:10:1000, default=500, show_value=true))
"""

# ╔═╡ be80dd80-d39d-4338-bb76-c644ea9d970c
cuts = [shift - t_pre: shift + t_post for shift in shifts]

# ╔═╡ aa23d53f-8333-4c94-8932-fac99ce9196d
template_data = [view(series, cut) for (series, cut) in zip(values(data), cuts)]

# ╔═╡ 4ef908f4-d8dc-4159-97ad-fc8d4621a652
let
	template_aligned = OffsetVector.(template_data, shifts)
	xlim = minimum(firstindex(series, 1) for series in template_aligned), maximum(lastindex(series, 1) for series in template_aligned)
	plots = []
	for (nch, series) = zip(keys(data), template_aligned)
		ylim = minimum(series), maximum(series)
		push!(plots,
			  plot(axes(series, 1), series;
				   xlim,
				   ylim,
			       title="Channel $nch",
				   titlealign=:left,
				   showaxis=false,
				   ticks=false,
				   label=nothing))
	end
	numchannels = length(template_data)
	plot(plots..., layout=(numchannels, 1), size=(800, 100numchannels), dpi=300)
end

# ╔═╡ 368f7fc1-0afc-439f-905a-76af04e8ca8d
md" ## Computing the cross-correlation"

# ╔═╡ 0d122aef-1767-4389-b043-97d07ecef1fe
md"""
tolerance $(@bind tolerance Slider(0:10, default=5, show_value=true))
"""

# ╔═╡ a7a589fc-4dcc-4137-a423-b8d2e7d41a0e
signal = correlatetemplate(collect(values(data)), template_data, offsets, tolerance)

# ╔═╡ 91924a12-6ee4-4d1b-8c49-9b54142cf733
let
	template_ref = templates[template_num, :sample] - t_pre
	window = t_pre + t_post - 1
	a = max(firstindex(signal, 1), template_ref - div(window, 2))
	b = min(lastindex(signal, 1), template_ref + div(window, 2))
	plot(a:b, signal[a:b], label=nothing, dpi=300)
end

# ╔═╡ abb8b867-39c2-40f2-9d7a-f40bacf8f550
md" ## Finding peaks"

# ╔═╡ 7d1f1903-8342-4715-8e88-3e3d063ea5db
md"""
correlation threshold $(@bind threshold Slider(0.0:0.05:1, default=0.4, show_value=true))

distance (in unit of template length) $(@bind reldistance Slider(0:0.5:5, default=2, show_value=true))
"""

# ╔═╡ 666c8d5f-c1c5-4cb4-9364-bf997b18a683
peaks, heights = findpeaks(signal, threshold, reldistance * (t_pre + t_post - 1))

# ╔═╡ 96deaa3d-fc40-466f-bf1b-d8d3dd58dcd4
md"""
Peak number $(@bind peak_num Slider(axes(peaks, 1), default=1, show_value=true)) / $(length(peaks))
"""

# ╔═╡ cc0c3e0a-4208-4339-bea1-75396b2778bd
let
	plt = plot(ylims=(0.0, 1.0))
	scale = round(Int, reldistance * (t_post + t_pre - 1))
	a = max(firstindex(signal, 1), peaks[peak_num] - scale)
	b = min(lastindex(signal, 1), peaks[peak_num] + scale)
	plot!(plt, a:b, signal[a:b], label=nothing)
	plot!(plt, peaks[a .< peaks .< b], heights[a.< peaks .< b],
		seriestype=:scatter, label=nothing, dpi=300)
	plt
end

# ╔═╡ b6848d59-e97e-4662-a2ee-78612694d6d8
md"## Relocalization"

# ╔═╡ 89b036db-5c47-4a17-8466-9062e8d5d452
toas = estimatetoa.(collect(values(data)), template_data, peaks[peak_num] .+ offsets, tolerance)

# ╔═╡ 6d6d27c9-3b1a-4646-8420-61a775f7842f
let
	peak = peaks[peak_num]
	peak_start = peak - div(t_pre, 2)
	peak_stop = peak + (t_pre + t_post - 1) + div(t_post, 2)
	data_start = minimum(series -> firstindex(series, 1), values(data))
	data_stop = maximum(series -> lastindex(series, 1), values(data))
	xlim = max(data_start, peak_start), min(data_stop, peak_stop)

	plots = []
	for (n, channel_num) = enumerate(keys(data))
		series_start = 	max(firstindex(data[channel_num], 1), peak_start)
		series_stop = min(lastindex(data[channel_num], 1), peak_stop)
		data_series = OffsetVector(view(data[channel_num], series_start:series_stop), series_start:series_stop)

		toa_f, toa_i = modf(toas[n][1])
		template_interpolated = subsampleshift(template_data[n], toa_f)
		template_aligned = OffsetVector(template_interpolated, Int(toa_i))

		data_inf, data_sup = extrema(data_series)
		data_center = mean(data_series)
		template_inf, template_sup = extrema(template_interpolated)
		template_center = mean(template_aligned)
		β = (data_sup - data_inf) / (template_sup - template_inf)
		push!(plots,
			  plot([axes(data_series, 1), axes(template_aligned, 1)],
				   [data_series, data_center .+ β .* (template_aligned .- template_center)];
				   xlim,
				   ylim=(data_inf, data_sup),
				   title="Channel $channel_num",
				   titlealign=:left,
				   showaxis=false,
				   ticks=false,
				   label=nothing))
	end
	numchannels = length(template_data)
	plot(plots..., layout=(numchannels, 1), size=(800, 100numchannels), dpi=300)
end

# ╔═╡ 2911437d-4775-4a1b-819e-adf6e91acd61
measurements = let correlation_threshold = 0.5
	measurements = copy(sensors_positions)
	measurements.toas = [(sample + t_pre) / samplefreq for (sample, _) in toas]
	measurements.cc = [cc for (_, cc) in toas]
	filter(row -> row.cc > correlation_threshold, measurements)
end

# ╔═╡ 2e6b7b35-9512-406d-9a97-f1f94bea8732
if nrow(measurements) > 4
	x0 = Vector(templates[template_num, [:north, :east, :up]])
	t0 = (peaks[peak_num] + t_pre) / samplefreq
	sensors_readings = eachrow(Matrix(measurements[!, [:north, :east, :up, :toas]]))
	candidate = locate(sensors_readings, v_p, [x0; t0])
end

# ╔═╡ 32ab23a5-60b6-4f07-8c0f-cfa2e0a9d7e1
candidate.minimizer

# ╔═╡ 91f9aee4-6f80-4b02-ab09-1ace3759ea22
let dx = 1, dy = 1, npx = 100, npy = 100
	y0, x0, z0, t0 = candidate.minimizer
	xrange = range(x0 - dx, x0 + dx, npx)
	yrange = range(y0 - dy, y0 + dy, npy)
	sensors_readings = eachrow(Matrix(measurements[!, [:north, :east, :up, :toas]]))
	residues = reshape(
		[residue_rms([y, x, z0, t0], sensors_readings, v_p)
			for x = xrange, y = yrange],
		npx, npy)
	contour(xrange, yrange, residues, fill=:true, c=:heat)
	y1, x1 = templates[template_num, [:north, :east]]
	scatter!([x1], [y1], markershape=:circle, label="Template")
	scatter!([x0], [y0], markershape=:x, label="Event")
	plot!(title="X-Y Section")
end

# ╔═╡ 5c40564b-bb44-47f4-a98e-e09f161e2d68
let dx = 1, dz = 1, npx = 100, npz = 100
	y0, x0, z0, t0 = candidate.minimizer
	xrange = range(x0 - dx, x0 + dx, npx)
	zrange = range(z0 - dz, z0 + dz, npz)
	sensors_readings = eachrow(Matrix(measurements[!, [:north, :east, :up, :toas]]))
	residues = reshape(
		[TemplateMatching.residue_rms([y0, x, z, t0], sensors_readings, v_p)
			for x = xrange, z = zrange],
		npx, npz)
	contour(xrange, zrange, residues, fill=:true, c=:heat)
	z1, x1 = templates[template_num, [:up, :east]]
	scatter!([x1], [z1], markershape=:circle, label="Template")
	scatter!([x0], [z0], markershape=:x, label="Event")
	plot!(title="X-Z Section")
end

# ╔═╡ 102e4534-6be4-4436-be36-69726a393253
md" ## Statistics of the detections "

# ╔═╡ 33bc3279-c65f-47b5-acf4-60853c2bdb91
begin
	catalogue = DataFrame()
	samples = peaks .+ t_pre
	catalogue.datetimes = @. starttime + Microsecond(samples)
	catalogue.sample = samples
	catalogue.correlation = heights
	catalogue.rel_mag = [magnitude(values(data), template_data, peak .+ offsets) for peak in peaks]
	catalogue
end

# ╔═╡ 471f21ba-23d3-42f4-8d62-12fe2f36124f
histogram(collect(signal), label=nothing)

# ╔═╡ 6e2e9172-4cf5-45f3-a239-12c09b6ca3da
summarystats(collect(signal))

# ╔═╡ feebcd5f-91c4-4a70-a064-e642e3172614
histogram(catalogue[!, :correlation], label=nothing)

# ╔═╡ 476101af-0108-4628-9ac1-f233a456a869
summarystats(catalogue[!, :correlation])

# ╔═╡ 3be237b0-9f19-4e9a-aa13-c6355771b8e8
histogram(catalogue[!, :rel_mag], label=nothing)

# ╔═╡ Cell order:
# ╠═aa854f02-46bf-4a2b-b1e9-5fb52b99925c
# ╟─c2f22858-6d3c-4f92-a05b-659deae566f6
# ╠═8b4426b8-588b-49ad-83c5-a79ed698d704
# ╠═d47c1cf4-4844-4c66-8a21-b62ad7740dcf
# ╟─a86b6c8a-9b90-4388-8208-d257ec6951d2
# ╠═f3f17c65-dab6-46da-83a9-87aad1181524
# ╠═8271f7e6-d8ac-472a-9538-75d0afacf42e
# ╟─57b35604-727b-4eac-bcbd-ea302e4c79a2
# ╠═3c07c33b-7f91-4e2d-a27a-09a8924f328c
# ╠═b59dd2c0-4c3e-4c3d-a4d4-f8da9dcbe286
# ╠═76e2e659-c708-4a08-8263-def6c3dfa930
# ╟─b62d7318-3a6d-4442-b6cc-3166b87dff8a
# ╠═9de113e1-36ce-4251-9fa5-5baa27941a85
# ╠═96cf6e3b-1c2b-49f6-979d-0fce28a55a06
# ╠═f3139f19-e386-4bb7-8aca-fa54062d4093
# ╠═5d5f3ca7-260d-4dbb-bde1-59f2bd58c0c9
# ╟─6ac6cce3-09ae-421a-b502-b5c8807bf555
# ╠═be090485-7917-464a-8e8c-aec4a2c009d3
# ╟─6ab4a8da-4613-4004-807f-c36876285ebd
# ╠═41b0f95e-d806-4167-8199-7392717782f9
# ╠═95f2bd0b-7d8b-43b9-b2f9-cfa6628ec485
# ╠═356cc4c2-e565-4847-9286-cfc1f835654d
# ╠═730340bb-8176-4282-8559-8d60b1d28d57
# ╠═f07ab63f-dfb4-493c-a40c-6b3363a858c3
# ╟─fe42e7e6-2c56-4f11-aae6-ec06ad7a219f
# ╠═298010b3-c34c-40ba-9e1b-b79f7378f1d4
# ╠═a1fc5e84-48df-45ea-b6a1-0073ed9f8340
# ╠═c7f0582a-8997-4641-b607-0b95a7f73a6c
# ╠═074983c0-2747-4b09-a0c0-4d4adba24c0b
# ╠═cd5f912a-fafb-4947-bd70-318151c2c49c
# ╟─033473d3-6958-4e7b-8c95-a9e15d882118
# ╠═be80dd80-d39d-4338-bb76-c644ea9d970c
# ╠═aa23d53f-8333-4c94-8932-fac99ce9196d
# ╟─4ef908f4-d8dc-4159-97ad-fc8d4621a652
# ╟─368f7fc1-0afc-439f-905a-76af04e8ca8d
# ╠═342cbcff-5a05-4199-b5b5-6cc0fa70c202
# ╟─0d122aef-1767-4389-b043-97d07ecef1fe
# ╠═a7a589fc-4dcc-4137-a423-b8d2e7d41a0e
# ╟─91924a12-6ee4-4d1b-8c49-9b54142cf733
# ╟─abb8b867-39c2-40f2-9d7a-f40bacf8f550
# ╟─7d1f1903-8342-4715-8e88-3e3d063ea5db
# ╠═666c8d5f-c1c5-4cb4-9364-bf997b18a683
# ╟─96deaa3d-fc40-466f-bf1b-d8d3dd58dcd4
# ╟─cc0c3e0a-4208-4339-bea1-75396b2778bd
# ╟─b6848d59-e97e-4662-a2ee-78612694d6d8
# ╠═7ddcab3c-a100-44be-baa3-af28cc08b2d7
# ╠═341b9066-dcb3-4d99-b985-82e24c11f28e
# ╠═89b036db-5c47-4a17-8466-9062e8d5d452
# ╠═6d6d27c9-3b1a-4646-8420-61a775f7842f
# ╠═2911437d-4775-4a1b-819e-adf6e91acd61
# ╠═2e6b7b35-9512-406d-9a97-f1f94bea8732
# ╠═32ab23a5-60b6-4f07-8c0f-cfa2e0a9d7e1
# ╟─91f9aee4-6f80-4b02-ab09-1ace3759ea22
# ╟─5c40564b-bb44-47f4-a98e-e09f161e2d68
# ╟─102e4534-6be4-4436-be36-69726a393253
# ╠═33bc3279-c65f-47b5-acf4-60853c2bdb91
# ╠═471f21ba-23d3-42f4-8d62-12fe2f36124f
# ╠═6e2e9172-4cf5-45f3-a239-12c09b6ca3da
# ╠═feebcd5f-91c4-4a70-a064-e642e3172614
# ╠═476101af-0108-4628-9ac1-f233a456a869
# ╠═3be237b0-9f19-4e9a-aa13-c6355771b8e8


text:
# This file was generated, do not modify it. # hide
v[1] < v[2]


text:
@info "Test singleton pattern..."

println("Singleton object id: ", objectid(Singleton()))
println("Singleton object id: ", objectid(Singleton()))


text:
module test_4

using AdventOfCode2018.Day4
using Test
using Dates

test_input = """[1518-11-01 00:00] Guard #10 begins shift
[1518-11-01 00:05] falls asleep
[1518-11-01 00:25] wakes up
[1518-11-01 00:30] falls asleep
[1518-11-01 00:55] wakes up
[1518-11-01 23:58] Guard #99 begins shift
[1518-11-02 00:40] falls asleep
[1518-11-02 00:50] wakes up
[1518-11-03 00:05] Guard #10 begins shift
[1518-11-03 00:24] falls asleep
[1518-11-03 00:29] wakes up
[1518-11-04 00:02] Guard #99 begins shift
[1518-11-04 00:36] falls asleep
[1518-11-04 00:46] wakes up
[1518-11-05 00:03] Guard #99 begins shift
[1518-11-05 00:45] falls asleep
[1518-11-05 00:55] wakes up
"""

parsed_input = [
	(DateTime(1518, 11, 01, 00, 00), 10, false, false),
	(DateTime(1518, 11, 01, 00, 05), -1, true, false),
	(DateTime(1518, 11, 01, 00, 25), -1, false, true),
	(DateTime(1518, 11, 01, 00, 30), -1, true, false),
	(DateTime(1518, 11, 01, 00, 55), -1, false, true),
	(DateTime(1518, 11, 01, 23, 58), 99, false, false),
	(DateTime(1518, 11, 02, 00, 40), -1, true, false),
	(DateTime(1518, 11, 02, 00, 50), -1, false, true),
	(DateTime(1518, 11, 03, 00, 05), 10, false, false),
	(DateTime(1518, 11, 03, 00, 24), -1, true, false),
	(DateTime(1518, 11, 03, 00, 29), -1, false, true),
	(DateTime(1518, 11, 04, 00, 02), 99, false, false),
	(DateTime(1518, 11, 04, 00, 36), -1, true, false),
	(DateTime(1518, 11, 04, 00, 46), -1, false, true),
	(DateTime(1518, 11, 05, 00, 03), 99, false, false),
	(DateTime(1518, 11, 05, 00, 45), -1, true, false),
	(DateTime(1518, 11, 05, 00, 55), -1, false, true),
]

@testset "day 3, part 1" begin

	parsed = parse_input(test_input)
	@test parsed == parsed_input

	map = build_sleep_map(parsed)
	@test part_1(map) == 240
	@test part_2(map) == 4455

end

end


text:
module PSPlot

# load inkscape into path
if !Sys.isapple()
    error("PSPlot currently only support OS X")
else
    using Homebrew
    brew_path = joinpath(Homebrew.prefix(), "bin")
    ENV["PATH"] = string(ENV["PATH"], ":", brew_path)
end

using Reexport
ENV["MATPLOTLIBRC"] = joinpath(dirname(pathof(PSPlot)
       ), "assets")
@reexport using PyPlot # makes all of PyPlot available

export printfig
include("figure.jl")

end # module


text:
#################### Gelman, Rubin, and Brooks Diagnostics ####################

function _gelmandiag(
    psi::AbstractArray{<:Real,3};
    alpha::Real = 0.05
)
    niters, nparams, nchains = size(psi)
    nchains > 1 || error("Gelman diagnostic requires at least 2 chains")

    rfixed = (niters - 1) / niters
    rrandomscale = (nchains + 1) / (nchains * niters)

    S2 = mapslices(cov, psi, dims = [1, 2])
    W = dropdims(mapslices(mean, S2, dims = [3]), dims = 3)

    psibar = reshape(mapslices(mean, convert(Array, psi), dims = [1]), nparams, nchains)'
    B = niters .* cov(psibar)

    w = diag(W)
    b = diag(B)
    s2 = reshape(mapslices(diag, S2, dims = [1, 2]), nparams, nchains)'
    psibar2 = vec(mapslices(mean, psibar, dims = [1]))

    var_w = vec(mapslices(var, s2, dims = [1])) ./ nchains
    var_b = (2 / (nchains - 1)) .* b.^2
    var_wb = (niters / nchains) .*
        (diag(cov(s2, psibar.^2)) .- 2 .* psibar2 .* diag(cov(s2, psibar)))

    V = @. rfixed * w + rrandomscale * b
    var_V = rfixed^2 * var_w + rrandomscale^2 * var_b + 2 * rfixed * rrandomscale * var_wb

    df = @. 2 * V^2 / var_V

    B_df = nchains - 1
    W_df = @. 2 * w^2 / var_w

    estimates = Array{Float64}(undef, nparams)
    upperlimits = Array{Float64}(undef, nparams)

    q = 1 - alpha / 2
    for i in 1:nparams
        correction = (df[i] + 3) / (df[i] + 1)
        rrandom = rrandomscale * b[i] / w[i]

        estimates[i] = sqrt(correction * (rfixed + rrandom))

        if !isnan(rrandom)
            rrandom *= quantile(FDist(B_df, W_df[i]), q)
        end
        upperlimits[i] = sqrt(correction * (rfixed + rrandom))
    end

    return estimates, upperlimits, W, B
end

"""
    gelmandiag(chains::AbstractArray{<:Real,3}; kwargs...)

Gelman, Rubin and Brooks diagnostics.
"""
function gelmandiag(chains::AbstractArray{<:Real,3}; kwargs...)
    estimates, upperlimits = _gelmandiag(chains; kwargs...)

    return (psrf = estimates, psrfci = upperlimits)
end

"""
    gelmandiag_multivariate(chains::AbstractArray{<:Real,3}; kwargs...)

Multivariate Gelman, Rubin and Brooks diagnostics.
"""
function gelmandiag_multivariate(
    chains::AbstractArray{<:Real,3};
    kwargs...
)
    niters, nparams, nchains = size(chains)
    if nparams < 2
        error("computation of the multivariate potential scale reduction factor requires ",
              "at least two variables")
    end

    estimates, upperlimits, W, B = _gelmandiag(chains; kwargs...)

    rfixed = (niters - 1) / niters
    rrandomscale = (nchains + 1) / (nchains * niters)
    multivariate = rfixed + rrandomscale * LinearAlgebra.eigmax(W \ B)

    return (psrf = estimates, psrfci = upperlimits, psrfmultivariate = multivariate)
end

function gelmandiag(
    chains::Chains{<:Real};
    sections = _default_sections(chains),
    transform = false,
    alpha = 0.05,
    kwargs...
)
    # Subset the chain.
    _chains = Chains(chains, _clean_sections(chains, sections))

    # Compute the potential scale reduction factor.
    psi = transform ? link(_chains) : _chains.value.data
    results = gelmandiag(psi; alpha = alpha, kwargs...)

    # Create a named tuple with the results.
    sym = Symbol(100 * (1 - alpha / 2), :%)
    nt = (; :parameters => names(_chains), :psrf => results.psrf, sym => results.psrfci)

    return ChainDataFrame("Gelman, Rubin, and Brooks Diagnostic", nt)
end

function gelmandiag_multivariate(
    chains::Chains{<:Real};
    sections = _default_sections(chains),
    transform = true,
    alpha = 0.05,
    kwargs...
)
    # Subset the chain.
    _chains = Chains(chains, _clean_sections(chains, sections))

    # Compute the potential scale reduction factor.
    psi = transform ? link(_chains) : _chains.value.data
    results = gelmandiag_multivariate(psi; alpha = alpha, kwargs...)

    # Create a named tuple with the results.
    sym = Symbol(100 * (1 - alpha / 2), :%)
    nt = (; :parameters => names(_chains), :psrf => results.psrf, sym => results.psrfci)

    return ChainDataFrame("Gelman, Rubin, and Brooks Diagnostic", nt),
        results.multivariate
end


text:
using InterfaceMechanics
using Test

# @test hello("Julia") == "Hello, Julia"
# @test domath(2.0) ≈ 7.0

@testset "Test InterfaceMechanics.jl" begin
    @testset "test Couloumb interface" begin
        include("test_couloumb.jl")
    end
end


text:
using Compat, Compat.Test
push!(LOAD_PATH, joinpath(pwd(), "src"))
using JuML
using Test
using Random

data = bitrand(10)
boolvar = BoolVariate("boolvar", data)
@test length(boolvar) == 10
@test getname(boolvar) == "boolvar"

data2 = convert(Vector{Bool}, boolvar)
@test isa(data2, Vector{Bool})

max_stars_repo_name: Statfactory/JuML

max_issues_repo_name: Statfactory/JuML

max_forks_repo_name: Statfactory/JuML


text:
@defcomp LGcycle begin
    e_globalLGemissions=Parameter(index=[time],unit="Mtonne/year")
    e_0globalLGemissions=Parameter(unit="Mtonne/year", default=557.2112715473608)
    c_LGconcentration=Variable(index=[time],unit="ppbv")
    pic_preindustconcLG=Parameter(unit="ppbv", default=0.)
    exc_excessconcLG=Variable(unit="ppbv")
    c0_LGconcbaseyr=Parameter(unit="ppbv", default=0.11)
    re_remainLG=Variable(index=[time],unit="Mtonne")
    nte_natLGemissions=Variable(index=[time],unit="Mtonne/year")
    air_LGfractioninatm=Parameter(unit="%", default=100.)
    tea_LGemissionstoatm=Variable(index=[time],unit="Mtonne/year")
    teay_LGemissionstoatm=Variable(index=[time],unit="Mtonne/t")
    y_year=Parameter(index=[time],unit="year")
    y_year_0=Parameter(unit="year")
    res_LGatmlifetime=Parameter(unit="year", default=1000.)
    den_LGdensity=Parameter(unit="Mtonne/ppbv", default=100000.)
    stim_LGemissionfeedback=Parameter(unit="Mtonne/degreeC", default=0.)
    rtl_g0_baselandtemp=Parameter(unit="degreeC", default=0.9258270139190647)
    rtl_g_landtemperature=Parameter(index=[time],unit="degreeC")
    re_remainLGbase=Variable(unit="Mtonne")

    function run_timestep(p, v, d, t)
        if is_first(t)
            #eq.3 from Hope (2006) - natural emissions (carbon cycle) feedback, using global temperatures calculated in ClimateTemperature component
            nte0=p.stim_LGemissionfeedback*p.rtl_g0_baselandtemp
            v.nte_natLGemissions[t]=p.stim_LGemissionfeedback*p.rtl_g0_baselandtemp
            #eq.6 from Hope (2006) - emissions to atmosphere depend on the sum of natural and anthropogenic emissions
            tea0=(p.e_0globalLGemissions+nte0)*p.air_LGfractioninatm/100
            v.tea_LGemissionstoatm[t]=(p.e_globalLGemissions[t]+v.nte_natLGemissions[t])*p.air_LGfractioninatm/100
            v.teay_LGemissionstoatm[t]=(v.tea_LGemissionstoatm[t]+tea0)*(p.y_year[t]-p.y_year_0)/2
            #adapted from eq.1 in Hope(2006) - calculate excess concentration in base year
            v.exc_excessconcLG=p.c0_LGconcbaseyr-p.pic_preindustconcLG
            #Eq. 2 from Hope (2006) - base-year remaining emissions
            v.re_remainLGbase=v.exc_excessconcLG*p.den_LGdensity
            v.re_remainLG[t]=v.re_remainLGbase*exp(-(p.y_year[t]-p.y_year_0)/p.res_LGatmlifetime)+
                v.teay_LGemissionstoatm[t]*p.res_LGatmlifetime*(1-exp(-(p.y_year[t]-p.y_year_0)/p.res_LGatmlifetime))/(p.y_year[t]-p.y_year_0)
        else
            #eq.3 from Hope (2006) - natural emissions (carbon cycle) feedback, using global temperatures calculated in ClimateTemperature component
            #Here assume still using area-weighted average regional temperatures (i.e. land temperatures) for natural emissions feedback
            v.nte_natLGemissions[t]=p.stim_LGemissionfeedback*p.rtl_g_landtemperature[t-1]
            #eq.6 from Hope (2006) - emissions to atmosphere depend on the sum of natural and anthropogenic emissions
            v.tea_LGemissionstoatm[t]=(p.e_globalLGemissions[t]+v.nte_natLGemissions[t])*p.air_LGfractioninatm/100
            #eq.7 from Hope (2006) - average emissions to atm over time period
            v.teay_LGemissionstoatm[t]=(v.tea_LGemissionstoatm[t]+v.tea_LGemissionstoatm[t-1])*(p.y_year[t]-p.y_year[t-1])/2
            #eq.10 from Hope (2006) - remaining emissions in atmosphere
            v.re_remainLG[t]=v.re_remainLG[t-1]*exp(-(p.y_year[t]-p.y_year[t-1])/p.res_LGatmlifetime)+
                v.teay_LGemissionstoatm[t]*p.res_LGatmlifetime*(1-exp(-(p.y_year[t]-p.y_year[t-1])/p.res_LGatmlifetime))/(p.y_year[t]-p.y_year[t-1])
        end

    #eq.11 from Hope(2006) - LG concentration
        v.c_LGconcentration[t]=p.pic_preindustconcLG+v.exc_excessconcLG*v.re_remainLG[t]/v.re_remainLGbase

    end
end


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   Modia3D.ForceElements (Modia3D/ForceElements/_module.jl)
#

using Base.Meta:quot, isexpr

"""
    @forceElement forceName(arguments) begin ... end - Generate a function that instantiates a new forceElement

`forceName(arguments) begin ... end` is treated as the constructor function of a new
(mutable) struct that consists of all left-hand-side (scalar or vector) symbols
present between `begin ... end`.

# Examples

```
julia> include("$(Modia3D.path)/examples/dynamics/Simulate_PendulumWithDamper.jl");
```

```julia
using ModiaMath
using Modia3D

@forceElement Damper(; d=1.0) begin
    w   = ModiaMath.RealScalar("w",   causality=ModiaMath.Input,  numericType=ModiaMath.WR)
    tau = ModiaMath.RealScalar("tau", causality=ModiaMath.Output, numericType=ModiaMath.WR)
end;
function Modia3D.computeTorque(damper::Damper, sim::ModiaMath.SimulationState)
    damper.tau.value = -damper.d*damper.w.value
end;
@assembly PendulumWithDamper(;Lx = 1.0, Ly=0.2*Lx, Lz=0.2*Lx, m=1.0, g=9.81) begin
   world  = Modia3D.Object3D(Modia3D.CoordinateSystem(0.5*Lx))
   body   = Modia3D.Object3D(Modia3D.Solid(Modia3D.SolidBeam(Lx,Ly,Lz), Modia3D.MassProperties(m=m), vmat1))
   frame1 = Modia3D.Object3D(body; r=[-Lx/2, 0.0, 0.0])
   cyl    = Modia3D.Object3D(frame1,Modia3D.Cylinder(Ly/2,1.2*Ly; material=vmat2))

   # Connect pendulum to world with a damper in the joint
   rev    = Modia3D.Revolute(world, frame1)
   d      = Damper(d=0.2)
   damper = Modia3D.AdaptorForceElementToFlange(w=d.w, tau=d.tau)
   Modia3D.connect(damper, rev)
end
pendulum = PendulumWithDamper()
model = Modia3D.SimulationModel( pendulum )
result = ModiaMath.simulate!(model, stopTime=5.0, interval=0.1, tolerance=1e-4, log=true)
ModiaMath.plot(result, ["rev.phi", "rev.w", "rev.a", "rev.tau"])
```
"""
macro forceElement(head, top_ex)
   # @forceElement Name(arguments) begin ... end
   #    is mapped to:
   # mutable struct Name <: ModiaMath.AbstractComponentWithVariables
   #    < declarations >
   #    function Name(<arguments>)
   #       < equations top_ex >
   #       _internal  = ModiaMath.ComponentInternal(:Name, nothing, computeVariables!)
   #       _force = new(_internal, <declarations>)
   #       ModiaMath.initComponent(_internal, name1, :name1)
   #       ModiaMath.initComponent(_internal, name2, :name2)
   #         ...
   #       return _force
   #    end
   # end

   if typeof(head) == Symbol
      force_symbol = head
      fc = Expr(:call, head)
   else
      @assert isexpr(head, :call)
      force_symbol = head.args[1]::Symbol
      fc = head
   end

   # Build:  internalComponentName = ModiaMath.Component(component_name)
   # internalComponentName = gensym()
   # push!(eq, :( a = ModiaMath.Component($component_name)) )

   # Inspect top_ex find
   #   symbol = < statements >
   # and store "symbol" in "names"
   this = gensym("this")
   eqInternal = []
   eq      = []
   push!(eqInternal,:( $this = new( ModiaMath.ComponentInternal($(quot(force_symbol))) ) ))
   namesSet = DataStructures.OrderedSet{Symbol}() # Set of left-hand side symbols
   helpArray = []

   # Add keyword arguments to namesSet
   if length(fc.args) <= 1
      # No keyword arguments
   elseif typeof(fc.args[2]) == Symbol
      # Positional arguments present, but no keyword arguments
   elseif fc.args[2].head == :parameters
      # Keyword arguments present
      kwargs = fc.args[2].args
      # println("... keyword arguments = ", kwargs)
      for kw in kwargs
         # println("... kw = ", kw, ", kw.head = ", kw.head, ", kw.args = ", kw.args, ", name = ", kw.args[1])
         symbol = kw.args[1]
         if in(symbol, namesSet)
            error("\nfrom @forceElement ", force_symbol, ": keyword argument",
                  "\n   ", kw,
                  "\nis defined twice.")
         end
         push!(namesSet, symbol)
         push!(helpArray, :($this.$symbol = $symbol))
      end
      # push!(fcp.args[2].args, Expr(:kw, :sceneOptions, :nothing))
   else
      error("Unknown error in forces.jl: typeof(fc.args[2]) = ", typeof(fc.args[2]))
   end

   for ex in top_ex.args
      # dump(ex)
      push!(eq, ex)
      if isexpr(ex, :(=), 2) && typeof(ex.args[1]) == Symbol
         symbol = ex.args[1]
         if in(symbol, namesSet)
            error("\nfrom @forceElement ", force_symbol, ": Left hand side symbol in",
                  "\n   ", ex,
                  "\nis defined twice.")
         end
         push!(namesSet, symbol)
         push!(eq, :( ModiaMath.initComponent!($this,$symbol,$(quot(symbol)) ) ) )
     end
   end
   names = collect(namesSet)

   # Build:
   #   mutable struct force_symbol <: ModiaMath.AbstractComponentWithVariables
   #      _internal::ModiaMath.ComponentInternal
   #      names[1]
   #      names[2]
   #       ...
   #      function force_symbol(...)
   #         $(eq...)
   #         ModiaMath.initComponentInternal!( new(...) )
   #      end
   #   end
   code = :(mutable struct $force_symbol <: Modia3D.AbstractForceTorque; end)
   # println("code.args = ", code.args)
   code.args[3].args = vcat( :( _internal::ModiaMath.ComponentInternal ), names )

   code2 = code.args[3].args
   push!(code2, :(function fc(); end) )
   i = length(code2)
   code2[i].args[1] = fc

   # Generate constructor function
   code2[i].args[2] = quote
      $(eqInternal...)
      $(helpArray...)
      $(eq...)
      return $this
   end

   # Print generated function (without "#....\n")
   # println( replace(sprint(showcompact,code) , r"# .*\n", "\n") )
   # println( sprint(showcompact,bcode) )

   return esc(code)
end


text:
# ---
# title: 1553. Minimum Number of Days to Eat N Oranges
# id: problem1553
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/minimum-number-of-days-to-eat-n-oranges/description/>
# hidden: true
# ---
#
# There are `n` oranges in the kitchen and you decided to eat some of these
# oranges every day as follows:
#
#   * Eat one orange.
#   * If the number of remaining oranges (`n`) is divisible by 2 then you can eat  n/2 oranges.
#   * If the number of remaining oranges (`n`) is divisible by 3 then you can eat  2*(n/3) oranges.
#
# You can only choose one of the actions per day.
#
# Return the minimum number of days to eat `n` oranges.
#
#
#
# **Example 1:**
#
#
#
#     Input: n = 10
#     Output: 4
#     Explanation: You have 10 oranges.
#     Day 1: Eat 1 orange,  10 - 1 = 9.
#     Day 2: Eat 6 oranges, 9 - 2*(9/3) = 9 - 6 = 3. (Since 9 is divisible by 3)
#     Day 3: Eat 2 oranges, 3 - 2*(3/3) = 3 - 2 = 1.
#     Day 4: Eat the last orange  1 - 1  = 0.
#     You need at least 4 days to eat the 10 oranges.
#
#
# **Example 2:**
#
#
#
#     Input: n = 6
#     Output: 3
#     Explanation: You have 6 oranges.
#     Day 1: Eat 3 oranges, 6 - 6/2 = 6 - 3 = 3. (Since 6 is divisible by 2).
#     Day 2: Eat 2 oranges, 3 - 2*(3/3) = 3 - 2 = 1. (Since 3 is divisible by 3)
#     Day 3: Eat the last orange  1 - 1  = 0.
#     You need at least 3 days to eat the 6 oranges.
#
#
# **Example 3:**
#
#
#
#     Input: n = 1
#     Output: 1
#
#
# **Example 4:**
#
#
#
#     Input: n = 56
#     Output: 6
#
#
#
#
# **Constraints:**
#
#   * `1 <= n <= 2*10^9`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
# differentation rule for the iterative algorithm in the inverse of `PlanarLayer`
ChainRulesCore.@scalar_rule(
    find_alpha(wt_y::Real, wt_u_hat::Real, b::Real),
    @setup(
        x = inv(1 + wt_u_hat * sech(Ω + b)^2),
    ),
    (x, - tanh(Ω + b) * x, x - 1),
)


text:
export pretty_print_txt

"""
    pretty_print_txt( bell_games :: Vector{BellGame}, filename :: String )

Prints a set of `BellGame`'s to `filename.txt` in a human-readable form.
"""
function pretty_print_txt(bell_games :: Vector{BellGame}, filename :: String)
    ext = occursin(r"\.txt$", filename) ? "" : ".txt"
    filepath = filename * ext
    open(filepath, "w") do io
        for id in 1:length(bell_games)
            game = bell_games[id]
            print(io, "(", id ,")  ", game.β, " ≥ ")
            print(io, "\n")

            for game_row in eachrow(game)
                println(io, "\t", game_row)
            end

            print(io, "\n")
        end
    end
end


text:
module Manifolds

import Base:
    Array,
    +,
    -,
    *,
    \,
    /,
    ^,
    ==,
    angle,
    axes,
    convert,
    copy,
    copyto!,
    dataids,
    eltype,
    exp,
    getindex,
    identity,
    inv,
    isapprox,
    length,
    log,
    one,
    promote_rule,
    setindex!,
    similar,
    show,
    similar,
    size,
    transpose,
    zero
import Distributions: _rand!, support
import LinearAlgebra: cross, det, Diagonal, dot, mul!, norm, I, UniformScaling
import ManifoldsBase:
    allocate,
    allocate_result,
    allocate_result_type,
    array_value,
    base_manifold,
    check_manifold_point,
    check_tangent_vector,
    distance,
    exp,
    exp!,
    geodesic,
    injectivity_radius,
    inner,
    isapprox,
    is_manifold_point,
    is_tangent_vector,
    is_decorator_manifold,
    inverse_retract,
    inverse_retract!,
    log,
    log!,
    manifold_dimension,
    norm,
    number_eltype,
    project_point,
    project_point!,
    project_tangent,
    project_tangent!,
    representation_size,
    retract,
    retract!,
    shortest_geodesic,
    vector_transport_along,
    vector_transport_along!,
    vector_transport_direction,
    vector_transport_direction!,
    vector_transport_to,
    vector_transport_to!,
    zero_tangent_vector,
    zero_tangent_vector!
import Random: rand
import Statistics: mean, mean!, median, median!, std, var
import StatsBase: kurtosis, mean_and_std, mean_and_var, moment, skewness

using Base.Iterators: repeated
using Distributions
using Einsum: @einsum
using FiniteDifferences
using HybridArrays
using LinearAlgebra
using LightGraphs
using LightGraphs: AbstractGraph
using ManifoldsBase: CoTVector, Manifold, MPoint, TVector, DefaultManifold
using ManifoldsBase:
    ArrayCoTVector, ArrayManifold, ArrayMPoint, ArrayTVector, ArrayCoTVector
using ManifoldsBase: AbstractRetractionMethod, ExponentialRetraction
using ManifoldsBase: AbstractInverseRetractionMethod, LogarithmicInverseRetraction
using ManifoldsBase: AbstractVectorTransportMethod, ParallelTransport, ProjectionTransport
using Markdown: @doc_str
using Random: AbstractRNG
using Requires
using SimpleWeightedGraphs: AbstractSimpleWeightedGraph, get_weight
using StaticArrays
using StatsBase: AbstractWeights, ProbabilityWeights, values, varcorrection
using UnsafeArrays

"""
    AbstractEstimationMethod

Abstract type for defining statistical estimation methods.
"""
abstract type AbstractEstimationMethod end

@doc doc"""
    hat(M::Manifold, x, vⁱ)

Given a basis $e_i$ on the tangent space at a point $x$ and tangent
component vector $v^i$, compute the equivalent vector representation
$v=v^i e_i$, where Einstein summation notation is used:

````math
\wedge: v^i \mapsto v^i e_i
````

For array manifolds, this converts a vector representation of the tangent
vector to an array representation. The [`vee`](@ref) map is the `hat` map's
inverse.
"""
function hat(M::Manifold, x, vⁱ)
    v = allocate_result(M, hat, x, vⁱ)
    return hat!(M, v, x, vⁱ)
end

function hat!(M::Manifold, v, x, vⁱ)
    is_decorator_manifold(M) === Val(true) && return hat!(base_manifold(M), v, x, vⁱ)
    error("hat! operator not defined for manifold $(typeof(M)), array $(typeof(v)), point $(typeof(x)), and vector $(typeof(vⁱ))")
end

@doc doc"""
    vee(M::Manifold, x, v)

Given a basis $e_i$ on the tangent space at a point $x$ and tangent
vector $v$, compute the vector components $v^i$, such that $v = v^i e_i$, where
Einstein summation notation is used:

````math
\vee: v^i e_i \mapsto v^i
````

For array manifolds, this converts an array representation of the tangent
vector to a vector representation. The [`hat`](@ref) map is the `vee` map's
inverse.
"""
function vee(M::Manifold, x, v)
    vⁱ = allocate_result(M, vee, x, v)
    return vee!(M, vⁱ, x, v)
end

function vee!(M::Manifold, vⁱ, x, v)
    is_decorator_manifold(M) === Val(true) && return vee!(base_manifold(M), vⁱ, x, v)
    error("vee! operator not defined for manifold $(typeof(M)), vector $(typeof(vⁱ)), point $(typeof(x)), and array $(typeof(v))")
end

function allocate_result(M::Manifold, f::typeof(vee), x, v)
    T = allocate_result_type(M, f, (x, v))
    return allocate(x, T, manifold_dimension(M))
end
function allocate_result(M::Manifold, f::typeof(vee), x::StaticArray, v)
    T = allocate_result_type(M, f, (x, v))
    return allocate(x, T, Size(manifold_dimension(M)))
end

"""
    PolarRetraction <: AbstractRetractionMethod

Retractions that are based on singular value decompositions of the matrix / matrices
for point and tangent vector on a [`Manifold`](@ref)
"""
struct PolarRetraction <: AbstractRetractionMethod end

"""
    ProjectionRetraction <: AbstractRetractionMethod

Retractions that are based on projection and usually addition in the embedding.
"""
struct ProjectionRetraction <: AbstractRetractionMethod end

"""
    QRRetraction <: AbstractRetractionMethod

Retractions that are based on a QR decomposition of the
matrix / matrices for point and tangent vector on a [`Manifold`](@ref)
"""
struct QRRetraction <: AbstractRetractionMethod end

"""
    PolarInverseRetraction <: AbstractInverseRetractionMethod

Inverse retractions that are based on a singular value decomposition of the
matrix / matrices for point and tangent vector on a [`Manifold`](@ref)
"""
struct PolarInverseRetraction <: AbstractInverseRetractionMethod end

"""
    ProjectionInverseRetraction <: AbstractInverseRetractionMethod

Inverse retractions that are based on a projection (or its inversion).
"""
struct ProjectionInverseRetraction <: AbstractInverseRetractionMethod end

"""
    QRInverseRetraction <: AbstractInverseRetractionMethod

Inverse retractions that are based on a QR decomposition of the
matrix / matrices for point and tangent vector on a [`Manifold`](@ref)
"""
struct QRInverseRetraction <: AbstractInverseRetractionMethod end

include("utils.jl")
include("numbers.jl")
include("orthonormal_bases.jl")
include("autodiff.jl")
include("SizedAbstractArray.jl")

include("manifolds/VectorBundle.jl")

include("distributions.jl")
include("projected_distribution.jl")
include("product_representations.jl")

include("manifolds/MetricManifold.jl")
include("manifolds/ProductManifold.jl")
include("manifolds/PowerManifold.jl")
include("manifolds/GraphManifold.jl")

include("manifolds/CholeskySpace.jl")
include("manifolds/Circle.jl")
include("manifolds/Euclidean.jl")
include("manifolds/FixedRankMatrices.jl")
include("manifolds/Grassmann.jl")
include("manifolds/Hyperbolic.jl")
include("manifolds/Rotations.jl")
include("manifolds/Stiefel.jl")
include("manifolds/Sphere.jl")
include("manifolds/Symmetric.jl")
include("manifolds/SymmetricPositiveDefinite.jl")
include("manifolds/SymmetricPositiveDefiniteLinearAffine.jl")
include("manifolds/SymmetricPositiveDefiniteLogCholesky.jl")
include("manifolds/SymmetricPositiveDefiniteLogEuclidean.jl")
include("manifolds/Torus.jl")

include("groups/group.jl")
include("groups/group_action.jl")
include("groups/group_operation_action.jl")
include("groups/array_manifold.jl")
include("groups/product_group.jl")
include("groups/semidirect_product_group.jl")

include("groups/translation_group.jl")
include("groups/special_orthogonal.jl")
include("groups/circle_group.jl")

include("groups/translation_action.jl")
include("groups/rotation_action.jl")

include("groups/special_euclidean.jl")

include("statistics.jl")

function __init__()
    @require ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" begin
        using .ForwardDiff
        include("forward_diff.jl")
    end

    @require OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" begin
        using .OrdinaryDiffEq: ODEProblem, AutoVern9, Rodas5, solve
        include("ode.jl")
    end
end
#
export CoTVector, Manifold, MPoint, TVector, Manifold
export Euclidean,
    CholeskySpace,
    Circle,
    FixedRankMatrices,
    Grassmann,
    Hyperbolic,
    Rotations,
    Sphere,
    Stiefel,
    SymmetricMatrices,
    SymmetricPositiveDefinite,
    Torus
export SVDMPoint, UMVTVector, AbstractNumbers, ℝ, ℂ, ℍ
# decorator manifolds
export ArrayManifold, ArrayMPoint, ArrayTVector, ArrayCoTVector
export CotangentBundle,
    CotangentSpaceAtPoint, CotangentBundleFibers, CotangentSpace, FVector
export AbstractPowerManifold, AbstractPowerRepresentation,
    MultidimentionalArrayPowerRepresentation, NestedPowerRepresentation, PowerManifold
export ProductManifold
export GraphManifold, GraphManifoldType, VertexManifold, EdgeManifold
export ProjectedPointDistribution, ProductRepr, TangentBundle, TangentBundleFibers
export TangentSpace, TangentSpaceAtPoint, VectorSpaceAtPoint, VectorSpaceType, VectorBundle
export VectorBundleFibers
export AbstractVectorTransportMethod, ParallelTransport, ProjectedPointDistribution
export Metric,
    RiemannianMetric,
    LorentzMetric,
    MinkowskiMetric,
    EuclideanMetric,
    MetricManifold,
    LinearAffineMetric,
    LogEuclideanMetric,
    LogCholeskyMetric,
    PowerMetric,
    ProductMetric
export AbstractVectorTransportMethod, ParallelTransport, ProjectionTransport
export AbstractRetractionMethod, QRRetraction, PolarRetraction, ProjectionRetraction
export AbstractInverseRetractionMethod,
    QRInverseRetraction, PolarInverseRetraction, ProjectionInverseRetraction
export AbstractEstimationMethod,
    GradientDescentEstimation,
    CyclicProximalPointEstimation,
    GeodesicInterpolation,
    GeodesicInterpolationWithinRadius
export allocate,
    allocate_result,
    base_manifold,
    bundle_projection,
    christoffel_symbols_first,
    christoffel_symbols_second,
    christoffel_symbols_second_jacobian,
    complex_dot,
    det_local_metric,
    distance,
    einstein_tensor,
    exp,
    exp!,
    flat,
    flat!,
    gaussian_curvature,
    geodesic,
    hat,
    hat!,
    incident_log,
    injectivity_radius,
    inner,
    inverse_local_metric,
    inverse_retract,
    inverse_retract!,
    isapprox,
    is_decorator_manifold,
    is_default_metric,
    is_manifold_point,
    is_tangent_vector,
    isapprox,
    inner,
    kurtosis,
    local_metric,
    local_metric_jacobian,
    log,
    log!,
    log_local_metric_density,
    manifold_dimension,
    metric,
    mean,
    mean!,
    mean_and_var,
    mean_and_std,
    median,
    median!,
    moment,
    norm,
    normal_tvector_distribution,
    number_eltype,
    one,
    project_point,
    project_point!,
    project_tangent,
    project_tangent!,
    projected_distribution,
    real_dimension,
    ricci_curvature,
    ricci_tensor,
    representation_size,
    retract,
    retract!,
    riemann_tensor,
    sharp,
    sharp!,
    shortest_geodesic,
    skewness,
    std,
    sym_rem,
    submanifold,
    submanifold_component,
    submanifold_components,
    var,
    vector_space_dimension,
    vector_transport_along,
    vector_transport_along!,
    vector_transport_direction,
    vector_transport_direction!,
    vector_transport_to,
    vector_transport_to!,
    vee,
    vee!,
    zero_vector,
    zero_vector!,
    zero_tangent_vector,
    zero_tangent_vector!,
    ×
# Lie group types & functions
export AbstractGroupAction,
    AbstractGroupOperation,
    AbstractGroupManifold,
    ActionDirection,
    AdditionOperation,
    CircleGroup,
    MultiplicationOperation,
    GroupManifold,
    GroupOperationAction,
    Identity,
    LeftAction,
    ProductGroup,
    ProductOperation,
    RightAction,
    RotationAction,
    SemidirectProductGroup,
    SpecialEuclidean,
    SpecialOrthogonal,
    TranslationGroup,
    TranslationAction
export apply,
    apply!,
    apply_diff,
    apply_diff!,
    base_group,
    center_of_orbit,
    compose,
    compose!,
    direction,
    g_manifold,
    identity,
    identity!,
    inv,
    inv!,
    inverse_apply,
    inverse_apply!,
    inverse_apply_diff,
    inverse_apply_diff!,
    inverse_translate,
    inverse_translate!,
    inverse_translate_diff,
    inverse_translate_diff!,
    optimal_alignment,
    optimal_alignment!,
    switch_direction,
    translate,
    translate!,
    translate_diff,
    translate_diff!
# Orthonormal bases
export AbstractBasis,
    AbstractOrthonormalBasis,
    AbstractPrecomputedOrthonormalBasis,
    ArbitraryOrthonormalBasis,
    DiagonalizingOrthonormalBasis,
    PrecomputedDiagonalizingOrthonormalBasis,
    PrecomputedOrthonormalBasis,
    PrecomputedProductOrthonormalBasis,
    ProjectedOrthonormalBasis
export get_basis, get_coordinates, get_vector, get_vectors, number_system
end # module


text:
function find_texfiles(rootdir::String)
    texfiles = String[]
    for dir in readdir(rootdir)
        isdir(dir) || continue
        files = find(readdir(dir)) do file
            endswith(file, ".tex")
        end
        if length(files) == 1
            f = joinpath(rootdir, dir, files[1])
            println(f)
            push!(texfiles, f)
        end
    end
    texfiles
end

rootdir = ""
find_texfiles(rootdir)

max_stars_repo_name: paperai/LightJATS

max_issues_repo_name: paperai/LightJATS

max_forks_repo_name: paperai/LightJATS


text:
using NBodyIPs, JuLIP, BenchmarkTools, StaticArrays
using Test

using NBodyIPs: ClusterBLDesc, transform, transform_d, fcut, fcut_d,
                invariants, invariants_d, invariants_ed
using NBodyIPs.Polys: NBPoly

@info("Setting up the test systems ...")
r0 = rnn(:W)
# TRANSFORM = "r -> exp( - 3 * ((r/$r0) - 1))"
TRANSFORM = "r -> 1/r"
rcut2 = 3.1 * r0
D2 = ClusterBLDesc(TRANSFORM, (:cos, 0.66*rcut2, rcut2) )
rcut3 = 2.3 * r0
D3 = ClusterBLDesc(TRANSFORM, (:cos, 0.66*rcut3, rcut3) )
rcut4 = 1.9 * r0
D4 = ClusterBLDesc(TRANSFORM, (:cos, 0.66*rcut4, rcut4) )

DD = [nothing, D2, D3, D4]

# TEST 1:
B2 = nbpolys(2, D2, 23)
c2 = rand(length(B2))
B3 = nbpolys(3, D3, 12)
c3 = rand(length(B3))
B4 = nbpolys(4, D4, 8)
c4 = rand(length(B4))

B4l = nbpolys(4, D4, 14)
c4l = rand(length(B4l))

IP2 = NBodyIP(B2, c2)
IP3 = NBodyIP([B2;B3], [c2;c3])
IP4 = NBodyIP([B2;B3;B4], [c2;c3;c4])
IP4l = NBodyIP([B2;B3;B4l], [c2;c3;c4l])

IP2f = fast(IP2)
IP3f = fast(IP3)
IP4f = fast(IP4)
IP4lf = fast(IP4l)

at = bulk(:W, cubic=true)*4
@btime forces($IP2f, $at)
@btime forces($IP3f, $at)
@btime forces($IP4f, $at)
@btime forces($IP4lf, $at)

# println((@elapsed forces(IP2f, at)) / length(at) * 1_000)
# println((@elapsed forces(IP3f, at)) / length(at) * 1_000)
# println((@elapsed forces(IP4f, at)) / length(at) * 1_000)
# println((@elapsed forces(IP4lf, at)) / length(at) * 1_000)

println(9.7490 / length(at))
println(20.673 / length(at))
println(28.727 / length(at))
println(37.033 / length(at))

max_stars_repo_name: cortner/ManyBodyIPs

max_issues_repo_name: cortner/ManyBodyIPs

max_forks_repo_name: cortner/ManyBodyIPs


text:
using HTTP
HTTP.download("https://raw.githubusercontent.com/joan2937/pigpio/master/pigpio.h", "pigpio.h")

using Clang
Clang.init(
    headers        = ["pigpio.h"],
    output_file    =  "pigpio.jl",
    common_file    =  "pigpio_common.jl",
    clang_includes = [CLANG_INCLUDE],
    header_wrapped = (root, current)->root == current,
    cursor_wrapped = (cursorname, cursor) -> !contains(string(cursor), "Thread")
) |> run


text:
using NeXLCore
"""
An abstract structure that implements

    NeXLCore.bremsstrahlung(::Type{<:NeXLBremsstrahlung}, e::AbstractFloat, e0::AbstractFloat, elm::Element; kwargs...)
"""
abstract type NeXLBremsstrahlung <: NeXLAlgorithm end

"""
@article{kramers1923,
    author = { Kramers, H. A. },
    title = {XCIII. On the theory of X-ray absorption and of the continuous X-ray spectrum},
    journal = {The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science},
    volume = {46},
    number = {275},
    pages = {836-871},
    year  = {1923},
    publisher = {Taylor & Francis},
    doi = {10.1080/14786442308565244},
    URL = { https://doi.org/10.1080/14786442308565244 },
    eprint = { https://doi.org/10.1080/14786442308565244}
}
"""
struct Kramers1923 <: NeXLBremsstrahlung end

"""
@InProceedings{ lifshin1974,
    author={Lifshin, E.},
    journal="Proceedings Ninth National Conference on Electron Probe Analysis, Ottowa, CA",
    year={1974},
    volume={174},
    number={3},
}
"""
struct Lifshin1974 <: NeXLBremsstrahlung end

"""
@article{reed1975shape,
  title={The shape of the continuous X-ray spectrum and background corrections for energy-dispersive electron microprobe analysis},
  author={Reed, SJB},
  journal={X-Ray Spectrometry},
  volume={4},
  number={1},
  pages={14--17},
  year={1975},
  publisher={Wiley Online Library}
}
"""
struct Reed1975 <: NeXLBremsstrahlung end

"""
@article{ smithgoldtomlinson1975
    title="The atomic number dependence of the X‐ray continuum intensity and the practical calculation of background in energy dispersive electron microprobe analysis",
    author = {Smith, D. G. W.  and Gold, C. M. and Tomlinson, D. A.},
    journal={X-ray Spectrometry},
    volume={4},
    pages={149-156},
    year={1975},
    doi={10.1002/xrs.1300040311},
}
"""
struct Smith1975 <: NeXLBremsstrahlung end

"""
@article{small1987modeling,
  title={Modeling of the bremsstrahlung radiation produced in pure-element targets by 10--40 keV electrons},
  author={Small, John A and Leigh, Stefan D and Newbury, Dale E and Myklebust, Robert L},
  journal={Journal of applied physics},
  volume={61},
  number={2},
  pages={459--469},
  year={1987},
  publisher={American Institute of Physics}
}
"""
struct Small1987 <: NeXLBremsstrahlung end

"""
@article{trincavelli1998model,
  title={Model for the bremsstrahlung spectrum in EPMA. Application to standardless quantification},
  author={Trincavelli, Jorge and Castellano, Gustavo and Riveros, J Alberto},
  journal={X-Ray Spectrometry: An International Journal},
  volume={27},
  number={2},
  pages={81--86},
  year={1998},
  publisher={Wiley Online Library}
}
"""
struct Trincavelli1997 <: NeXLBremsstrahlung end

"""
@article{castellano2004analytical,
  title={Analytical model for the bremsstrahlung spectrum in the 0.25--20 keV photon energy range},
  author={Castellano, Gustavo and Osan, Janos and Trincavelli, Jorge},
  journal={Spectrochimica Acta Part B: Atomic Spectroscopy},
  volume={59},
  number={3},
  pages={313--319},
  year={2004},
  publisher={Elsevier}
}
"""
struct Castellano2004a <: NeXLBremsstrahlung end
struct Castellano2004b <: NeXLBremsstrahlung end

bremsstrahlung(::Type{Kramers1923}, e::AbstractFloat, e0::AbstractFloat, elm::Element) =
    e < e0 ? z(elm) * (e0 - e) / e : 0.0

bremsstrahlung(
    ::Type{Lifshin1974},
    e::AbstractFloat,
    e0::AbstractFloat,
    elm::Element;
    a = 0.0001,
) = bremsstrahlung(Kramers1923, e, e0, elm) * (1.0 - 1.0e-3 * a * (e0 - e))

bremsstrahlung(
    ::Type{Reed1975},
    e::AbstractFloat,
    e0::AbstractFloat,
    elm::Element;
    b = 0.0001,
) = e < e0 ? z(elm) * (1.0e-3 * (e0 - e)) / ((1.0e-3 * e)^(1.0 + b)) : 0.0

function bremsstrahlung(
    ::Type{Smith1975},
    e::AbstractFloat,
    e0::AbstractFloat,
    elm::Element,
)
    n(e, z) = 1.0e-3 * e * (0.0739 - 0.0051 * log(z)) + 1.6561 - 0.115 * log(z)
    x(e, z) = 1.76 - 0.00145 * z / (1.0e-3 * e)
    zz = Float64(z(elm))
    return ((e0 - e) / e)^x(e, zz) * zz^n(e, zz)
end

function bremsstrahlung(
    ::Type{Small1987},
    e::AbstractFloat,
    e0::AbstractFloat,
    elm::Element,
)
    M(e0) = 0.00599e-3 * e0 + 1.05
    B(e0) = -0.0322e-3 * e0
    return e < e0 ? (z(elm) * ((e0 - e) / e))^M(e0) * (1.0e-3 * e)^B(e0) : 0.0
end

bremsstrahlung(::Type{Trincavelli1997}, e::AbstractFloat, e0::AbstractFloat, elm::Element) =
    e < e0 ?
    (sqrt(z(elm)) * (e0 - e) / e) * (
        -54.86 - 1.072e-3 * e +
        0.2835e-3 * e0 +
        30.4 * log(z(elm)) +
        875.0 / (z(elm)^2 * (1.0e-3 * e0)^0.08)
    ) : 0.0

bremsstrahlung(::Type{Castellano2004a}, e::AbstractFloat, e0::AbstractFloat, elm::Element) =
    e < e0 ?
    (sqrt(z(elm)) * (e0 - e) / e) *
    (
        -73.90 - 1.2446 * (1.0e-3 * e) +
        36.502 * log(z(elm)) +
        (148.5 * (1.0e-3 * e0)^0.1239) / z(elm)
    ) * #
    (1.0 + (-0.006624 + 0.0002906 * (1.0e-3 * e0)) * z(elm) / (1.0e-3 * e)) : 0.0

function bremsstrahlung(
    ::Type{Castellano2004b},
    e::AbstractFloat,
    e0::AbstractFloat,
    elm::Element,
)
    zz, E0, E = z(elm), e0 / 1000.0, e / 1000.0
    return e < e0 ?
           (
        (-E + E0) *
        sqrt(zz) *
        (1.0 + ((-0.006626 + 0.0002906 * E0) * zz) / E) *
        (-77.28317356370013 + (148.5 * E0^0.1293) / zz + 36.502 * log(zz))
    ) / E : 0.0
end

bremsstrahlung(e::AbstractFloat, e0::AbstractFloat, elm::Element) =
    bremsstrahlung(Castellano2004a, e, e0, elm)

function bremsstrahlung(
    ty::Type{<:NeXLBremsstrahlung},
    e::AbstractFloat,
    e0::AbstractFloat,
    mat::Material;
    kwargs...,
)
    af = atomicfraction(mat) # According to Trincavelli1997
    return mapreduce(
        elm -> af[elm] * bremsstrahlung(ty, e, e0, elm; kwargs...),
        +,
        keys(af),
    )
end

# Evaluating these models
# The weave script "NeXLSpectrum/weave/continuummodel.jmd" can be used to compare the NeXLBremsstrahlung algorithms
# against a K412 spectrum and the associated standard spectra.  In addition to the continuum model, it is also necessary
# to provide a detector efficiency model and matrix correction model.  Often a conventional Z⋅A matrix correction
# model is adapted by replacing the edge energy with the the Bremsstrahlung emission energy.  This is correct to
# the degree which the ϕ(ρz)-curve models share the same emission profile as the continuum emission.

# Evaluating the models I find that Castellano2004a, Trincavelli1997 work well with the Riveros1993 matrix correction
# algorithm and the AP33Tabulation window.  Smith1975 works surprisigly well with the CitZAF matrix correction model.
# Other old models based on Si(Li) data tend to not do too well at lower energies.  This shouldn't surprise anyone
# as these models were often based on data from Be window detectors.  Castellano2004a and Trincavelli1997 were
# designed around the Riveros1993 matrix correction model and don't perform well using CitZAF.

# My current recommendation is Castellano2004a & Riveros1993.


text:
using IntelVectorMath, SpecialFunctions
using AcuteBenchmark

cd(@__DIR__)
rootPath = pwd()
################################################################
# Dimensions
oneArgDims = [10 20 50 100 200 500 700 1000 5000 10000]

twoArgDims = [10 20 50 100 200 500 700 1000 5000 10000; # arg1
              10 20 50 100 200 500 700 1000 5000 10000] # arg2

################################################################
## Real Functions

typesReal=[Float32, Float64]

configsRealBase = FunbArray([

    # oneArgDims
    Funb(acos, [(-1, 1)], typesReal, oneArgDims),
    Funb(asin, [(-1, 1)], typesReal, oneArgDims),
    Funb(atan, [(-50, 50)], typesReal, oneArgDims),
    Funb(cos, [(-1000, 1000)], typesReal, oneArgDims),
    Funb(sin, [(-1000, 1000)], typesReal, oneArgDims),
    Funb(tan, [(-1000, 1000)], typesReal, oneArgDims),
    Funb(acosh, [(1, 1000)], typesReal, oneArgDims),
    Funb(asinh, [(-1000, 1000)], typesReal, oneArgDims),
    Funb(atanh, [(-1, 1)], typesReal, oneArgDims),
    Funb(cosh, [(0, 89.415985f0)], typesReal, oneArgDims),
    Funb(sinh, [(-89.415985f0, 89.415985f0)], typesReal, oneArgDims),
    Funb(tanh, [(-8.66434f0, 8.66434f0)], typesReal, oneArgDims),
    Funb(cbrt, [(-1000, 1000)], typesReal, oneArgDims),
    Funb(sqrt, [(0, 1000)], typesReal, oneArgDims),
    Funb(exp, [(-88.72284f0, 88.72284f0)], typesReal, oneArgDims),
    Funb(expm1, [(-88.72284f0, 88.72284f0)], typesReal, oneArgDims),
    Funb(log, [(0, 1000)], typesReal, oneArgDims),
    # Funb(log10, [(0, 1000)], typesReal, oneArgDims), # faulty
    Funb(abs, [(-1000, 1000)], typesReal, oneArgDims),
    Funb(abs2, [(-1000, 1000)], typesReal, oneArgDims),
    Funb(ceil, [(-1000, 1000)], typesReal, oneArgDims),
    Funb(floor, [(-1000, 1000)], typesReal, oneArgDims),
    Funb(round, [(-1000, 1000)], typesReal, oneArgDims),
    Funb(trunc, [(-1000, 1000)], typesReal, oneArgDims),
    Funb(cis, [(-1000, 1000)], typesReal, oneArgDims),
    Funb(SpecialFunctions.erf, [(-3.8325067f0, 3.8325067f0)], typesReal, oneArgDims),
    Funb(SpecialFunctions.erfc, [(-3.7439213f0, 10.019834f0)], typesReal, oneArgDims),
    Funb(SpecialFunctions.erfinv, [(-1, 1)], typesReal, oneArgDims),
    Funb(SpecialFunctions.erfcinv, [(0, 2)], typesReal, oneArgDims),
    # Funb(SpecialFunctions.lgamma, [(0, 1000)], typesReal, oneArgDims),
    Funb(SpecialFunctions.gamma, [(0, 36)], typesReal, oneArgDims),

    # twoArgDims
    Funb(atan, [(-1, 1), (-1, 1)], typesReal, twoArgDims),
    Funb(hypot, [(-1000, 1000), (-1000, 1000)], typesReal, twoArgDims),
    # Funb(/, [(-1000, 1000), (-1000, 1000)], typesReal, twoArgDims),
    # Funb(^, [(0, 100), (-5, 20)], typesReal, twoArgDims),
])

configsRealIVM = FunbArray([

    # oneArgDims
    Funb(IVM.acos, [(-1, 1)], typesReal, oneArgDims),
    Funb(IVM.asin, [(-1, 1)], typesReal, oneArgDims),
    Funb(IVM.atan, [(-50, 50)], typesReal, oneArgDims),
    Funb(IVM.cos, [(-1000, 1000)], typesReal, oneArgDims),
    Funb(IVM.sin, [(-1000, 1000)], typesReal, oneArgDims),
    Funb(IVM.tan, [(-1000, 1000)], typesReal, oneArgDims),
    Funb(IVM.acosh, [(1, 1000)], typesReal, oneArgDims),
    Funb(IVM.asinh, [(-1000, 1000)], typesReal, oneArgDims),
    Funb(IVM.atanh, [(-1, 1)], typesReal, oneArgDims),
    Funb(IVM.cosh, [(0, 89.415985f0)], typesReal, oneArgDims),
    Funb(IVM.sinh, [(-89.415985f0, 89.415985f0)], typesReal, oneArgDims),
    Funb(IVM.tanh, [(-8.66434f0, 8.66434f0)], typesReal, oneArgDims),
    Funb(IVM.cbrt, [(-1000, 1000)], typesReal, oneArgDims),
    Funb(IVM.sqrt, [(0, 1000)], typesReal, oneArgDims),
    Funb(IVM.exp, [(-88.72284f0, 88.72284f0)], typesReal, oneArgDims),
    Funb(IVM.expm1, [(-88.72284f0, 88.72284f0)], typesReal, oneArgDims),
    Funb(IVM.log, [(0, 1000)], typesReal, oneArgDims),
    # Funb(IVM.log10, [(0, 1000)], typesReal, oneArgDims), # faulty
    Funb(IVM.abs, [(-1000, 1000)], typesReal, oneArgDims),
    Funb(IVM.abs2, [(-1000, 1000)], typesReal, oneArgDims),
    Funb(IVM.ceil, [(-1000, 1000)], typesReal, oneArgDims),
    Funb(IVM.floor, [(-1000, 1000)], typesReal, oneArgDims),
    Funb(IVM.round, [(-1000, 1000)], typesReal, oneArgDims),
    Funb(IVM.trunc, [(-1000, 1000)], typesReal, oneArgDims),
    Funb(IVM.cis, [(-1000, 1000)], typesReal, oneArgDims),
    Funb(IVM.erf, [(-3.8325067f0, 3.8325067f0)], typesReal, oneArgDims),
    Funb(IVM.erfc, [(-3.7439213f0, 10.019834f0)], typesReal, oneArgDims),
    Funb(IVM.erfinv, [(-1, 1)], typesReal, oneArgDims),
    Funb(IVM.erfcinv, [(0, 2)], typesReal, oneArgDims),
    # Funb(IVM.lgamma, [(0, 1000)], typesReal, oneArgDims), # faulty
    Funb(IVM.gamma, [(0, 36)], typesReal, oneArgDims),

    # twoArgDims
    Funb(IVM.atan, [(-1, 1), (-1, 1)], typesReal, twoArgDims),
    Funb(IVM.hypot, [(-1000, 1000), (-1000, 1000)], typesReal, twoArgDims),
    # Funb(IVM.(/), [(-1000, 1000), (-1000, 1000)], typesReal, twoArgDims),
    # Funb(IVM.(^), [(0, 100), (-5, 20)], typesReal, twoArgDims),
])

################################################################
# Complex Functions

typesComplex=[Complex{Float32}, Complex{Float64}]

configsComplexBase = FunbArray([
    # oneArgDims
    Funb(acos, [(-1, 1)], typesComplex, oneArgDims),
    Funb(asin, [(-1, 1)], typesComplex, oneArgDims),
    # Funb(atan, [(-50, 50)], typesComplex, oneArgDims),
    # Funb(cos, [(-10, 10)], typesComplex, oneArgDims),
    # Funb(sin, [(-10, 10)], typesComplex, oneArgDims),
    # Funb(tan, [(-10, 10)], typesComplex, oneArgDims),
    Funb(acosh, [(1, 1000)], typesComplex, oneArgDims),
    Funb(asinh, [(-1000, 1000)], typesComplex, oneArgDims),
    # Funb(atanh, [(-1, 1)], typesComplex, oneArgDims),
    # Funb(cosh, [(0, 89.415985f0)], typesComplex, oneArgDims),
    # Funb(sinh, [(-89.415985f0, 89.415985f0)], typesComplex, oneArgDims),
    # Funb(tanh, [(-8.66434f0, 8.66434f0)], typesComplex, oneArgDims),
    Funb(sqrt, [(0, 1000)], typesComplex, oneArgDims),
    Funb(exp, [(-88.72284f0, 88.72284f0)], typesComplex, oneArgDims),
    Funb(log, [(0, 1000)], typesComplex, oneArgDims),
    # Funb(log10, [(0, 1000)], typesComplex, oneArgDims),
    Funb(abs, [(-1000, 1000)], typesComplex, oneArgDims),
    Funb(angle, [(-1000, 1000)], typesComplex, oneArgDims),
    Funb(conj, [(-1000, 1000)], typesComplex, oneArgDims),
    # Funb(cis, [(-1000, 1000)], typesComplex, oneArgDims),

    # twoArgDims
    # Funb(/, [(-1000, 1000), (-1000, 1000)], typesComplex, twoArgDims),
    # Funb(^, [(0, 100), (-2, 10)], typesComplex, twoArgDims),

    # missing?
    ## oneArgDims

    # Funb(cbrt, [(-1000, 1000)], typesComplex, oneArgDims),
    # Funb(expm1, [(-88.72284f0, 88.72284f0)], typesComplex, oneArgDims),
    # Funb(abs2, [(-1000, 1000)], typesComplex, oneArgDims),
    # Funb(ceil, [(-1000, 1000)], typesComplex, oneArgDims),
    # Funb(SpecialFunctions.erf, [(-3.8325067f0, 3.8325067f0)], typesComplex, oneArgDims),
    # Funb(SpecialFunctions.erfc, [(-3.7439213f0, 10.019834f0)], typesComplex, oneArgDims),
    # Funb(SpecialFunctions.erfinv, [(-1, 1)], typesComplex, oneArgDims),
    # Funb(SpecialFunctions.erfcinv, [(0, 2)], typesComplex, oneArgDims),
    # Funb(SpecialFunctions.lgamma, [(0, 1000)], typesComplex, oneArgDims),
    # Funb(SpecialFunctions.gamma, [(0, 36)], typesComplex, oneArgDims),
    # Funb(floor, [(-1000, 1000)], typesComplex, oneArgDims),
    # Funb(round, [(-1000, 1000)], typesComplex, oneArgDims),
    # Funb(trunc, [(-1000, 1000)], typesComplex, oneArgDims),

    ## twoArgDims
    # Funb(atan, [(-1, 1), (-1, 1)], typesComplex, twoArgDims),
    # Funb(hypot, [(-1000, 1000), (-1000, 1000)], typesComplex, twoArgDims),
])

configsComplexIVM = FunbArray([
    # oneArgDims
    Funb(IVM.acos, [(-1, 1)], typesComplex, oneArgDims),
    Funb(IVM.asin, [(-1, 1)], typesComplex, oneArgDims),
    # Funb(IVM.atan, [(-50, 50)], typesComplex, oneArgDims),
    # Funb(IVM.cos, [(-10, 10)], typesComplex, oneArgDims),
    # Funb(IVM.sin, [(-10, 10)], typesComplex, oneArgDims),
    # Funb(IVM.tan, [(-10, 10)], typesComplex, oneArgDims),
    Funb(IVM.acosh, [(1, 1000)], typesComplex, oneArgDims),
    Funb(IVM.asinh, [(-1000, 1000)], typesComplex, oneArgDims),
    # Funb(IVM.atanh, [(-1, 1)], typesComplex, oneArgDims),
    # Funb(IVM.cosh, [(0, 89.415985f0)], typesComplex, oneArgDims),
    # Funb(IVM.sinh, [(-89.415985f0, 89.415985f0)], typesComplex, oneArgDims),
    # Funb(IVM.tanh, [(-8.66434f0, 8.66434f0)], typesComplex, oneArgDims),
    Funb(IVM.sqrt, [(0, 1000)], typesComplex, oneArgDims),
    Funb(IVM.exp, [(-88.72284f0, 88.72284f0)], typesComplex, oneArgDims),
    Funb(IVM.log, [(0, 1000)], typesComplex, oneArgDims),
    # Funb(IVM.log10, [(0, 1000)], typesComplex, oneArgDims),
    Funb(IVM.abs, [(-1000, 1000)], typesComplex, oneArgDims),
    Funb(IVM.angle, [(-1000, 1000)], typesComplex, oneArgDims),
    Funb(IVM.conj, [(-1000, 1000)], typesComplex, oneArgDims),
    # Funb(IVM.cis, [(-1000, 1000)], typesComplex, oneArgDims),

    # twoArgDims
    # Funb(IVM.(/), [(-1000, 1000), (-1000, 1000)], typesComplex, twoArgDims),
    # Funb(IVM.(^), [(0, 100), (-2, 10)], typesComplex, twoArgDims),

    # missing?
    ## oneArgDims

    # Funb(IVM.cbrt, [(-1000, 1000)], typesComplex, oneArgDims),
    # Funb(IVM.expm1, [(-88.72284f0, 88.72284f0)], typesComplex, oneArgDims),
    # Funb(IVM.abs2, [(-1000, 1000)], typesComplex, oneArgDims),
    # Funb(IVM.ceil, [(-1000, 1000)], typesComplex, oneArgDims),
    # Funb(IVM.erf, [(-3.8325067f0, 3.8325067f0)], typesComplex, oneArgDims),
    # Funb(IVM.erfc, [(-3.7439213f0, 10.019834f0)], typesComplex, oneArgDims),
    # Funb(IVM.erfinv, [(-1, 1)], typesComplex, oneArgDims),
    # Funb(IVM.erfcinv, [(0, 2)], typesComplex, oneArgDims),
    # Funb(IVM.lgamma, [(0, 1000)], typesComplex, oneArgDims),
    # Funb(IVM.gamma, [(0, 36)], typesComplex, oneArgDims),
    # Funb(IVM.floor, [(-1000, 1000)], typesComplex, oneArgDims),
    # Funb(IVM.round, [(-1000, 1000)], typesComplex, oneArgDims),
    # Funb(IVM.trunc, [(-1000, 1000)], typesComplex, oneArgDims),

    ## twoArgDims
    # Funb(IVM.atan, [(-1, 1), (-1, 1)], typesComplex, twoArgDims),
    # Funb(IVM.hypot, [(-1000, 1000), (-1000, 1000)], typesComplex, twoArgDims),
])

################################################################
# Performing Benchmarks
# Real
println("\nBenchmarking Base Real Functions\n")
benchmark!(configsRealBase)

println("\nBenchmarking IntelVectorMath Real Functions\n")
benchmark!(configsRealIVM)

cd(joinpath(rootPath,"Real"))

bar(configsRealBase => configsRealIVM, uniqueType = true, dimAnnotation = false, uniqueDim = true, configName =  "Base" => "IntelVectorMath")

dimplot([configsRealBase,configsRealIVM],["Base", "IntelVectorMath"])
################################################################
# Complex
println("\nBenchmarking Base Complex Functions\n")
benchmark!(configsComplexBase)

println("\nBenchmarking IntelVectorMath Complex Functions\n")
benchmark!(configsComplexIVM)

cd(joinpath(rootPath,"Complex"))

bar(configsComplexBase => configsComplexIVM, uniqueType = true, dimAnnotation = false, uniqueDim = true, configName = "Base" => "IntelVectorMath")

dimplot([configsComplexBase,configsComplexIVM],["Base", "IntelVectorMath"])


text:
# Autogenerated wrapper script for obstack_jll for x86_64-linux-musl
export libobstack

JLLWrappers.@generate_wrapper_header("obstack")
JLLWrappers.@declare_library_product(libobstack, "libobstack.so.1")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libobstack,
        "lib/libobstack.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module FEFunctionsTests

using FillArrays
using Test
using Gridap.Arrays
using Gridap.TensorValues
using Gridap.ReferenceFEs
using Gridap.Geometry
using Gridap.FESpaces
using Gridap.CellData
using LinearAlgebra

order = 1
domain =(0,1,0,1)
partition = (3,3)
model = CartesianDiscreteModel(domain,partition)
grid_topology = get_grid_topology(model)
polytopes = get_polytopes(grid_topology)
reffes = [LagrangianRefFE(Float64,p,order) for p in polytopes]
dirichlet_tags = ["tag_1","tag_6"]
V = GradConformingFESpace(reffes,model,dirichlet_tags)

vh = FEFunction(V,rand(num_free_dofs(V)))
@test is_a_fe_function(vh)
test_fe_function(vh)

cellids = [1,3,5,2]
cell_vals = get_cell_values(vh,cellids)
@test cell_vals == reindex(get_cell_values(vh),cellids)

cellidsL = cellids
cellidsR = [2,4,3,1]
cellidsS = SkeletonPair(cellidsL,cellidsR)
cell_vals = get_cell_values(vh,cellidsS)
@test isa(cell_vals[1],BlockArrayCoo)

end # module