text:
# in this overload we will work on some meta/debug information management
function CC.typeinf(analyzer::AbstractAnalyzer, frame::InferenceState)
    linfo = frame.linfo

    #= logging start =#
    local sec::Float64, depth::Int
    logger_activated = !isnothing(JETLogger(analyzer).inference_logger)
    if logger_activated
        sec = time()
        with_inference_logger(analyzer, ==(DEBUG_LOGGER_LEVEL)) do io
            depth = get_depth(analyzer)

            print_rails(io, depth)
            printstyled(io, "┌ @ "; color = RAIL_COLORS[(depth+1)%N_RAILS+1])
            print(io, linfo)
            if is_constant_propagated(frame)
                printstyled(io, " (constant prop': ", frame.result.argtypes, ')'; color = NOERROR_COLOR)
            end
            file, line = get_file_line(linfo)
            print(io, ' ', file, ':', line)
            println(io)
            set_depth!(analyzer, get_depth(analyzer) + 1) # manipulate this only in debug mode
        end
    end
    #= logging end =#

    prev_frame = get_current_frame(analyzer)
    set_current_frame!(analyzer, frame)

    ret = @invoke typeinf(analyzer::AbstractInterpreter, frame::InferenceState)

    empty!(get_throw_locs(analyzer))
    set_current_frame!(analyzer, prev_frame)

    #= logging start =#
    if logger_activated
        sec = round(time() - sec; digits = 3)
        with_inference_logger(analyzer, ==(INFO_LOGGER_LEVEL)) do io
            println(io, "inference on $linfo finished in $sec sec")
        end
        with_inference_logger(analyzer, ==(DEBUG_LOGGER_LEVEL)) do io
            print_rails(io, depth)
            printstyled(io, "└─→ "; color = RAIL_COLORS[(depth+1)%N_RAILS+1])
            result = get_result(frame)
            isa(result, InferenceState) || printstyled(io, result; color = TYPE_ANNOTATION_COLOR)
            println(io, " (",
                        join(filter(!isnothing, (
                             linfo,
                             ret ? nothing : "in cycle",
                             "$(length(get_reports(analyzer))) reports",
                             "$sec sec"
                             )), ", "),
                        ')')
            set_depth!(analyzer, get_depth(analyzer) - 1) # manipulate this only in debug mode
        end
    end
    #= logging end =#

    return ret
end

# TODO: disable optimization for better performance, only do necessary analysis work by ourselves

# in this overload we can work on `frame.src::CodeInfo` (and also `frame::InferenceState`)
# where type inference (and also optimization if applied) already ran on
function CC._typeinf(analyzer::AbstractAnalyzer, frame::InferenceState)
    linfo = frame.linfo
    parent = frame.parent
    isentry = parent === nothing
    iscp = is_constant_propagated(frame)

    # some methods like `getproperty` can't propagate accurate types without actual values,
    # and constant prop' plays a somewhat critical role in those cases by overwriteing the
    # previous non-constant inference result (under the current design constant prop' always
    # happens after inference with non-constant abstract elements)
    # JET also needs that in order to reduce false positive reports, and here we will
    # throw-away previously-collected error reports that are "lineage" of this frame,
    # when it is being re-inferred with constants
    # NOTE `frame.linfo` is the exactly same object as that of the previous non-constant inference
    # IDEA we may still want to keep some "serious" error reports like `GlobalUndefVarErrorReport`
    # even when constant prop' reveals it never happens given the current constant arguments
    iscp && !isentry && filter!(!is_from_same_frame(parent.linfo, linfo), get_reports(analyzer))

    reports_before             = Set(get_reports(analyzer))
    uncaught_exceptions_before = Set(get_uncaught_exceptions(analyzer))

    ret = @invoke _typeinf(analyzer::AbstractInterpreter, frame::InferenceState)

    # report pass for (local) undef var error
    ReportPass(analyzer)(LocalUndefVarErrorReport, analyzer, frame, frame.src.code)

    # XXX this is a dirty fix for performance problem, we need more "proper" fix
    # https://github.com/aviatesk/JET.jl/issues/75
    unique!(report_identity_key, get_reports(analyzer))

    # report pass for uncaught `throw` calls
    ReportPass(analyzer)(UncaughtExceptionReport, analyzer, frame, frame.src.code)

    reports_after = Set(get_reports(analyzer))
    uncaught_exceptions_after = Set(get_uncaught_exceptions(analyzer))

    # compute JET analysis results that should be cached for this linfo
    this_caches = union!(setdiff!(reports_after, reports_before),
                         setdiff!(uncaught_exceptions_after, uncaught_exceptions_before))

    if !isempty(this_caches)
        if iscp
            result = frame.result
            argtypes = result.argtypes
            cache = get_cache(analyzer)
            @static JET_DEV_MODE && @assert jet_cache_lookup(linfo, argtypes, cache) === nothing "invalid local caching $linfo, $argtypes"
            local_cache = InferenceErrorReportCache[]
            for report in this_caches
                cache_report!(local_cache, report)
                # TODO make this holds when the `analyzer` hooks into `finish` or `optimize`
                # more generally, handle cycles correctly
                @static JET_DEV_MODE && if isa(analyzer, JETAnalyzer)
                    actual, expected = first(report.vst).linfo, linfo
                    @assert actual === expected "invalid local caching detected, expected $expected but got $actual"
                end
            end
            # branching on https://github.com/JuliaLang/julia/pull/39972
            given_argtypes, overridden_by_const = @static if VERSION ≥ v"1.7.0-DEV.705"
                def = result.linfo.def
                va_overwride = isa(def, Method) && def.is_for_opaque_closure
                matching_cache_argtypes(linfo, argtypes, va_overwride)
            else
                matching_cache_argtypes(linfo, argtypes)
            end
            push!(cache, AnalysisResult(linfo, given_argtypes, overridden_by_const, local_cache))
        elseif frame.cached # only cache when `NativeInterpreter` does
            cache = jet_report_cache(analyzer)
            @static JET_DEV_MODE && @assert !haskey(cache, linfo) || isentry "invalid global caching $linfo"
            global_cache = InferenceErrorReportCache[]
            for report in this_caches
                cache_report!(global_cache, report)
                # TODO make this holds when the `analyzer` hooks into `finish` or `optimize`
                # more generally, handle cycles correctly
                @static JET_DEV_MODE && if isa(analyzer, JETAnalyzer)
                    actual, expected = first(report.vst).linfo, linfo
                    @assert actual === expected "invalid global caching detected, expected $expected but got $actual"
                end
            end
            cache[linfo] = global_cache
        end
    end

    set_to_be_updated!(analyzer, this_caches)

    if !iscp && !isentry
        # refinement for this `linfo` may change analysis result for parent frame
        # XXX: is this okay from performance perspective ?
        add_backedge!(linfo, parent)
    end

    return ret
end

"""
    is_from_same_frame(parent_linfo::MethodInstance, current_linfo::MethodInstance) ->
        (report::InferenceErrorReport) -> Bool

Returns a function that checks if a given `InferenceErrorReport` is generated from `current_linfo`.
It also checks `current_linfo` is a "lineage" of `parent_linfo` (i.e. entered from it).

This function is supposed to be used to filter out reports collected from analysis on `current_linfo`
  without using constants when entering into the constant analysis. As such, this function
  assumes that when a report should be filtered out, the first elment of its virtual stack
  frame `st` is for `parent_linfo` and the second element of that is for `current_linfo`.

Example: Assume `linfo2` will produce a report for some reason.
```
entry
└─ linfo1
   ├─ linfo2 (report1: linfo2)
   ├─ linfo3 (report1: linfo1->linfo2, report2: linfo3->linfo2)
   │  └─ linfo2 (report1: linfo1->linfo2, report2: linfo2)
   └─ linfo3′ (report1: linfo1->linfo2, ~~report2: linfo1->linfo3->linfo2~~)
```
In the example analysis above, `report2` will be filtered out on re-entering into `linfo3′`
  (i.e. we're analyzing `linfo3` with constants argument), because
  `is_from_same_frame(linfo1, linfo3)(report2)` returns `true`.
Note that `report1` is still kept there because of the lineage check, i.e.
  `is_from_same_frame(linfo1, linfo3)(report1)` returns `false`.
"""
function is_from_same_frame(parent_linfo::MethodInstance,
                            current_linfo::MethodInstance,
                            )
    function (report::InferenceErrorReport)
        @inbounds begin
            vst = report.vst
            length(vst) > 1 || return false
            vst[1].linfo === parent_linfo || return false
            return vst[2].linfo === current_linfo
        end
    end
end

is_unreachable(@nospecialize(x)) = isa(x, ReturnNode) && !isdefined(x, :val)

@reportdef struct LocalUndefVarErrorReport <: InferenceErrorReport
    name::Symbol
end
get_msg(T::Type{LocalUndefVarErrorReport}, analyzer::AbstractAnalyzer, state, name::Symbol) =
    return "local variable $(name) is not defined"

# these report passes use `:throw_undef_if_not` and `:(unreachable)` introduced by the native
# optimization pass, and thus supposed to only work on post-optimization code
(::SoundPass)(::Type{LocalUndefVarErrorReport}, analyzer::AbstractAnalyzer, frame::InferenceState, stmts::Vector{Any}) =
    report_undefined_local_slots!(analyzer, frame, stmts, false)
(::BasicPass)(::Type{LocalUndefVarErrorReport}, analyzer::AbstractAnalyzer, frame::InferenceState, stmts::Vector{Any}) =
    report_undefined_local_slots!(analyzer, frame, stmts, true)

function report_undefined_local_slots!(analyzer::AbstractAnalyzer, frame::InferenceState, stmts::Vector{Any}, unsound::Bool)
    for (idx, stmt) in enumerate(stmts)
        if isa(stmt, Expr) && stmt.head === :throw_undef_if_not
            sym = stmt.args[1]::Symbol

            # slots in toplevel frame may be a abstract global slot
            istoplevel(analyzer, frame) && is_global_slot(analyzer, sym) && continue

            if unsound
                next_idx = idx + 1
                if checkbounds(Bool, stmts, next_idx) && is_unreachable(@inbounds stmts[next_idx])
                    # the optimization so far has found this statement is never "reachable";
                    # JET reports it since it will invoke undef var error at runtime, or will just
                    # be dead code otherwise
                    add_new_report!(LocalUndefVarErrorReport(analyzer, (frame, idx), sym), analyzer)
                else
                    # by excluding this pass, this analysis accepts some false negatives and
                    # some undefined variable error may happen in actual execution (thus unsound)
                end
            else
                add_new_report!(LocalUndefVarErrorReport(analyzer, (frame, idx), sym), analyzer)
            end
        end
    end
end

@withmixedhash struct ReportIdentityKey
    T::Type{<:InferenceErrorReport}
    sig::Vector{Any}
    # entry_frame::VirtualFrame
    error_frame::VirtualFrame
end

report_identity_key(report::T) where {T<:InferenceErrorReport} =
    ReportIdentityKey(T, report.sig, #=first(report.vst),=# last(report.vst))

"""
    UncaughtExceptionReport <: InferenceErrorReport

Represents general `throw` calls traced during inference.
This is reported only when it's not caught by control flow.
"""
@reportdef struct UncaughtExceptionReport <: InferenceErrorReport
    throw_calls::Vector{Tuple{Int,Expr}} # (pc, call)
end
function UncaughtExceptionReport(analyzer::AbstractAnalyzer, sv::InferenceState, throw_calls::Vector{Tuple{Int,Expr}})
    vf = get_virtual_frame(analyzer, sv.linfo)
    msg = length(throw_calls) == 1 ? "may throw" : "may throw either of"
    sig = Any[]
    ncalls = length(throw_calls)
    for (i, (pc, call)) in enumerate(throw_calls)
        call_sig = _get_sig(analyzer, (sv, pc), call)
        append!(sig, call_sig)
        i ≠ ncalls && push!(sig, ", ")
    end
    return UncaughtExceptionReport([vf], msg, sig, throw_calls)
end

# this error report is really special, and might not be eligible for possible overloads since:
# - this report does not only report error points but also "clean up" caught error points
# - this pass is tightly bound to that of `SeriousExceptionReport`
function (::SoundBasicPass)(::Type{UncaughtExceptionReport}, analyzer::AbstractAnalyzer, frame::InferenceState, stmts::Vector{Any})
    # report `throw` calls "appropriately"
    if get_result(frame) === Bottom
        # if the return type here is `Bottom` annotated, this _may_ mean there're uncaught
        # `throw` calls
        # XXX: well, it's possible that the `throw` calls within them are all caught but the
        # other critical errors make the return type `Bottom`
        # NOTE: to reduce the false positive `UncaughtExceptionReport`s described above, we count
        # `throw` calls here after optimization, since it may have eliminated "unreachable"
        # `throw` calls
        codelocs    = frame.src.codelocs
        linetable   = frame.src.linetable::Vector
        throw_locs  = get_throw_locs(analyzer)
        throw_calls = Tuple{Int,Expr}[]
        for (pc, stmt) in enumerate(stmts)
            is_throw_call_expr(analyzer, frame, stmt) || continue
            # if this `throw` is already reported, don't duplciate
            linetable[codelocs[pc]]::LineInfoNode in throw_locs && continue
            push!(throw_calls, (pc, stmt))
        end
        if !isempty(throw_calls)
            add_new_report!(UncaughtExceptionReport(analyzer, frame, throw_calls), analyzer)
        end
    else
        # the non-`Bottom` result here may mean `throw` calls from the children frames
        # (if exists) are caught and not propagated here;
        # we don't want to cache `UncaughtExceptionReport`s for those calls for this frame
        # and its parents, so just filter them away
        empty!(get_uncaught_exceptions(analyzer))
    end
end

# basically same as `is_throw_call`, but also toplevel module handling added
function is_throw_call_expr(analyzer::AbstractAnalyzer, frame::InferenceState, @nospecialize(e))
    if isa(e, Expr)
        if e.head === :call
            f = e.args[1]
            if istoplevel(analyzer, frame) && isa(f, Symbol)
                f = GlobalRef(get_toplevelmod(analyzer), f)
            end
            if isa(f, GlobalRef)
                ff = CC.abstract_eval_global(f.mod, f.name)
                if isa(ff, Const) && ff.val === Core.throw
                    return true
                end
            end
        end
    end
    return false
end


text:
using Test
using Optim
using OptimTestProblems
using OptimTestProblems.MultivariateProblems
const MVP = MultivariateProblems

import PositiveFactorizations: Positive, cholesky # for the IPNewton tests
using Random

import LineSearches
import NLSolversBase
import NLSolversBase: clear!
import LinearAlgebra: norm, diag, I, Diagonal, dot, eigen, issymmetric, mul!
import SparseArrays: normalize!, spdiagm

debug_printing = false

special_tests = [
    "bigfloat/initial_convergence",
]
special_tests = map(s->"./special/"*s*".jl", special_tests)

general_tests = [
    "api",
    "callables",
    "callbacks",
    "convergence",
    "default_solvers",
    "deprecate",
    "initial_convergence",
    "objective_types",
    "Optim",
    "optimize",
    "type_stability",
#    "types",
    "counter",
    "maximize",
]
general_tests = map(s->"./general/"*s*".jl", general_tests)

univariate_tests = [
    # optimize
    "optimize/interface",
    "optimize/optimize",
    # solvers
    "solvers/golden_section",
    "solvers/brent",
    # "initial_convergence",
    "dual",
]
univariate_tests = map(s->"./univariate/"*s*".jl", univariate_tests)

multivariate_tests = [
    ## optimize
    "optimize/interface",
    "optimize/optimize",
    "optimize/inplace",
    ## solvers
    ## constrained
    "solvers/constrained/fminbox",
    "solvers/constrained/ipnewton/constraints",
    "solvers/constrained/ipnewton/counter",
    "solvers/constrained/ipnewton/ipnewton_unconstrained",
    "solvers/constrained/samin",
    ## first order
    "solvers/first_order/accelerated_gradient_descent",
    "solvers/first_order/bfgs",
    "solvers/first_order/cg",
    "solvers/first_order/gradient_descent",
    "solvers/first_order/l_bfgs",
    "solvers/first_order/momentum_gradient_descent",
    "solvers/first_order/ngmres",
    ## second order
    "solvers/second_order/newton",
    "solvers/second_order/newton_trust_region",
    "solvers/second_order/krylov_trust_region",
    ## zeroth order
    "solvers/zeroth_order/grid_search",
    "solvers/zeroth_order/nelder_mead",
    "solvers/zeroth_order/particle_swarm",
    "solvers/zeroth_order/simulated_annealing",
    ## other
    "array",
    "extrapolate",
    "lsthrow",
    "precon",
    "manifolds",
    "complex",
    "fdtime",
    "arbitrary_precision",
    "successive_f_tol",
    "f_increase",
]
multivariate_tests = map(s->"./multivariate/"*s*".jl", multivariate_tests)

input_tuple(method, prob) = ((MVP.objective(prob),),)
input_tuple(method::Optim.FirstOrderOptimizer, prob) = ((MVP.objective(prob),), (MVP.objective(prob), MVP.gradient(prob)))
input_tuple(method::Optim.SecondOrderOptimizer, prob) = ((MVP.objective(prob),), (MVP.objective(prob), MVP.gradient(prob)), (MVP.objective(prob), MVP.gradient(prob), MVP.hessian(prob)))

function run_optim_tests(method; convergence_exceptions = (),
                         minimizer_exceptions = (),
                         minimum_exceptions = (),
                         f_increase_exceptions = (),
                         iteration_exceptions = (),
                         skip = (),
                         show_name = false,
                         show_trace = false,
                         show_res = false,
                         show_itcalls = false)
    # Loop over unconstrained problems
    for (name, prob) in MultivariateProblems.UnconstrainedProblems.examples
        if !isfinite(prob.minimum) || !any(isfinite, prob.solutions)
            debug_printing && println("$name has no registered minimum/minimizer. Skipping ...")
            continue
        end

        show_name && printstyled("Problem: ", name, "\n", color=:green)
        # Look for name in the first elements of the iteration_exceptions tuples
        iter_id = findall(n->n[1] == name, iteration_exceptions)
        # If name wasn't found, use default 1000 iterations, else use provided number
        iters = length(iter_id) == 0 ? 1000 : iteration_exceptions[iter_id[1]][2]
        # Construct options
        allow_f_increases = (name in f_increase_exceptions)
        dopts = Optim.default_options(method)
        if haskey(dopts, :allow_f_increases)
            allow_f_increases = allow_f_increases || dopts[:allow_f_increases]
            delete!(dopts, :allow_f_increases)
        end
        options = Optim.Options(allow_f_increases = allow_f_increases,
                                iterations = iters, show_trace = show_trace;
                                dopts...)

        # Use finite difference if it is not differentiable enough
        if  !(name in skip)
            for (i, input) in enumerate(input_tuple(method, prob))
                if (!prob.isdifferentiable && i > 1) || (!prob.istwicedifferentiable && i > 2)
                    continue
                end

                # Loop over appropriate input combinations of f, g!, and h!
                results = Optim.optimize(input..., prob.initial_x, method, options)
                @test isa(summary(results), String)
                show_res && println(results)
                show_itcalls && printstyled("Iterations: $(Optim.iterations(results))\n", color=:red)
                show_itcalls && printstyled("f-calls: $(Optim.f_calls(results))\n", color=:red)
                show_itcalls && printstyled("g-calls: $(Optim.g_calls(results))\n", color=:red)
                show_itcalls && printstyled("h-calls: $(Optim.h_calls(results))\n", color=:red)
                if !((name, i) in convergence_exceptions)
                    @test Optim.converged(results)
                    # Print on error, easier to debug CI
                    if !(Optim.converged(results))
                        printstyled(name, " did not converge with i = ", i, "\n", color=:red)
                        printstyled(results, "\n", color=:red)
                    end
                end
                if !((name, i) in minimum_exceptions)
                    @test Optim.minimum(results) < prob.minimum + sqrt(eps(typeof(prob.minimum)))
                end
                if !((name, i) in minimizer_exceptions)
                    @test norm(Optim.minimizer(results) - prob.solutions) < 1e-2
                end
            end
        else
            debug_printing && printstyled("Skipping $name\n", color=:blue)
        end
    end
end

function run_optim_tests_constrained(method; convergence_exceptions = (),
                         minimizer_exceptions = (),
                         minimum_exceptions = (),
                         f_increase_exceptions = (),
                         iteration_exceptions = (),
                         skip = (),
                         show_name = false,
                         show_trace = false,
                         show_res = false,
                         show_itcalls = false)
    # TODO: Update with constraint problems too?
    # Loop over unconstrained problems
    for (name, prob) in MVP.UnconstrainedProblems.examples
        if !isfinite(prob.minimum) || !any(isfinite, prob.solutions)
            debug_printing && println("$name has no registered minimum/minimizer. Skipping ...")
            continue
        end
        show_name && printstyled("Problem: ", name, "\n", color=:green)
        # Look for name in the first elements of the iteration_exceptions tuples
        iter_id = findall(n->n[1] == name, iteration_exceptions)
        # If name wasn't found, use default 1000 iterations, else use provided number
        iters = length(iter_id) == 0 ? 1000 : iteration_exceptions[iter_id[1]][2]
        # Construct options
        allow_f_increases = (name in f_increase_exceptions)
        options = Optim.Options(iterations = iters, show_trace = show_trace; Optim.default_options(method)...)

        # Use finite difference if it is not differentiable enough
        if  !(name in skip) && prob.istwicedifferentiable
            # Loop over appropriate input combinations of f, g!, and h!
            df = TwiceDifferentiable(MVP.objective(prob), MVP.gradient(prob),
                                     MVP.objective_gradient(prob), MVP.hessian(prob), prob.initial_x)
            infvec = fill(Inf, size(prob.initial_x))
            constraints = TwiceDifferentiableConstraints(-infvec, infvec)
            results = optimize(df,constraints,prob.initial_x, method, options)
            @test isa(Optim.summary(results), String)
            show_res && println(results)
            show_itcalls && printstyled("Iterations: $(Optim.iterations(results))\n", color=:red)
            show_itcalls && printstyled("f-calls: $(Optim.f_calls(results))\n", color=:red)
            show_itcalls && printstyled("g-calls: $(Optim.g_calls(results))\n", color=:red)
            show_itcalls && printstyled("h-calls: $(Optim.h_calls(results))\n", color=:red)
            if !(name in convergence_exceptions)
                @test Optim.converged(results)
                # Print on error
                if !(Optim.converged(results))
                    printstyled(name, "did not converge\n", color=:red)
                    printstyled(results, "\n", color=:red)
                end
            end
            if !(name in minimum_exceptions)
                @test Optim.minimum(results) < prob.minimum + sqrt(eps(typeof(prob.minimum)))
            end
            if !(name in minimizer_exceptions)
                @test norm(Optim.minimizer(results) - prob.solutions) < 1e-2
            end
        else
            debug_printing && printstyled("Skipping $name\n", color=:blue)
        end
    end
end

@testset "special" begin
    for my_test in special_tests
        println(my_test)
        @time include(my_test)
    end
end
@testset "general" begin
    for my_test in general_tests
        println(my_test)
        @time include(my_test)
    end
end
@testset "univariate" begin
    for my_test in univariate_tests
        println(my_test)
        @time include(my_test)
    end
end
@testset "multivariate" begin
    for my_test in multivariate_tests
        println(my_test)
        @time include(my_test)
    end
end

println("Literate examples")
@time include("examples.jl")

@testset "show method for options" begin
    o = Optim.Options()
    @test occursin(" = ", sprint(show, o))
end

ext: jl


text:
"""
    NormLLData
        y::Vector
        X::Matrix
        Z::Matrix

Holds data for hierarchical linear model, with response `y`, fixed effects
design matrix `X`, and random effects design matrix `Z`.
"""
struct NormLLData{T}
    y::Vector{T}
    X::Matrix{T}
    Z::Matrix{T}
end
"""
    NormLLData(::NormLLData, grp::Integer)

Outer constructor that returns data subset by random effects membership. Used
when integrating out random effect parameters one at a time using quadrature.
"""
function NormLLData(nld::NormLLData, grp::Integer)
    idx = Z[:, grp] .== 1
    NormLLData(nld.y[idx], nld.X[idx, :], nld.Z[idx, grp:grp])
end

"""
    NormLL
        β
        u
        σobs
        σgrps
        nld::NormLLData

Holds parameter values and data for hierarchical linear model. Also a functor to
return the complete data likelihood. Fixed effects parameters are in the vector
`β`, random effects in `u`, observation standard deviation in σobs, and group
standard deviation σgrp. `nld` holds the associated data object.
"""
struct NormLL{Tβ, Tu, Tσobs, Tσgrp, Tnld<:NormLLData} <: Function
    β::Tβ
    u::Tu
    σobs::Tσobs
    σgrp::Tσgrp
    nld::Tnld
end
# Outer constructor to make changing `u` easy
"""
    NormLL(::NormLL, u)

Outer constructor to make swapping random effects vector easy. Makes both
Laplace approximation and quadrature.
"""
function NormLL(NLL::NormLL, u)
    NormLL(NLL.β, u, NLL.σobs, NLL.σgrp, NLL.nld)
end
"""
Negative log-likelihood function. Full data likelihood.
"""
function (NLL::NormLL)()
    pred = NLL.nld.X * NLL.β + NLL.nld.Z * NLL.u
    datalik = -sum(logpdf.(Normal.(pred, NLL.σobs), NLL.nld.y))
    randlik = -sum(logpdf.(Normal(0.0, NLL.σgrp), NLL.u))
    datalik + randlik
end
"""
Single argument version to make Laplace an quadrature simple.
"""
function (NLL::NormLL)(u)
    NormLL(NLL, u)()
end


text:
__precompile__() # this module is safe to precompile
module PyCallCoolProp

export PropsSI, CoolProp

using PyCall

const CoolProp = PyNULL()
const PropsSI  = PyNULL()

function __init__()
    copy!(CoolProp, pyimport_conda("CoolProp.CoolProp", "CoolProp"))
    copy!(PropsSI, pyimport_conda("CoolProp.CoolProp", "CoolProp").PropsSI)
end

end


text:
@testset "multithreads" begin
    @test Threads.nthreads() > 1

    img = testimage("lighthouse")
    @testset "jpeg_encode" begin
        ref = jpeg_encode(img)

        out = [fill(zero(UInt8), size(ref)) for _ in 1:Threads.nthreads()]
        Threads.@threads for i in 1:Threads.nthreads()
            out[i] = jpeg_encode(img)
        end
        @test all(out .== Ref(ref))
    end

    @testset "jpeg_decode" begin
        out = [similar(img) for _ in 1:Threads.nthreads()]
        tmpdir = mktempdir()

        for i in 1:Threads.nthreads()
            jpeg_encode(joinpath(tmpdir, "$i.jpg"), img)
        end
        ref = jpeg_decode(joinpath(tmpdir, "1.jpg"))
        Threads.@threads for i in 1:Threads.nthreads()
            out[i] = jpeg_decode(joinpath(tmpdir, "$i.jpg"))
        end
        @test all(out .== Ref(ref))
    end
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "umat_binaries"
version = v"0.5.1"

# Collection of sources required to build umat_binaries
sources = [
    "https://raw.githubusercontent.com/KratosMultiphysics/Kratos/80d68b9a7e89d97438fd04726e5198084b7be43e/applications/constitutive_laws_application/custom_external_libraries/umat/mises_umat.f" =>
    "bfb0f46ae0c9cef0cc91c1dbf78943dd7afd24a5a757dd3224ed66cdd71e7ec6",

    "https://raw.githubusercontent.com/KratosMultiphysics/Kratos/80d68b9a7e89d97438fd04726e5198084b7be43e/applications/constitutive_laws_application/custom_external_libraries/umat/ABA_PARAM.INC" =>
    "b7d74a332dda559e06720db8b7f907aedb72f58b8ed957c7c67ba7007a8e93a8",

    "https://raw.githubusercontent.com/KratosMultiphysics/Kratos/80d68b9a7e89d97438fd04726e5198084b7be43e/applications/constitutive_laws_application/custom_external_libraries/umat/xit.f" =>
    "14f74acd4a20ad680b0c354416f302e5f2e078f3a5750df60888b9b41cb887a5",

    #"https://sourceforge.net/code-snapshots/svn/p/pa/parafem/code/parafem-code-r2281-trunk-parafem-src-umats-dp.zip" =>
    #"114d5773e806de0070eee7379131ceae25c3516fb49b8e3fb8f0f9e1926717bf"

    # The above is the original but for some reason sourceforge blogged the trafic
    "https://github.com/TeroFrondelius/umat_binaries_builder/releases/download/v0.3.1/parafem-code-r2281-trunk-parafem-src-umats-dp.zip" =>
    "114d5773e806de0070eee7379131ceae25c3516fb49b8e3fb8f0f9e1926717bf",

    # Gurson model from UMAT.jl/umat_models
    "https://github.com/JuliaFEM/UMAT.jl.git" =>
    "db0046a44b2c8dfd5ca5b3212cbd6fee37f6baa3"

]

# Bash recipe for building across all platforms
prefix_script = """
export name=$name
export version=$version
"""
script = prefix_script * raw"""
cd $WORKSPACE/srcdir
cat >CMakeLists.txt <<EOL
cmake_minimum_required(VERSION 3.5)
project($name)
set(VERSION $version)
enable_language(Fortran)
set(SOURCE_FILES mises_umat.f xit.f)
set(LIBRARY_NAME mises_umat)

add_library(\${LIBRARY_NAME} SHARED \${SOURCE_FILES})
add_library(elastic SHARED elastic.f)
add_library(isotropic_plast_exp SHARED isotropic_plast_exp.f)
add_library(isotropic_plast_imp SHARED isotropic_plast_imp.f)

install(TARGETS mises_umat DESTINATION lib)
install(TARGETS elastic DESTINATION lib)
install(TARGETS isotropic_plast_exp DESTINATION lib)
install(TARGETS isotropic_plast_imp DESTINATION lib)
EOL

export dp=parafem-code-r2281-trunk-parafem-src-umats-dp
cp $dp/elasticity/elastic.f .
cp $dp/plasticity_exp/code_exp.f isotropic_plast_exp.f
cp $dp/plasticity_imp/code_imp.f isotropic_plast_imp.f

sed -i -e '1i\       IMPLICIT REAL*8(A-H,O-Z)' ABA_PARAM.INC

cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=/opt/$target/$target.toolchain
make
make install

################################################################################
# From this forward is UMAT.jl/umat_models building
################################################################################

cd UMAT.jl/umat_models
sed -i 's/CALL ROTSIG(/!CALL ROTSIG(/g' gurson_porous_plasticity.f90

if [[ ${nbits} == 64 ]]; then
    export OB=openblas64_
    sed -i 's/call dgesv(/call dgesv_64(/g' gurson_porous_plasticity.f90
else
    export OB=openblas
fi

cat >CMakeLists.txt <<EOL
cmake_minimum_required(VERSION 3.5)
project($name)
set(VERSION $version)
enable_language(Fortran)
set(CMAKE_Fortran_FLAGS "-fdefault-real-8 -fdefault-integer-8")
add_library(gurson_porous_plasticity SHARED gurson_porous_plasticity.f90)
target_link_libraries(gurson_porous_plasticity $OB)
install(TARGETS gurson_porous_plasticity DESTINATION lib)

add_library(drucker_prager_plasticity SHARED drucker_prager_plasticity.f90 drucker_prager_plasticity_deriv.f90 ../../xit.f)
install(TARGETS drucker_prager_plasticity DESTINATION lib)
EOL

cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=/opt/$target/$target.toolchain
make
make install

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:i686, libc=:glibc),
    Windows(:x86_64),
    MacOS(:x86_64),
    Windows(:i686),
    Linux(:x86_64, libc=:glibc),
    Linux(:aarch64, libc=:glibc),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf),
    Linux(:powerpc64le, libc=:glibc),
    Linux(:i686, libc=:musl),
    Linux(:x86_64, libc=:musl),
    Linux(:aarch64, libc=:musl),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf),
    FreeBSD(:x86_64)
]

platforms = expand_gcc_versions(platforms)

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libmises_umat", :mises_umat),
    LibraryProduct(prefix, "libelastic", :elastic),
    LibraryProduct(prefix, "libisotropic_plast_exp", :isotropic_plast_exp),
    LibraryProduct(prefix, "libisotropic_plast_imp", :isotropic_plast_imp),
    LibraryProduct(prefix, "libgurson_porous_plasticity", :gurson_porous_plasticity),
    LibraryProduct(prefix, "libdrucker_prager_plasticity", :drucker_prager_plasticity)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    "https://github.com/JuliaLinearAlgebra/OpenBLASBuilder/releases/download/v0.3.0-3/build_OpenBLAS.v0.3.0.jl"
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
########################################################################################################################
###                                                     Packages                                                     ###
########################################################################################################################
using Dash
using DashHtmlComponents
using DashCoreComponents
using DashCytoscape
using DashTable
using JSON
using LightGraphs
using LinearAlgebra

########################################################################################################################
###                                                     Auxiliary                                                    ###
########################################################################################################################
#! Inverse Key-Value Function
function reverselookup(d, v)
    for k in keys(d)
        if d[k] == v
            return(k)
        end
    end
    error("Failed to reverse the Key-Value pair.")
end

########################################################################################################################
###                                                     Variables                                                    ###
########################################################################################################################
#! Default Graph Elements
graphElements = [
    Dict("data" => Dict("id" => "v4", "label" => "4"),
    "position" => Dict("x" => 0, "y" => 200)),
    Dict("data" => Dict("id" => "v5", "label" => "5"),
    "position" => Dict("x" => 100, "y" => 200)),
    Dict("data" => Dict("id" => "v1", "label" => "1"),
    "position" => Dict("x" => 0, "y" => 100)),
    Dict("data" => Dict("id" => "v2", "label" => "2"),
    "position" => Dict("x" => 100, "y" => 100)),
    Dict("data" => Dict("id" => "v3", "label" => "3"),
    "position" => Dict("x" => 50, "y" => 25)),
    Dict("data" => Dict("id" => "e12", "source" => "v1", "target" => "v2"),),
    Dict("data" => Dict("id" => "e14", "source" => "v1", "target" => "v4"),),
    Dict("data" => Dict("id" => "e13", "source" => "v1", "target" => "v3"),),
    Dict("data" => Dict("id" => "e23", "source" => "v2", "target" => "v3"),),
    Dict("data" => Dict("id" => "e25", "source" => "v2", "target" => "v5"),),
    Dict("data" => Dict("id" => "e45", "source" => "v4", "target" => "v5"),)
    ]

########################################################################################################################
###                                                  Initialization                                                  ###
########################################################################################################################
app = dash(external_stylesheets = ["https://codepen.io/chriddyp/pen/bWLwgP.css"],
suppress_callback_exceptions = true)

########################################################################################################################
###                                                      Layout                                                      ###
########################################################################################################################
app.layout = html_div(style = Dict("maxWidth" => "100%",)) do
    #! Header
    html_header(style = Dict("maxWidth" => "100%")) do
        html_img(src = "https://i.imgur.com/3Xw9SYD.png",
            style = Dict("display" => "block", "maxWidth" => "100%",
            "height" => "auto", "marginLeft" => "auto", "marginRight" => "auto"))
        end,
    #! Body
    dcc_tabs(id = "Tabs", value = "Tab1") do
        dcc_tab(label="Slide 1", value = "Tab1"),
        dcc_tab(label="Slide 2", value = "Tab2")
    end,
    #! Selected Slide
    html_div(id = "selectedTab"),
    #! Empty Dash DataTable Box
    #? Due to a bug, you need to declare a dash_datatable in the initial layout.
    #? See: https://github.com/plotly/dash/issues/1010
    html_div() do
        html_div(id = "dashtableContent"),
        dcc_location(id = "dashtableLocation", refresh = false),
        html_div(dash_datatable(id = "dtInitial"), style = Dict("display" => "none"))
    end,
    #! Footer
    html_footer() do
        html_br(),
        html_img(src = "https://i.creativecommons.org/l/by/4.0/88x31.png",
        style = Dict("display" => "inline", "align" => "textRight",
        "verticalAlign" => "middle", "paddingRight" => 5)),
        html_div("Copyleft. All rights reversed. XIII Jornada
        de Iniciação Científica da Escola Nacional de Ciência Estatísticas (ENCE).",
        style = Dict("display" => "inline", "textAlign" => "justify"))
    end
end

########################################################################################################################
###                                             Callback: Selected Slide                                             ###
########################################################################################################################
#! Slide Selection
callback!(app,
Output("selectedTab", "children"),
Input("Tabs", "value")) do value
    #! Tab 1
    if value == "Tab1"
        return(
            html_div() do
                #! Column 1
                html_div(className = "five columns", style = Dict("display" => "flex",
                "border" => "2px solid", "height" => "840px", "borderRadius" => "25px", "width" => "30%",
                "marginTop" => "15px", "flexDirection" => "column", "marginLeft" => "5px")) do
                    #! Input: Refresh Info
                    html_div(html_button("Refresh", id = "refreshInfo",
                    style = Dict("fontWeight" => "bold", "fontSize" => "14px")),
                    style = Dict("display" => "block",
                    "padding" => "10px", "marginLeft" => "auto",
                    "marginRight" => "auto", "width" => "auto")),
                    #! Input: Graph Property
                    html_div(style = Dict("display" => "block", "padding" => "10px", "fontWeight" => "bold")) do
                        dcc_dropdown(id = "propertySelection",
                        options = [
                            (label = "Adjacency Matrix", value = "selectAdjMatrix"),
                            (label = "Laplacian Matrix", value = "selectLapMatrix"),
                            (label = "Eigenvalues", value = "selectAlgCon"),
                            (label = "Fiedler's Vector", value = "selectFiedler"),
                            (label = "Reliability", value = "selectReliability")
                            ],
                        value = "selectAdjMatrix")
                    end,
                    #! Output: Interactive DataTable
                    html_div(id = "graphTabOutput",
                    style = Dict("textAlign" => "center", "padding" => "1px",
                    "fontWeight" => "bold", "fontSize" => "20px"))
                end,
                #! Column 2
                html_div(className = "seven columns", style = Dict("display" => "flex",
                "width" => "65%", "flexDirection" => "column")) do
                    #! Input: Graph Manipulation
                    html_div([
                        html_button("Add Vertex",
                        id = "addVtx", n_clicks = nothing,
                        style = Dict("fontWeight" => "bold", "fontSize" => "14px")),
                        html_button("Remove Vertex",
                        id = "removeVtx", n_clicks = nothing,
                        style = Dict("fontWeight" => "bold", "fontSize" => "14px")),
                        html_button("Add Edge",
                        id = "addEdg", n_clicks = nothing,
                        style = Dict("fontWeight" => "bold", "fontSize" => "14px")),
                        html_button("Remove Edge",
                        id = "removeEdg", n_clicks = nothing,
                        style = Dict("fontWeight" => "bold", "fontSize" => "14px"))
                        ],
                    style = Dict("display" => "inlineBlock", "width" => "100%",
                    "marginLeft" => "auto", "fontWeight" => "bold",
                    "marginRight" => "auto", "marginTop" => "20px")),
                    #! Output: Interactive Graph SL2
                    cyto_cytoscape(id = "graphTab",
                    style = Dict("width" => "100%", "height" => "780px"),
                    layout = Dict("name" => "preset",),
                    autoRefreshLayout = false,
                    stylesheet = [
                        Dict("selector" => "node",
                        "style" => Dict("label" => "data(label)",
                        "textHalign" => "center", "textValign" => "center"))
                        ],
                    elements = graphElements)
                end
            end,
            #! Output: Graph Element Info
            html_div(className = "row") do
                html_pre(id = "selectedElement",
                style = Dict("border" => "thin lightgrey solid", "width" => "100%",
                "marginTop" => "15px", "overflowX" => "scroll", "display" => "inlineBlock"))
            end
        )
    #! Tab 2
    elseif value == "Tab2"
        return(
            html_div() do
                dcc_markdown("""
                    Under construction.
                    """,
                    style = Dict("textAlign" => "justify", "fontWeight" => "bold", "fontSize" => "18px",
                    "paddingLeft" => "5px", "paddingRight" => "5px")),
                html_br()
            end
        )
    else
        return(html_div(html_h3("Error: Couldn't generate the tabs")))
    end
end

########################################################################################################################
###                                                Callback: Slide 1                                                 ###
########################################################################################################################
#! Callback: Graph Element Info
callback!(app,
Output("selectedElement", "children"),
[Input("graphTab", "selectedNodeData"),
Input("graphTab", "selectedEdgeData")]) do info_vtx, info_edg
    return([JSON.json(info_vtx), JSON.json(info_edg)])
end

#! Callback: Interactive Graph
callback!(app,
Output("graphTab", "elements"),
[Input("addVtx", "n_clicks"),
Input("removeVtx", "n_clicks"),
Input("addEdg", "n_clicks"),
Input("removeEdg", "n_clicks")],
[State("graphTab", "elements"),
State("graphTab", "selectedNodeData"),
State("graphTab", "selectedEdgeData")]) do mkvtx, rmvtx, mkedg, rmedg, elements, info_vtx, info_edg
    ctx = callback_context().triggered
    if isassigned(ctx, 1)
        #! Graph Manipulation: Add Vertex
        if ctx[1][:prop_id] == "addVtx.n_clicks"
            num_id = [parse(Int, lstrip(x[:data][:id], 'v')) for x ∈ elements if occursin("v", x[:data][:id])]
            n = minimum(setdiff(1:(maximum(num_id) + 1), num_id))
            push!(elements, Dict("data" => Dict("id" => "v$(n)",
            "label" => "$(n)"), "position" => Dict("x" => 50, "y" => 50)))
            return(elements)
        end
        if !isnothing(info_vtx)
            #! Graph Manipulation: Remove Vertex
            if ctx[1][:prop_id] == "removeVtx.n_clicks"
                if length(info_vtx) > 0
                    id_vtx = [vtx[:id] for vtx ∈ info_vtx]
                    rlt_edg = [element[:data][:id] for element ∈ elements if (occursin('e', element[:data][:id]) && (element[:data][:source] ∈ id_vtx || element[:data][:target] ∈ id_vtx))]
                    return(filter(x -> x[:data][:id] ∉ vcat(id_vtx, rlt_edg), elements))
                else
                    println("At least one vertex must be selected in order to remove the element!")
                    return(elements)
                end
            end
            #! Graph Manipulation: Add Edge
            if ctx[1][:prop_id] == "addEdg.n_clicks"
                if length(info_vtx) == 2
                    id_vtx = sort([vtx[:label] for vtx ∈ info_vtx])
                    newedg = "e"*join(id_vtx)
                    if length(filter(x -> x[:data][:id] == newedg, elements)) == 0
                        push!(elements, Dict("data" => Dict("id" => "$(newedg)",
                        "source" => "v$(id_vtx[1])",
                        "target" => "v$(id_vtx[2])"),))
                        return(elements)
                    else
                        println("The selected vertices cannot be neighbours in order to add an edge!")
                        return(elements)
                    end
                else
                    println("Two vertices must be selected in order to add an edge!")
                    return(elements)
                end
            end
        end
        if !isnothing(info_edg)
            #! Graph Manipulation: Remove Edge
            if ctx[1][:prop_id] == "removeEdg.n_clicks"
                if length(info_edg) ≥ 1
                    id_edg = [edg[:id] for edg ∈ info_edg]
                    return(filter(x -> x[:data][:id] ∉ id_edg, elements))
                else
                    println("At least one edge must be selected in order to remove the element!")
                    return(elements)
                end
            end
        end
        return(elements)
    end
    return(elements)
end

#! Callback: Interactive DataTable
callback!(app,
Output("graphTabOutput", "children"),
[Input("refreshInfo", "n_clicks")],
[State("propertySelection", "value"),
State("graphTab", "elements")]) do refresh, property, elements
    ctx = callback_context().triggered
    if isassigned(ctx, 1) && ctx[1][:prop_id] == "refreshInfo.n_clicks"
        vertices = []
        edges = []
        for element ∈ elements
            if occursin("v", element[:data][:id])
                push!(vertices,
                parse(Int, lstrip(element[:data][:id], 'v')))
            else
                push!(edges,
                parse.(Int, lstrip.([element[:data][:source], element[:data][:target]], 'v')))
            end
        end
        sort!(vertices)
        vtxdict = Dict{Int64, Int64}(enumerate(vertices))
        #! Generating graph
        graph = SimpleGraph(length(vertices))
        for edge ∈ edges
        add_edge!(graph, reverselookup(vtxdict, edge[1]), reverselookup(vtxdict, edge[2]))
        end
        #! Graph Property: Adjacency Matrix
        if property == "selectAdjMatrix"
            aux = Matrix(adjacency_matrix(graph))
            return(
                dash_datatable(id = "datatableGraph",
                columns = [Dict("name" => string(vtxdict[i]), "id" => string(i)) for i ∈ 1:size(aux, 2)],
                data = [Dict(string(j) => aux[i,j] for j ∈ 1:size(aux, 2)) for i ∈ 1:size(aux, 1)],
                page_size = 20,
                sort_action = "native",
                style_table = Dict("height" => "auto",
                "width" => "auto", "overflowY" => "auto",
                "overflowX" => "auto"),
                style_cell = Dict("minWidth" => "50px", "fontSize" => "16px"),
                style_header = Dict("fontWeight" => "bold",
                "backgroundColor" => "rgb(230, 230, 230)"))
                )
        #! Graph Property: Laplacian Matrix
        elseif property == "selectLapMatrix"
            aux = Matrix(laplacian_matrix(graph))
            return(
                dash_datatable(id = "datatableGraph",
                columns = [Dict("name" => string(vtxdict[i]), "id" => string(i)) for i ∈ 1:size(aux, 2)],
                data = [Dict(string(j) => aux[i,j] for j ∈ 1:size(aux, 2)) for i ∈ 1:size(aux, 1)],
                page_size = 20,
                sort_action = "native",
                style_table = Dict("height" => "auto",
                "width" => "auto", "overflowY" => "auto",
                "overflowX" => "auto"),
                style_cell = Dict("minWidth" => "50px", "fontSize" => "16px"),
                style_header = Dict("fontWeight" => "bold",
                "backgroundColor" => "rgb(230, 230, 230)"))
                )
        #! Graph Property: Eigenvalues
        elseif property == "selectAlgCon"
            aux = round.(laplacian_spectrum(graph), digits = 4)
            return([
                dash_datatable(id = "datatableGraph",
                columns = [Dict("name" => string(vtxdict[i]), "id" => string(i)) for i ∈ 1:length(aux)],
                data = [Dict(string(i) => aux[i] for i ∈ 1:length(aux))],
                page_size = 20,
                sort_action = "native",
                style_table = Dict("height" => "auto",
                "width" => "auto", "overflowY" => "auto",
                "overflowX" => "auto"),
                style_cell = Dict("minWidth" => "50px", "fontSize" => "16px"),
                style_header = Dict("fontWeight" => "bold",
                "backgroundColor" => "rgb(230, 230, 230)")),
                html_br(),
                html_br(),
                html_b("Laplacian Eigenvalues:",
                style = Dict("color" => "Tomato")),
                dcc_markdown("""
                The second smallest eigenvalue is the *Algebraic Connectivity* of the graph.
                """,
                style = Dict("textAlign" => "justify"))
            ])
        #! Graph Property: Fiedler's Vector
        elseif property == "selectFiedler"
            aux = eigen(Symmetric(Matrix(laplacian_matrix(graph))), 2:2)
            return([
                dash_datatable(id = "datatableGraph",
                columns = [Dict("name" => string(vtxdict[i]), "id" => string(i)) for i ∈ 1:length(aux.vectors)],
                data = [Dict(string(i) => round(aux.vectors[i], digits = 4) for i ∈ 1:length(aux.vectors))],
                page_size = 20,
                sort_action = "native",
                style_table = Dict("height" => "auto",
                "width" => "auto", "overflowY" => "auto",
                "overflowX" => "auto"),
                style_cell = Dict("minWidth" => "50px", "fontSize" => "16px"),
                style_header = Dict("fontWeight" => "bold",
                "backgroundColor" => "rgb(230, 230, 230)")),
                html_br(),
                html_br(),
                html_b("Fidler's Vector:",
                style = Dict("color" => "Tomato")),
                dcc_markdown("""
                The eigenvector associated with the second smallest eigenvalue is the *Fiedler's Vector* of the graph.
                """,
                style = Dict("textAlign" => "justify"))
            ])
        #! Graph Property: Reliability
        elseif property == "selectReliability"
            return([
                html_b("Reliability:",
                style = Dict("color" => "Tomato")),
                dcc_markdown("""
                The *Node Reliability* is given by the following function:.
                """,
                style = Dict("textAlign" => "justify")),
                html_img(src = "https://i.imgur.com/It7MGUC.png", width = "100%")
            ])
        else
            return("Select an option on the menu above.")
        end
    else
        return("Select an option on the menu above.")
    end
end

########################################################################################################################
###                                                     Server                                                       ###
########################################################################################################################
run_server(app, "0.0.0.0", 8080, debug = true)

max_stars_repo_name: fausto-mpj/EduGraph

max_issues_repo_name: fausto-mpj/EduGraph

max_forks_repo_name: fausto-mpj/EduGraph


text:
module RandomWalk

export randomwalk2d, dlaline!

"""
2D Random Walk

# Arguments

-`time::Integer`: the number of steps of the random walk
-`walks::Integer`: the number of walks, defaults to the value of `time`

# Returns

two tuples of two integers giving the x and y limits of the points and
an integer matrix counting the number of times the walker landed on each point.
"""
function randomwalk2d(time::Integer; walks::Integer=100*time)
    x, y = collect(Integer, -time:time), collect(Integer, -time:time)
    grid = zeros(Integer, 2 * time + 1, 2 * time + 1)
    steps = ((0, 1), (1, 0), (0, -1), (-1, 0))

    for _ in 1:walks
        xpos, ypos = time + 1, time + 1
        for _ in 1:time
            step = steps[rand(1:4)]
            xpos += step[1]
            ypos += step[2]
        end
        grid[xpos, ypos] += 1
    end

    return randomwalk2dcrop(x, y, grid)
end

"""
crop the output of `randomwalk2d`
"""
@inline function randomwalk2dcrop(x::Vector{<:Integer}, y::Vector{<:Integer},
        grid::Matrix{<:Integer})
    height, width = size(grid)
    bottom = left = 1
    top = height
    right = width

    while all(==(0), grid[top, :])
        top -= 1
    end
    while all(==(0), grid[bottom, :])
        bottom += 1
    end
    while all(==(0), grid[:, left])
        left += 1
    end
    while all(==(0), grid[:, right])
        right -= 1
    end

    return x[left:right], y[bottom:top], grid[bottom:top, left:right]
end

"""
Applies diffusion-limited aggregation on a 2D surface, assuming the "seed" is a line

the line is at the bottom and new particles appear from the top of the pile

# Arguments

- `surface::BitMatrix`: The simulation surface. Particles occupy points that are true.
- `particles::Integer`: the number of particles (walkers) simulated
- `escapethresh::Integer`: if the particle gets this much higher from the pile, it will
    "escape". It will be removed from the simulation and a new particle will be generated
    in its place. Defaults to 5.

# Returns
the `surface`
"""
function dlaline!(surface::BitMatrix, particles::Integer; escapethresh::Integer=5)
    height, width = size(surface)
    maxheight = height
    while true ∉ surface[maxheight - 1, :]
        maxheight -= 1
    end
    steps = ((0, 1), (1, 0), (0, -1), (-1, 0))

    while particles > 0
        ypos, xpos = maxheight, rand(1:width)

        while true
            xprev, yprev = xpos, ypos
            step = steps[rand(1:4)]
            ypos += step[1]
            xpos += step[2]
            # periodic boundary conditions (if xpos = width + 1, then it will be shifted
            # left by width, and if xpos = 0, it will be shifted right by width)
            xpos += (width + 1 - 2xpos) ÷ (width + 1) * width

            if ypos > min(maxheight + escapethresh, height)
                break
            elseif surface[ypos, xpos]
                surface[yprev, xprev] = true
                particles -= 1
                if yprev > maxheight
                    maxheight = yprev
                end
                break
            end
        end
    end
    return surface # for convenience
end

end

max_stars_repo_name: slhshamloo/comp-phys

max_issues_repo_name: slhshamloo/comp-phys

max_forks_repo_name: slhshamloo/comp-phys


text:
# -----------------------------------------------------
# -----------------------------------------------------

# see: http://www.research.att.com/export/sites/att_labs/groups/infovis/res/legacy_papers/DBLP-journals-camwa-Koren05.pdf
# also: http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.3.2055&rep=rep1&type=pdf

function spectral_graph(adjmat::AbstractMatrix; node_weights::AbstractVector = ones(size(adjmat,1)), kw...)
    positions = NetworkLayout.Spectral.layout(adjmat; node_weights=convert(Vector{Float64}, node_weights))

    ([p[1] for p in positions], [p[2] for p in positions], [p[3] for p in positions])
end

function spectral_graph(source::AbstractVector{Int}, destiny::AbstractVector{Int}, weights::AbstractVector; kw...)
    spectral_graph(get_adjacency_matrix(source, destiny, weights); kw...)
end

function spring_graph(adjmat::AbstractMatrix;
                      dim = 2,
                      x = rand(size(adjmat)[1]),
                      y = rand(size(adjmat)[1]),
                      z = rand(size(adjmat)[1]),
                      maxiter = 100, initialtemp = 2.0, C = 2.0, kw...
                      )
    @assert dim == 2 || dim == 3
    T = Float64
    adjmat = make_symmetric(adjmat)
    startpostions = if dim == 2
        [Point(T(x[i]), T(y[i])) for i in 1:length(x)]
    elseif dim == 3
        [Point(T(x[i]), T(y[i]), T(z[i])) for i in 1:length(x)]
    end

    positions = NetworkLayout.Spring.layout(adjmat, GeometryTypes.Point{dim, Float32};
                                            iterations = maxiter, initialtemp = initialtemp,
                                            C = C, startpositions = startpostions)
    if dim == 2
        ([p[1] for p in positions], [p[2] for p in positions], nothing)
    else
        ([p[1] for p in positions], [p[2] for p in positions], [p[3] for p in positions])
    end
end

function spring_graph(source::AbstractVector{Int}, destiny::AbstractVector{Int}, weights::AbstractVector; kw...)
    spring_graph(get_adjacency_matrix(source, destiny, weights); kw...)
end

function sfdp_graph(adjmat::AbstractMatrix;
                    dim = 2,
                    x = rand(size(adjmat)[1]),
                    y = rand(size(adjmat)[1]),
                    z = rand(size(adjmat)[1]),
                    maxiter = 100, tol = 1e-10, C = 1.0, K = 1.0, kw...
                    )
    @assert dim == 2 || dim == 3
    adjmat = make_symmetric(adjmat)
    T = Float64
    startpostions = if dim == 2
        [Point(T(x[i]), T(y[i])) for i in 1:length(x)]
    elseif dim == 3
        [Point(T(x[i]), T(y[i]), T(z[i])) for i in 1:length(x)]
    end

    positions = NetworkLayout.SFDP.layout(adjmat, dim,
                                          iterations = maxiter, tol = tol, C = C, K = K)
    if dim == 2
        ([p[1] for p in positions], [p[2] for p in positions], nothing)
    else
        ([p[1] for p in positions], [p[2] for p in positions], [p[3] for p in positions])
    end
end

function sfdp_graph(source::AbstractVector{Int}, destiny::AbstractVector{Int}, weights::AbstractVector; kw...)
    sfpd_graph(get_adjacency_matrix(source, destiny, weights); kw...)
end

function circular_graph(adjmat::AbstractMatrix;
                        dim = 2,
                        x = rand(size(adjmat)[1]),
                        y = rand(size(adjmat)[1]),
                        z = rand(size(adjmat)[1]),
                        kw...
                        )
    @assert dim == 2

    positions = NetworkLayout.Circular.layout(adjmat)

    ([p[1] for p in positions], [p[2] for p in positions], nothing)
end

function circular_graph(source::AbstractVector{Int}, destiny::AbstractVector{Int}, weights::AbstractVector; kw...)
    circular_graph(get_adjacency_matrix(source, destiny, weights); kw...)
end

function shell_graph(adjmat::AbstractMatrix;
                     dim = 2,
                     x = rand(size(adjmat)[1]),
                     y = rand(size(adjmat)[1]),
                     z = rand(size(adjmat)[1]),
                     nlist = nothing, kw...
                     )
    @assert dim == 2
    positions = NetworkLayout.Shell.layout(adjmat, nlist=nlist)

    ([p[1] for p in positions], [p[2] for p in positions], nothing)
end

function shell_graph(source::AbstractVector{Int}, destiny::AbstractVector{Int}, weights::AbstractVector; kw...)
    shell_graph(get_adjacency_matrix(source, destiny, weights); kw...)
end

# -----------------------------------------------------

# Axis-by-Axis Stress Minimization -- Yehuda Koren and David Harel
# See: http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.437.3177&rep=rep1&type=pdf

# # NOTES:
# #   - dᵢⱼ = the "graph-theoretical distance between nodes i and j"
# #         = Aᵢⱼ
# #   - kᵢⱼ = dᵢⱼ⁻²
# #   - b̃ᵢ = ∑ᵢ≠ⱼ ((x̃ⱼ ≤ x̃ᵢ ? 1 : -1) / dᵢⱼ)
# #   - need to solve for x each iteration: Lx = b̃

# # Solve for one axis at a time while holding the others constant.
# # dims is 2 (2D) or 3 (3D).  free_dims is a vector of the dimensions to update (for example if you fix y and solve for x)
# function by_axis_stress_graph(adjmat::AbstractMatrix, node_weights::AbstractVector = ones(size(adjmat,1));
#                               dims = 2, free_dims = 1:dims,
#                               x = rand(length(node_weights)),
#                               y = rand(length(node_weights)),
#                               z = rand(length(node_weights)))
#     adjmat = make_symmetric(adjmat)
#     L, D = compute_laplacian(adjmat, node_weights)

#     n = length(node_weights)
#     maxiter = 100 # TODO: something else

#     @assert dims == 2

#     @show adjmat L

#     for _ in 1:maxiter
#         x̃ = x
#         b̃ = Float64[sum(Float64[(i==j || adjmat[i,j] == 0) ? 0.0 : ((x̃[j] <= x̃[i] ? 1.0 : -1.0) / adjmat[i,j]) for j=1:n]) for i=1:n]
#         @show x̃ b̃
#         x = L \ b̃

#         xdiff = x - x̃
#         @show norm(xdiff)
#         if norm(xdiff) < 1e-4
#             info("converged. norm(xdiff) = $(norm(xdiff))")
#             break
#         end
#     end
#     @show x y
#     x, y, z
# end

norm_ij(X, i, j) = sqrt(sum(Float64[(v[i]-v[j])^2 for v in X]))
stress(X, dist, w, i, j) = w[i,j] * (norm_ij(X, i, j) - dist[i,j])^2
function stress(X, dist, w)
    tot = 0.0
    for i=1:size(X,1), j=1:i-1
        tot += stress(X, dist, w, i, j)
    end
    tot
end

# follows section 2.3 from http://link.springer.com/chapter/10.1007%2F978-3-540-31843-9_25#page-1
# Localized optimization, updates: x
function by_axis_local_stress_graph(adjmat::AbstractMatrix;
                              node_weights::AbstractVector = ones(size(adjmat,1)),
                              dim = 2, free_dims = 1:dim,
                              x = rand(length(node_weights)),
                              y = rand(length(node_weights)),
                              z = rand(length(node_weights)),
                              maxiter = 1000,
                              kw...)
    adjmat = make_symmetric(adjmat)
    n = length(node_weights)

    # graph-theoretical distance between node i and j (i.e. shortest path distance)
    # TODO: calculate a real distance
    dist = estimate_distance(adjmat)
    # @show dist

    # also known as kᵢⱼ in "axis-by-axis stress minimization".  the -2 could also be 0 or -1?
    w = dist .^ -2

    # in each iteration, we update one dimension/node at a time, reducing the total stress with each update
    X = dim == 2 ? (x, y) : (x, y, z)
    laststress = stress(X, dist, w)
    for k in 1:maxiter
        for p in free_dims
            for i=1:n
                numer, denom = 0.0, 0.0
                for j=1:n
                    i==j && continue
                    numer += w[i,j] * (X[p][j] + dist[i,j] * (X[p][i] - X[p][j]) / norm_ij(X, i, j))
                    denom += w[i,j]
                end
                if denom != 0
                    X[p][i] = numer / denom
                end
            end
        end

        # check for convergence of the total stress
        thisstress = stress(X, dist, w)
        if abs(thisstress - laststress) / abs(laststress) < 1e-6
            # info("converged. numiter=$k last=$laststress this=$thisstress")
            break
        end
        laststress = thisstress
    end

    dim == 2 ? (X..., nothing) : X
end

function by_axis_local_stress_graph(source::AbstractVector{Int}, destiny::AbstractVector{Int}, weights::AbstractVector; kw...)
    by_axis_local_stress_graph(get_adjacency_matrix(source, destiny, weights); kw...)
end

# -----------------------------------------------------

function buchheim_graph(adjlist::AbstractVector;
                    node_weights::AbstractVector = ones(length(adjlist)),
                    root::Symbol = :top,  # flow of tree: left, right, top, bottom
                    layers_scalar = 1.0,
                    layers = nothing,
                    dim = 2,
                    kw...)
    # @show adjlist typeof(adjlist)
    positions = NetworkLayout.Buchheim.layout(adjlist, nodesize = convert(Vector{Float64}, node_weights))
    Float64[p[1] for p in positions], Float64[p[2] for p in positions], nothing
end

# -----------------------------------------------------

function tree_graph(adjmat::AbstractMatrix; kw...)
    tree_graph(get_source_destiny_weight(adjmat)...; kw...)
end

function tree_graph(source::AbstractVector{Int}, destiny::AbstractVector{Int}, weights::AbstractVector;
                    node_weights::AbstractVector = ones(max(maximum(source), maximum(destiny))),
                    root::Symbol = :top,  # flow of tree: left, right, top, bottom
                    layers_scalar = 1.0,
                    layers = nothing,
                    positions = nothing,
                    dim = 2,
                    add_noise = true,
                    kw...)
    extrakw = Dict{Symbol,Any}(kw)
    # @show root layers positions dim add_noise extrakw
    n = length(node_weights)

    # @show root

    # TODO: compute layers, which get bigger as you go away from the root
    if layers == nothing
        # layers = rand(1:4, n)
        layers = compute_tree_layers2(source, destiny, n)
    end

    # reverse direction?
    if root in (:top, :right)
        layers = -layers
    end

    # add noise
    if add_noise
        layers = layers + 0.6rand(size(layers)...)
    end

    # TODO: normalize layers somehow so it's in line with distances
    layers .*= layers_scalar
    if dim == 2
        if root in (:top, :bottom)
            extrakw[:y] = layers
            extrakw[:free_dims] = if isnothing(positions)
                [1]
            else
                extrakw[:x] = positions
                Int[]
            end
        elseif root in (:left, :right)
            extrakw[:x] = layers
            # extrakw[:free_dims] = [2]
            extrakw[:free_dims] = if isnothing(positions)
                [2]
            else
                extrakw[:y] = positions
                Int[]
            end
        else
            error("unknown root: $root")
        end
    else
        error("3d not supported")
    end

    # now that we've fixed one dimension, let the stress algo solve for the other(s)
    by_axis_local_stress_graph(get_adjacency_matrix(source, destiny, weights);
                               node_weights=node_weights,
                               dim=dim,
                               extrakw...)
end

function adjlist_and_degrees(source, destiny, n)
    # build a list of children (adjacency list)
    alist = Vector{Int}[Int[] for i=1:n]
    indeg, outdeg = zeros(Int, n), zeros(Int, n)
    for (si,di) in zip(source, destiny)
        push!(alist[si], di)
        indeg[di] += 1
        outdeg[si] += 1
    end
    alist, indeg, outdeg
end

function compute_tree_layers(source, destiny, n)
    alist, indeg, outdeg = adjlist_and_degrees(source, destiny, n)

    # choose root to be the node with lots going out, but few coming in
    netdeg = outdeg - 50indeg
    idxs = sortperm(netdeg, rev=true)
    # rootidx = findmax(netdeg)
    # @show outdeg indeg netdeg idxs alist
    placed = Int[]

    layers = zeros(n)
    for i=1:n
        idx = shift!(idxs)

        # first, place this after its parents
        for j in placed
            if idx in alist[j]
                layers[idx] = max(layers[idx], layers[j] + 1)
            end
        end

        # next, shift its children lower
        for j in idxs
            if j in alist[idx]
                layers[j] = max(layers[j], layers[idx] + 1)
            end
        end

        push!(placed, idx)
    end
    layers
end

# an alternative algo to pick tree layers... generate a list of roots,
# and for each root, make a pass through the tree (without recurrency)
# and push the children below their parents
function compute_tree_layers2(source, destiny, n)
    alist, indeg, outdeg = adjlist_and_degrees(source, destiny, n)
    roots = filter(i->indeg[i]==0, 1:n)
    if isempty(roots)
        roots = [1]
    end

    layers = zeros(Int,n)
    for i in roots
        shift_children!(layers, alist, Int[], i)
    end

    # now that we've shifted children out, move parents closer to their closest children
    while true
        shifted = false
        for parent=1:n
            if !(isempty(alist[parent]))
                minidx = minimum(layers[child] for child in alist[parent])
                if layers[parent] < minidx - 1
                    shifted = true
                    layers[parent] = minidx - 1
                end
            end
        end
        shifted || break
    end

    layers
end

function shift_children!(layers, alist, placed, parent)
    for idx in alist[parent]
        if !(idx in placed) && layers[idx] <= layers[parent]
            layers[idx] = layers[parent] + 1
        end
    end
    for idx in alist[parent]
        if idx != parent && !(idx in placed)
            push!(placed, idx)
            shift_children!(layers, alist, placed, idx)
        end
    end
end

# -----------------------------------------------------

# TODO: maybe also implement Catmull-Rom Splines? http://www.mvps.org/directx/articles/catmull/

# -----------------------------------------------------

function arc_diagram(source::AbstractVector{Int},
                     destiny::AbstractVector{Int},
                     weights::AbstractVector;
                     kw...)
    X = unique(vcat(source, destiny))
    O = zeros(Int,length(X))
    X, O, O
end

# -----------------------------------------------------

function chord_diagram( source::AbstractVector{Int},
                        destiny::AbstractVector{Int},
                        weights::AbstractVector;
                        kw... )
    nodes = unique(vcat(source, destiny))
    N = length(nodes)
    δ = 2pi / N

    x = Array{Float64}(undef, N)
    y = Array{Float64}(undef, N)
    for i in 1:N
        x[i] = sin((i-1)*δ)
        y[i] = cos((i-1)*δ)
    end

    x, y, zeros(Int,length(x))
end


text:
"""
GenX: An Configurable Capacity Expansion Model
Copyright (C) 2021,  Massachusetts Institute of Technology
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
A complete copy of the GNU General Public License v2 (GPLv2) is available
in LICENSE.txt.  Users uncompressing this from an archive may not have
received this license file.  If not, see <http://www.gnu.org/licenses/>.
"""

@doc raw"""
	configure_cplex(solver_settings_path::String)

Reads user-specified solver settings from cplex\_settings.yml in the directory specified by the string solver\_settings\_path.

Returns a MathOptInterface OptimizerWithAttributes CPLEX optimizer instance to be used in the GenX.generate_model() method.

The CPLEX optimizer instance is configured with the following default parameters if a user-specified parameter for each respective field is not provided:

 - CPX\_PARAM\_EPRHS = 1e-6 (Constraint (primal) feasibility tolerances. See https://www.ibm.com/support/knowledgecenter/en/SSSA5P_12.5.1/ilog.odms.cplex.help/CPLEX/Parameters/topics/EpRHS.html)
 - CPX\_PARAM\_EPOPT = 1e-6 (Dual feasibility tolerances. See https://www.ibm.com/support/knowledgecenter/en/SSSA5P_12.5.1/ilog.odms.cplex.help/CPLEX/Parameters/topics/EpOpt.html)
 - CPX\_PARAM\_AGGFILL = 10 (Allowed fill during presolve aggregation. See https://www.ibm.com/support/knowledgecenter/en/SSSA5P_12.5.1/ilog.odms.cplex.help/CPLEX/Parameters/topics/AggFill.html)
 - CPX\_PARAM\_PREDUAL = 0 (Decides whether presolve should pass the primal or dual linear programming problem to the LP optimization algorithm. See https://www.ibm.com/support/knowledgecenter/en/SSSA5P_12.5.1/ilog.odms.cplex.help/CPLEX/Parameters/topics/PreDual.html)
 - CPX\_PARAM\_TILIM = 1e+75	(Limits total time solver. See https://www.ibm.com/support/knowledgecenter/en/SSSA5P_12.5.1/ilog.odms.cplex.help/CPLEX/Parameters/topics/TiLim.html)
 - CPX\_PARAM\_EPGAP = 1e-4	(Relative (p.u. of optimal) mixed integer optimality tolerance for MIP problems (ignored otherwise). See https://www.ibm.com/support/knowledgecenter/en/SSSA5P_12.5.1/ilog.odms.cplex.help/CPLEX/Parameters/topics/EpGap.html)
 - CPX\_PARAM\_LPMETHOD = 0 (Algorithm used to solve continuous models (including MIP root relaxation). See https://www.ibm.com/support/knowledgecenter/de/SSSA5P_12.7.0/ilog.odms.cplex.help/CPLEX/Parameters/topics/LPMETHOD.html)
 - CPX\_PARAM\_BAREPCOMP = 1e-8 (Barrier convergence tolerance (determines when barrier terminates). See https://www.ibm.com/support/knowledgecenter/en/SSSA5P_12.5.1/ilog.odms.cplex.help/CPLEX/Parameters/topics/BarEpComp.html)
 - CPX\_PARAM\_NUMERICALEMPHASIS = 0 (Numerical precision emphasis. See https://www.ibm.com/support/knowledgecenter/en/SSSA5P_12.5.1/ilog.odms.cplex.help/CPLEX/Parameters/topics/NumericalEmphasis.html)
 - CPX\_PARAM\_BAROBJRNG = 1e+75 (Sets the maximum absolute value of the objective function. See https://www.ibm.com/support/knowledgecenter/en/SSSA5P_12.5.1/ilog.odms.cplex.help/CPLEX/Parameters/topics/BarObjRng.html)
 - CPX\_PARAM\_SOLUTIONTYPE = 2 (Solution type for LP or QP. See https://www.ibm.com/support/knowledgecenter/hr/SSSA5P_12.8.0/ilog.odms.cplex.help/CPLEX/Parameters/topics/SolutionType.html)

"""
function configure_cplex(solver_settings_path::String)

    solver_settings = YAML.load(open(solver_settings_path))

    # Set solve to use CPLEX for MIP or LP problems
    # Optional setup parameters ############################################
    # Constraint (primal) feasibility tolerances. See https://www.ibm.com/support/knowledgecenter/en/SSSA5P_12.5.1/ilog.odms.cplex.help/CPLEX/Parameters/topics/EpRHS.html
    FeasibilityTol = 1e-6
    if (haskey(solver_settings, "Feasib_Tol"))
        FeasibilityTol = solver_settings["Feasib_Tol"]
    end

    # Dual feasibility tolerances. See https://www.ibm.com/support/knowledgecenter/en/SSSA5P_12.5.1/ilog.odms.cplex.help/CPLEX/Parameters/topics/EpOpt.html
    OptimalityTol = 1e-6
    if (haskey(solver_settings, "Optimal_Tol"))
        OptimalityTol = solver_settings["Optimal_Tol"]
    end

    # Allowed fill during presolve aggregation. See https://www.ibm.com/support/knowledgecenter/en/SSSA5P_12.5.1/ilog.odms.cplex.help/CPLEX/Parameters/topics/AggFill.html
    AggFill = 10
    if (haskey(solver_settings, "AggFill"))
        AggFill = solver_settings["AggFill"]
    end

    # Decides whether presolve should pass the primal or dual linear programming problem to the LP optimization algorithm. See https://www.ibm.com/support/knowledgecenter/en/SSSA5P_12.5.1/ilog.odms.cplex.help/CPLEX/Parameters/topics/PreDual.html
    PreDual = 0
    if (haskey(solver_settings, "PreDual"))
        PreDual = solver_settings["PreDual"]
    end

    # Limits total time solver. See https://www.ibm.com/support/knowledgecenter/en/SSSA5P_12.5.1/ilog.odms.cplex.help/CPLEX/Parameters/topics/TiLim.html
    TimeLimit = 1e+75
    if (haskey(solver_settings, "TimeLimit"))
        TimeLimit = solver_settings["TimeLimit"]
    end

    # Relative (p.u. of optimal) mixed integer optimality tolerance for MIP problems (ignored otherwise). See https://www.ibm.com/support/knowledgecenter/en/SSSA5P_12.5.1/ilog.odms.cplex.help/CPLEX/Parameters/topics/EpGap.html
    MIPGap = 1e-4
    if (haskey(solver_settings, "MIPGap"))
        MIPGap = solver_settings["MIPGap"]
    end

    # Algorithm used to solve continuous models (including MIP root relaxation). See https://www.ibm.com/support/knowledgecenter/de/SSSA5P_12.7.0/ilog.odms.cplex.help/CPLEX/Parameters/topics/LPMETHOD.html
    Method = 0
    if (haskey(solver_settings, "Method"))
        Method = solver_settings["Method"]
    end

    # Barrier convergence tolerance (determines when barrier terminates). See https://www.ibm.com/support/knowledgecenter/en/SSSA5P_12.5.1/ilog.odms.cplex.help/CPLEX/Parameters/topics/BarEpComp.html
    BarConvTol = 1e-8
    if (haskey(solver_settings, "BarConvTol"))
        BarConvTol = solver_settings["BarConvTol"]
    end

    # Numerical precision emphasis. See https://www.ibm.com/support/knowledgecenter/en/SSSA5P_12.5.1/ilog.odms.cplex.help/CPLEX/Parameters/topics/NumericalEmphasis.html
    NumericFocus = 0
    if (haskey(solver_settings, "NumericFocus"))
        NumericFocus = solver_settings["NumericFocus"]
    end

    # Sets the maximum absolute value of the objective function. See https://www.ibm.com/support/knowledgecenter/en/SSSA5P_12.5.1/ilog.odms.cplex.help/CPLEX/Parameters/topics/BarObjRng.html
    BarObjRng = 1e+75
    if (haskey(solver_settings, "BarObjRng"))
        BarObjRng = solver_settings["BarObjRng"]
    end

    # Solution type for LP or QP. See https://www.ibm.com/support/knowledgecenter/hr/SSSA5P_12.8.0/ilog.odms.cplex.help/CPLEX/Parameters/topics/SolutionType.html
    SolutionType = 2
    if (haskey(solver_settings, "SolutionType"))
        SolutionType = solver_settings["SolutionType"]
    end
    ########################################################################

    OPTIMIZER = optimizer_with_attributes(
        CPLEX.Optimizer,
        "CPX_PARAM_EPRHS" => FeasibilityTol,
        "CPX_PARAM_EPOPT" => OptimalityTol,
        "CPX_PARAM_AGGFILL" => AggFill,
        "CPX_PARAM_PREDUAL" => PreDual,
        "CPX_PARAM_TILIM" => TimeLimit,
        "CPX_PARAM_EPGAP" => MIPGap,
        "CPX_PARAM_LPMETHOD" => Method,
        "CPX_PARAM_BAREPCOMP" => BarConvTol,
        "CPX_PARAM_NUMERICALEMPHASIS" => NumericFocus,
        "CPX_PARAM_BAROBJRNG" => BarObjRng,
        "CPX_PARAM_SOLUTIONTYPE" => SolutionType,
    )
    return OPTIMIZER
end

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM


text:
# function Vnk_NGG_old(n, k, β, σ)
#     β_arb = RR(β)
#     σ_arb = RR(σ)
#     return exp(β_arb) * σ_arb^(k-1) // RR(gamma(n)) * sum([binom(n-1, i, RR) * (-1)^i * β_arb ^(i//σ_arb) * real(gamma(k-i//CC(σ_arb), CC(β_arb))) for i in 0:(n-1)])
# end

"""
    Vnk_NGG(n, k, β, σ)

Calculation of the Vnk with arbitrary precision. This function is less accurate for small values of k.

# Examples
```julia-repl
julia> GibbsTypePriors.Vnk_NGG(950,50, 0.5, 0.2)
[3.8352033538906132801013353906868979398003615464023850790896863425714903553677226925055088746374316383109460954810284197859760984821731337318749595774947223616262209073739394021758498274625163338449050351539476265420058089446771803706284551125841735045167180428032503440230421477960661771838407576902892779249509437591900661480877143065301185622812541265406061355208499807014555543564206395963786037047501112201458115571178196894484765277010637093746870221360851774075869753015052347546075947342163011788718201438701670521824067688930574380546890246338996488357288469652851508886233959379599363414674816239620509492570297419371380404876526047019213064592743056547790007323431419320067128517155931031446743616067264862092384235340693861382516829416661260795864740780157957535281924044152165190002256884051072596878645089791000193415582741991189356825795269080464284696216369289690437152742813909679764874302349498025317865584922699306108146943457154750027744867148116339573187411564950371482234809283153942745404958070388914430591794906329509545071131596480723325249062012958207542739071651236894164675019841047663212800147578131631060869316509586195876115342586677261181109184501899398551285655781119232061487943017880727759501463530477776e-2387 +/- 8.27e-3625]
```
"""
function Vnk_NGG(n, k, β, σ)
    β_arb = RR(β)
    σ_arb = RR(σ)
    n_m_1_Flint = FlintZZ(n-1)
    return exp(β_arb) * σ_arb^(k-1) // RR(gamma(n)) * sum(binomial(n_m_1_Flint, FlintZZ(i)) * (-1)^i * β_arb ^(i//σ_arb) * real(gamma(k-i//CC(σ_arb), CC(β_arb))) for i in 0:(n-1))
end
#
# function Vnk_NGG2(n, k, β, σ)
#     gamma_n::arb = gamma(n)
#     return exp(β) * σ^(k-1) // gamma_n * sum([binom(n-1, i, RR) * (-1)^i * β^(i//σ) * real(gamma(k-i//CC(σ), CC(β))) for i in 0:(n-1)])
# end
# function Vnk_NGG3(n, k, β, σ)
#     gamma_n::arb = gamma(n)
#     return exp(β) * σ^(k-1) // gamma_n * sum(binom(n-1, i, RR) * (-1)^i * β^(i//σ) * real(gamma(k-i//CC(σ), CC(β))) for i in 0:(n-1))
# end
# function Vnk_NGG4(n, k, β, σ)
#     gamma_n::arb = gamma(n)
#     n_m_1_Flint = FlintZZ(n-1)
#     return exp(β) * σ^(k-1) // gamma_n * sum(binomial(n_m_1_Flint, FlintZZ(i)) * (-1)^i * β^(i//σ) * real(gamma(k-i//CC(σ), CC(β))) for i in 0:(n-1))
# end

# @memoize function Vnk_NGG_rec(n, k, β, σ)
#     # println("n = $n, k = $k")
#     if n < 100 || k ≥ n/2
#         return Vnk_NGG(n, k, β, σ)
#     else
#         return (Vnk_NGG_rec(n-1, k, β, σ)-Vnk_NGG_rec(n, k+1, β, σ))/RR(n-1-σ*k)
#     end
# end

@memoize function Vnk_NGG_rec(n, k, β, σ)
    res = Vnk_NGG(n, k, β, σ)
    if has_reasonable_precision(res)
        return res
    else
        println("n = $n, k = $k, Accuracy = $(accuracy_bits(res))")
        return (Vnk_NGG_rec(n-1, k, β, σ)-Vnk_NGG_rec(n, k+1, β, σ))/RR(n-1-σ*k)
    end
end

"""
    Vnk_2PD(n, k, θ, σ)

Calculation of the Vnk for a 2-parameter Poisson-Dirichlet (also known as Pitman-Yor) process with arbitrary precision.

# Examples
```julia-repl
julia> GibbsTypePriors.Vnk_2PD(7, 6, 0.5, 0.01)
[1.9770514285714285829380091500954323693644674824124026984992171309957729383892396509600074280486078176779783843427919567657142584744842411272863181245428808578583616597575557259694224232057348431696299003094318048035523344218002322698832143887862395401395274358633297792959324407688665501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501165501e-5 +/- 2.28e-1510]
```
"""
function Vnk_2PD(n, k, θ, σ)
    θ_arb::arb = RR(θ)
    σ_arb::arb = RR(σ)
    if k==1
        num::arb = arb_1
    else
        num = prod([θ_arb + i * σ_arb for i in 1:(k-1)])
    end
    return num // risingfac(1+θ_arb,n-1)
end

max_stars_repo_name: konkam/GibbsTypePriors

max_issues_repo_name: konkam/GibbsTypePriors

max_forks_repo_name: konkam/GibbsTypePriors