text:
# File: nlopt.jl
# Author: Jannis Harder
# Description: Bindings to NLopt

let
    Option(T...) = Union((), T...)

    global NLoptConstraint
    type NLoptConstraint
        f::Function
        tolerance::Float64
        equality::Bool
    end

    global nlopt_equality
    function nlopt_equality(f, t)
        NLoptConstraint(f, t, true)
    end

    global nlopt_inequality
    function nlopt_inequality(f, t)
        NLoptConstraint(f, t, false)
    end
    nlopt_inequality(f) = nlopt_inequality(f, 0.0)

    global NLopt
    type NLopt
        algorithm::Option(Symbol)
        min_objective::Option(Function) # set only one of these two
        max_objective::Option(Function)
        lower_bounds::Option(Array{Float64,1}, Float64)
        upper_bounds::Option(Array{Float64,1}, Float64)
        constraints::Array{NLoptConstraint, 1}
        stopval::Option(Float64)
        ftol_rel::Option(Float64)
        ftol_abs::Option(Float64)
        xtol_rel::Option(Float64)
        xtol_abs::Option(Array{Float64,1}, Float64)
        maxeval::Option(Int)
        maxtime::Option(Float64)
        local_optimizer::Option(NLopt)
        initial_step::Option(Array{Float64,1}, Float64)
        population::Option(Int)
        srand::Option(Int)

        function NLopt(algorithm)
            new(algorithm, #algorithm
                (), (), (), (),
                Array(NLoptConstraint, 0), # constraint lists
                (), (), (), (), (), (), (), (), (), (), ())
        end
    end

    NLopt() = NLopt(())
    function NLopt(algorithm, args...)
        opt = NLopt(algorithm)
        if mod(length(args), 2) == 1
            if !isa(args[end], Symbol)
                error("nlopt: expected argument name")
            end
            error("nlopt: no value given for parameter ", args[end])
        end
        for i in 1:div(length(args), 2)
            sym = args[2*i-1]
            if !isa(sym, Symbol)
                error("nlopt: expected argument name")
            end
            expr = Expr(:., {:opt, Expr(:quote, {sym}, Any)}, Any)
            eval(:((opt, v) -> $expr = v))(opt, convert(fieldtype(opt, sym),args[2*i]))
        end
        opt
    end

    global nlopt_wrapper_version = (2,2,4,0)

    global nlopt_algorithms =
    [:gn_direct,
     :gn_direct_l,
     :gn_direct_l_rand,
     :gn_direct_noscal,
     :gn_direct_l_noscal,
     :gn_direct_l_rand_noscal,
     :gn_orig_direct,
     :gn_orig_direct_l,
     :gd_stogo,
     :gd_stogo_rand,
     :ld_lbfgs_nocedal,
     :ld_lbfgs,
     :ln_praxis,
     :ld_var1,
     :ld_var2,
     :ld_tnewton,
     :ld_tnewton_restart,
     :ld_tnewton_precond,
     :ld_tnewton_precond_restart,
     :gn_crs2_lm,
     :gn_mlsl,
     :gd_mlsl,
     :gn_mlsl_lds,
     :gd_mlsl_lds,
     :ld_mma,
     :ln_cobyla,
     :ln_newuoa,
     :ln_newuoa_bound,
     :ln_neldermead,
     :ln_sbplx,
     :ln_auglag,
     :ld_auglag,
     :ln_auglag_eq,
     :ld_auglag_eq,
     :ln_bobyqa,
     :gn_isres,
     :auglag,
     :auglag_eq,
     :g_mlsl,
     :g_mlsl_lds,
     :ld_slsqp]

    algorithm_map = HashTable()

    for i in 1:length(nlopt_algorithms)
        algorithm_map[nlopt_algorithms[i]] = i - 1
    end

    # makes wrapping of c code simpler
    function wrap_dlopen(name, prefix)
        dl = dlopen(name)
        function wrap(sym, ret, types...)
            symptr = dlsym(dl, symbol(strcat(prefix, sym)))
            nargs = length(types)
            nsyms = map(i -> symbol(strcat("a", i)), [1:length(types)])
            eval(:(($nsyms...,) -> ccall($symptr, $ret, ($types...,), $nsyms...)))
        end
        return wrap
    end

    isnull(x) = int(x) == 0
    isset(x) = () != x

    run(`make -C $(require_resource("")) -s`)
    wrapper_dl = wrap_dlopen(require_resource("nlopt_wrapper.so"), :nlopt_wrapper_)

    version = Array(Int32, 7)
    wrapper_dl(:version, Void, Ptr{Int32})(version)

    nlopt_wrapper_lib_version = tuple(int(version[4:])...)
    global nlopt_lib_version = tuple(int(version[1:3])...)

    if nlopt_wrapper_lib_version != nlopt_wrapper_version
        error("nlopt: mismatch between wrapper version ",
              nlopt_wrapper_version, " and wrapper lib version ",
              nlopt_wrapper_lib_version)
    end
    if nlopt_lib_version != nlopt_wrapper_version[1:3]
        println("Warning: nlopt: mismatch between wrapper version ",
                nlopt_wrapper_version, " and nlopt version ",
                nlopt_lib_version)
    end

    create = wrapper_dl(:create, Ptr{Void}, Int32, Int32)
    free = wrapper_dl(:free, Void, Ptr{Void})
    objective = wrapper_dl(:objective, Int32, Ptr{Void}, Int32)
    optimize_start = wrapper_dl(:optimize_start, Int32, Ptr{Void}, Ptr{Float64})
    optimize_callback = wrapper_dl(:optimize_callback, Int32, Ptr{Void}, Ptr{Float64}, Ptr{Float64}, Float64, Int32, Ptr{Int32})
    optimize_finalize = wrapper_dl(:optimize_finalize, Int32, Ptr{Void}, Ptr{Float64}, Ptr{Float64})
    dimopt = wrapper_dl(:dimopt, Void, Ptr{Void}, Ptr{Float64}, Int32)
    doubleopt = wrapper_dl(:doubleopt, Void, Ptr{Void}, Float64, Int32)
    intopt = wrapper_dl(:intopt, Void, Ptr{Void}, Int, Int32)
    add_constraint = wrapper_dl(:add_constraint, Void, Ptr{Void}, Int32, Float64, Int32)
    local_optimizer = wrapper_dl(:local_optimizer, Void, Ptr{Void}, Ptr{Void})

    const st_idle = 0
    const st_running = 1
    const st_value = 2
    const st_grad = 3
    const st_done = 4

    function ensure_width(n, x)
        if isa(x, Float64)
            return (0 * [1:n]) + x
        elseif length(x) != n
            error("nlopt: dimension mismatch (got ", length(x), ", expected ", n)
        else
            return x
        end
    end

    function set_options(w, opt, n)
        if isset(opt.lower_bounds)
            dimopt(w, ensure_width(n, opt.lower_bounds), 0)
        end
        if isset(opt.upper_bounds)
            dimopt(w, ensure_width(n, opt.upper_bounds), 1)
        end
        if isset(opt.xtol_abs)
            dimopt(w, ensure_width(n, opt.xtol_abs), 1)
        end
        if isset(opt.stopval)
            doubleopt(w, opt.stopval, 0)
        end
        if isset(opt.ftol_rel)
            doubleopt(w, opt.ftol_rel, 1)
        end
        if isset(opt.ftol_abs)
            doubleopt(w, opt.ftol_abs, 2)
        end
        if isset(opt.xtol_rel)
            doubleopt(w, opt.xtol_rel, 3)
        end
        if isset(opt.maxtime)
            doubleopt(w, opt.maxtime, 4)
        end
        if isset(opt.maxeval)
            intopt(w, opt.maxeval, 0)
        end
        if isset(opt.population)
            intopt(w, opt.population, 1)
        end
        if isset(opt.srand)
            intopt(w, opt.srand, 2)
        end
    end

    function create_handle(opt, n)
        if opt.algorithm == ()
            error("nlopt: no algorithm specified")
        elseif !has(algorithm_map, opt.algorithm)
            error("nlopt: unknown algorithm ", opt.algorithm)
        end
        algorithm_id = algorithm_map[opt.algorithm]
        w = create(algorithm_id, n)
        if isnull(w)
            error("nlopt: could not initialize")
        end
        return w
    end

    function ignorable_arg(f)
        wrapped = false
        function f2(x, d)
            return f(x)
        end
        function f1(x, d)
            local r
            try
                r = f(x, d)
                f1 = f
            catch err
                f1 = f2
                r = f1(x, d)
            end
            return r
        end
        return (x, d) -> f1(x, d)
    end

    global nlopt_optimize
    function nlopt_optimize(opt, x::Array{Float64, 1})
        n = length(x)
        w = create_handle(opt, n)
        ret = ()
        try
            if isset(opt.max_objective) && isset(opt.min_objective)
                error("nlopt: only one of max_objective and min_objective can be used")
            elseif isset(opt.max_objective)
                objective(w, 1)
                objective_f = opt.max_objective
            elseif isset(opt.min_objective)
                objective(w, 0)
                objective_f = opt.min_objective
            else
                error("nlopt: no objective function specified")
            end

            set_options(w, opt, n)

            if isset(opt.local_optimizer)
                local_w = create_handle(opt.local_optimizer, n)
                try
                    set_options(local_w, opt.local_optimizer, n)
                catch err
                    free(local_w)
                    throw(err)
                end
                local_optimizer(w, local_w)
                free(local_w)
            end

            functions = Array(Function, 1)
            functions[1] = ignorable_arg(objective_f)

            for c in opt.constraints
                push(functions, ignorable_arg(c.f))
                add_constraint(w, length(functions), c.tolerance, c.equality ? 1 : 0)
            end

            optimize_start(w, x)
            x = Array(Float64, n)
            grad = Array(Float64, n)
            f = 0.0
            function_id = Array(Int32, 1)
            function_id[1] = 1
            exception = ()
            wrapped = false
            force_stop = 0
            while true
                status = optimize_callback(w, x, grad, f, force_stop, function_id)
                if status == st_done
                    break
                elseif status == st_value || status == st_grad
                    is_grad = status == st_grad
                    r = 0.0
                    try
                        r = functions[function_id[1]](x, is_grad)
                    catch err
                        exception = err
                        force_stop = 1
                        continue
                    end
                    if is_grad
                        (f, grad) = r
                    else
                        if isa(r, Tuple)
                            f = r[1]
                        else
                            f = r
                        end
                    end
                else
                    assert(false, "invalid state")
                end
            end
            f_ptr = Array(Float64, 1)
            res_code = optimize_finalize(w, x, f_ptr)
            res = {-1 => :failure,
                   -2 => :invalid_args,
                   -3 => :out_of_memory,
                   -4 => :roundoff_limited,
                   -5 => :forced_stop,
                   1 => :success,
                   2 => :stopval_reached,
                   3 => :ftol_reached,
                   4 => :xtol_reached,
                   5 => :maxeval_reached,
                   6 => :maxtime_reached}[res_code]
            if res == :forced_stop
                if isa(exception, Exception)
                    throw(exception)
                end
                res = exception
            elseif res_code < 0
                error("nlopt: optimization failed, reason: ", res)
            end
            ret = (res, x, f_ptr[1])
        catch x
            free(w)
            throw(x)
        end
        free(w)
        return ret
    end
end

max_stars_repo_name: MetalNinjas/julia-nlopt

max_issues_repo_name: MetalNinjas/julia-nlopt

max_forks_repo_name: MetalNinjas/julia-nlopt


text:
function solve_ode(u0, tspan, p, time_interval=1)
        function SIR(du,u,p,t)
                du[1]=-p[1]*u[1]*u[2]/(u[1]+u[2]+u[3])
                du[2]=p[1]*u[1]*u[2]/(u[1]+u[2]+u[3]) - p[2]*u[2]
                du[3] = p[2]*u[2]
        end
        problem = ODEProblem(SIR, u0, tspan,p)
        sol=solve(problem, saveat=time_interval)
        return sol
end

function solve_ode_seasonal_forcing(u0, tspan, p,time_interval=1)
        function SIR_forcing(du,u,p,t)
                L, beta0, beta1, theta, gamma = p
                # 17 weeks (~ 4 months) seasonal?
                beta_t = beta0*(1 + beta1*cos( (2*pi*t - theta)/26 ) )
                # du[1]=(Ntot-u[1]-u[2])/L - beta_t*u[2]*u[1]/Ntot
                # du[2]=beta_t*u[1]*u[2]/Ntot - u[2]/D
                du[1] =  -beta_t*u[2]*u[1]/(u[1]+u[2]+u[3]) + u[3]/L
                du[2] = beta_t*u[1]*u[2]/(u[1]+u[2]+u[3]) - gamma*u[2]
                du[3] = gamma*u[2]

        end
        problem = ODEProblem(SIR_forcing, u0, tspan,p)
        sol=solve(problem, saveat=time_interval)
        return sol
end

function solve_seir(u0, tspan, p, time_interval=1)
        function SEIR(du,u,p,t)
                du[1]=p[1]*(u[1]+u[2]+u[3]) - p[1]*u[1] - p[2]*u[1]*u[3]/(u[1]+u[2]+u[3]+u[4])
                du[2]=p[2]*u[1]*u[3]/(u[1]+u[2]+u[3]+u[4]) - (p[1] + p[3])*u[2]
                du[3] = p[3]*u[2] - (p[4]+p[1])*u[3]
                du[4] = p[4]*u[3] - p[1]*u[4]
        end
        problem = ODEProblem(SEIR, u0, tspan,p)
        sol=solve(problem, saveat=time_interval)
        return sol
end

ext: jl


text:
function breadthfirstsearch!(g::AbstractGraph, start::Vertex)
    q = LinkedQueue{typeof(start)}()
    enqueue!(q, start)
    setmark!(g, start, 1)
    while length(q) != 0
        v = dequeue!(q)
        #previsit(g, v)
        let w = first(g, v)
            while w <= countvertices(g)
                if getmark(g, w) == 0
                    setmark!(g, w, 1)
                    enqueue!(q, w)
                    #process(g, v, w)
                    println(v, w)
                end
                w = next(g, v, w)
            end
        end
    end
end


text:
# ---
# title: 979. Distribute Coins in Binary Tree
# id: problem979
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Tree, Depth-first Search
# link: <https://leetcode.com/problems/distribute-coins-in-binary-tree/description/>
# hidden: true
# ---
#
# Given the `root` of a binary tree with `N` nodes, each `node` in the tree has
# `node.val` coins, and there are `N` coins total.
#
# In one move, we may choose two adjacent nodes and move one coin from one node
# to another.  (The move may be from parent to child, or from child to parent.)
#
# Return the number of moves required to make every node have exactly one coin.
#
#
#
# **Example 1:**
#
# **![](https://assets.leetcode.com/uploads/2019/01/18/tree1.png)**
#
#
#
#     Input: [3,0,0]
#     Output: 2
#     Explanation: From the root of the tree, we move one coin to its left child, and one coin to its right child.
#
#
# **Example 2:**
#
# **![](https://assets.leetcode.com/uploads/2019/01/18/tree2.png)**
#
#
#
#     Input: [0,3,0]
#     Output: 3
#     Explanation: From the left child of the root, we move two coins to the root [taking two moves].  Then, we move one coin from the root of the tree to the right child.
#
#
# **Example 3:**
#
# **![](https://assets.leetcode.com/uploads/2019/01/18/tree3.png)**
#
#
#
#     Input: [1,0,2]
#     Output: 2
#
#
# **Example 4:**
#
# **![](https://assets.leetcode.com/uploads/2019/01/18/tree4.png)**
#
#
#
#     Input: [1,0,0,null,3]
#     Output: 4
#
#
#
#
# **Note:**
#
#   1. `1<= N <= 100`
#   2. `0 <= node.val <= N`
#
#
## @lc code=start
using LeetCode

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


text:
using PSIS
using Test
using Random
using ReferenceTests
using Distributions: GeneralizedPareto, Normal, Cauchy, Exponential, TDist, logpdf
using LogExpFunctions: logsumexp, softmax
using Logging: SimpleLogger, with_logger
using AxisArrays: AxisArrays

@testset "PSISResult" begin
    @testset "vector log-weights" begin
        log_weights = randn(500)
        log_weights_norm = logsumexp(log_weights)
        tail_length = 100
        reff = 2.0
        tail_dist = PSIS.GeneralizedPareto(1.0, 1.0, 0.5)
        result = PSISResult(log_weights, log_weights_norm, reff, tail_length, tail_dist)
        @test result isa PSISResult{Float64}
        @test sort(propertynames(result)) == [
            :log_weights,
            :log_weights_norm,
            :nchains,
            :ndraws,
            :nparams,
            :pareto_shape,
            :reff,
            :tail_dist,
            :tail_length,
            :weights,
        ]
        @test result.log_weights == log_weights
        @test result.log_weights_norm == log_weights_norm
        @test result.weights ≈ softmax(log_weights)
        @test result.reff == reff
        @test result.nparams == 1
        @test result.ndraws == 500
        @test result.nchains == 1
        @test result.tail_length == tail_length
        @test result.tail_dist == tail_dist
        @test result.pareto_shape == 0.5
        @test result.ess ≈ ess_is(result)

        @testset "show" begin
            @test sprint(show, "text/plain", result) == """
                PSISResult with 1 parameters, 500 draws, and 1 chains
                Pareto shape (k) diagnostic values:
                                    Count       Min. ESS
                 (-Inf, 0.5]  good  1 (100.0%)  $(floor(Int, result.ess))"""
        end
    end

    @testset "array log-weights" begin
        log_weights = randn(3, 500, 4)
        log_weights_norm = dropdims(logsumexp(log_weights; dims=(2, 3)); dims=(2, 3))
        tail_length = [1600, 1601, 1602]
        reff = [0.8, 0.9, 1.1]
        tail_dist = [
            PSIS.GeneralizedPareto(1.0, 1.0, 0.5),
            PSIS.GeneralizedPareto(1.0, 1.0, 0.6),
            PSIS.GeneralizedPareto(1.0, 1.0, 0.7),
        ]
        result = PSISResult(log_weights, log_weights_norm, reff, tail_length, tail_dist)
        @test result isa PSISResult{Float64}
        @test result.log_weights == log_weights
        @test result.log_weights_norm == log_weights_norm
        @test result.weights ≈ softmax(log_weights; dims=(2, 3))
        @test result.reff == reff
        @test result.nparams == 3
        @test result.ndraws == 500
        @test result.nchains == 4
        @test result.tail_length == tail_length
        @test result.tail_dist == tail_dist
        @test result.pareto_shape == [0.5, 0.6, 0.7]

        @testset "show" begin
            proposal = Normal()
            target = TDist(7)
            rng = MersenneTwister(42)
            x = rand(rng, proposal, 30, 100)
            log_ratios = logpdf.(target, x) .- logpdf.(proposal, x)
            reff = [100; ones(29)]
            result = psis(log_ratios, reff)
            @test sprint(show, "text/plain", result) == """
                PSISResult with 30 parameters, 100 draws, and 1 chains
                Pareto shape (k) diagnostic values:
                                        Count       Min. ESS
                 (-Inf, 0.5]  good       1 (3.3%)   95
                  (0.5, 0.7]  okay       3 (10.0%)  97
                    (0.7, 1]  bad        5 (16.7%)  ——
                    (1, Inf)  very bad  20 (66.7%)  ——
                          ——  missing    1 (3.3%)   ——"""
        end
    end
end

@testset "psis/psis!" begin
    @testset "importance sampling tests" begin
        target = Exponential(1)
        x_target = 1  # 𝔼[x] with x ~ Exponential(1)
        x²_target = 2  # 𝔼[x²] with x ~ Exponential(1)
        # For θ < 1, the closed-form distribution of importance ratios with k = 1 - θ is
        # GeneralizedPareto(θ, θ * k, k), and the closed-form distribution of tail ratios is
        # GeneralizedPareto(5^k * θ, θ * k, k).
        # For θ < 0.5, the tail distribution has no variance, and estimates with importance
        # weights become unstable
        @testset "Exponential($θ) → Exponential(1)" for (θ, atol) in [
            (0.8, 0.05), (0.55, 0.2), (0.3, 0.7)
        ]
            proposal = Exponential(θ)
            k_exp = 1 - θ
            for sz in ((100_000,), (5, 100_000), (5, 100_000, 4))
                dims = length(sz) == 1 ? Colon() : 2:length(sz)
                rng = MersenneTwister(42)
                x = rand(rng, proposal, sz)
                logr = logpdf.(target, x) .- logpdf.(proposal, x)

                r = psis(logr)
                @test r isa PSISResult
                logw = r.log_weights
                @test logw isa typeof(logr)

                if length(sz) == 3
                    @test all(r.tail_length .== PSIS.tail_length(1, 400_000))
                else
                    @test all(r.tail_length .== PSIS.tail_length(1, 100_000))
                end

                k = r.pareto_shape
                @test k isa (length(sz) == 1 ? Number : AbstractVector)
                tail_dist = r.tail_dist
                if length(sz) == 1
                    @test tail_dist isa PSIS.GeneralizedPareto
                    @test tail_dist.k == k
                else
                    @test tail_dist isa Vector{<:PSIS.GeneralizedPareto}
                    @test map(d -> d.k, tail_dist) == k
                end

                w = r.weights
                @test all(x -> isapprox(x, k_exp; atol=0.15), k)
                @test all(x -> isapprox(x, x_target; atol=atol), sum(x .* w; dims=dims))
                @test all(
                    x -> isapprox(x, x²_target; atol=atol), sum(x .^ 2 .* w; dims=dims)
                )
            end
        end
    end

    @testset "keywords" begin
        @testset "sorted=true" begin
            x = randn(100)
            perm = sortperm(x)
            @test psis(x).log_weights ==
                invpermute!(psis(x[perm]; sorted=true).log_weights, perm)
            @test psis(x).pareto_shape == psis(x[perm]; sorted=true).pareto_shape
        end
    end

    @testset "warnings" begin
        io = IOBuffer()
        logr = randn(5)
        result = with_logger(SimpleLogger(io)) do
            psis(logr)
        end
        @test result.log_weights == logr
        @test ismissing(result.tail_dist)
        @test ismissing(result.pareto_shape)
        msg = String(take!(io))
        @test occursin(
            "Warning: 1 tail draws is insufficient to fit the generalized Pareto distribution.",
            msg,
        )

        skipnan(x) = filter(!isnan, x)
        io = IOBuffer()
        for logr in [
            [NaN; randn(100)],
            [Inf; randn(100)],
            fill(-Inf, 100),
            vcat(ones(50), fill(-Inf, 435)),
        ]
            result = with_logger(SimpleLogger(io)) do
                psis(logr)
            end
            @test skipnan(result.log_weights) == skipnan(logr)
            @test ismissing(result.tail_dist)
            @test ismissing(result.pareto_shape)
            msg = String(take!(io))
            @test occursin("Warning: Tail contains non-finite values.", msg)
        end

        io = IOBuffer()
        rng = MersenneTwister(42)
        x = rand(rng, Exponential(50), 1_000)
        logr = logpdf.(Exponential(1), x) .- logpdf.(Exponential(50), x)
        result = psis(logr)
        result = with_logger(SimpleLogger(io)) do
            psis(logr)
        end
        @test result.log_weights != logr
        @test result.pareto_shape > 0.7
        msg = String(take!(io))
        @test occursin(
            "Warning: Pareto shape k = 0.73 > 0.7. $(PSIS.BAD_SHAPE_SUMMARY)", msg
        )

        io = IOBuffer()
        with_logger(SimpleLogger(io)) do
            PSIS.check_pareto_shape(PSIS.GeneralizedPareto(0.0, 1.0, 1.1))
        end
        msg = String(take!(io))
        @test occursin(
            "Warning: Pareto shape k = 1.1 > 1. $(PSIS.VERY_BAD_SHAPE_SUMMARY)", msg
        )

        io = IOBuffer()
        with_logger(SimpleLogger(io)) do
            PSIS.check_pareto_shape(PSIS.GeneralizedPareto(0.0, 1.0, 0.8))
        end
        msg = String(take!(io))
        @test occursin(
            "Warning: Pareto shape k = 0.8 > 0.7. $(PSIS.BAD_SHAPE_SUMMARY)", msg
        )

        io = IOBuffer()
        with_logger(SimpleLogger(io)) do
            PSIS.check_pareto_shape(PSIS.GeneralizedPareto(0.0, 1.0, 0.69))
        end
        msg = String(take!(io))
        @test isempty(msg)

        tail_dist = [
            missing,
            PSIS.GeneralizedPareto(0, 1, 0.69),
            PSIS.GeneralizedPareto(0, 1, 0.71),
            PSIS.GeneralizedPareto(0, 1, 1.1),
        ]
        io = IOBuffer()
        with_logger(SimpleLogger(io)) do
            PSIS.check_pareto_shape(tail_dist)
        end
        msg = String(take!(io))
        @test occursin(
            "Warning: 1 parameters had Pareto shape values 0.7 < k ≤ 1. $(PSIS.BAD_SHAPE_SUMMARY)",
            msg,
        )
        @test occursin(
            "Warning: 1 parameters had Pareto shape values k > 1. $(PSIS.VERY_BAD_SHAPE_SUMMARY)",
            msg,
        )
        @test occursin(
            "Warning: For 1 parameters, the generalized Pareto distribution could not be fit to the tail draws.",
            msg,
        )
    end

    @testset "test against reference values" begin
        rng = MersenneTwister(42)
        proposal = Normal()
        target = Cauchy()
        sz = (5, 1_000, 4)
        x = rand(rng, proposal, sz)
        logr = logpdf.(target, x) .- logpdf.(proposal, x)
        @testset for r_eff in (0.7, 1.2), improved in (true, false)
            r_effs = fill(r_eff, sz[1])
            result = psis(logr, r_effs; improved=improved)
            logw = result.log_weights
            @test !isapprox(logw, logr)
            basename = "normal_to_cauchy_reff_$(r_eff)"
            if improved
                basename = basename * "_improved"
            end
            @test_reference(
                "references/$basename.jld2",
                Dict("log_weights" => logw, "pareto_shape" => result.pareto_shape),
                by =
                    (ref, x) ->
                        isapprox(ref["log_weights"], x["log_weights"]; rtol=1e-6) &&
                            isapprox(ref["pareto_shape"], x["pareto_shape"]; rtol=1e-6),
            )
        end
    end

    # https://github.com/arviz-devs/PSIS.jl/issues/27
    @testset "no failure for very low log-weights" begin
        psis(rand(1000) .- 1500)
    end

    @testset "compatibility with arrays with named axes/dims" begin
        param_names = [Symbol("x[$i]") for i in 1:10]
        iter_names = 101:200
        chain_names = 1:4
        x = randn(length(param_names), length(iter_names), length(chain_names))

        @testset "AxisArrays" begin
            logr = AxisArrays.AxisArray(
                x,
                AxisArrays.Axis{:param}(param_names),
                AxisArrays.Axis{:iter}(iter_names),
                AxisArrays.Axis{:chain}(chain_names),
            )
            result = psis(logr)
            @test result.log_weights isa AxisArrays.AxisArray
            @test AxisArrays.axes(result.log_weights) == AxisArrays.axes(logr)
            for k in (:pareto_shape, :tail_length, :tail_dist, :reff, :log_weights_norm)
                prop = getproperty(result, k)
                @test prop isa AxisArrays.AxisArray
                @test AxisArrays.axes(prop) == (AxisArrays.axes(logr, 1),)
            end
        end
    end
end


text:
module Jeeves
    using LinearAlgebra, Combinatorics
    using Optim, Distributions
    using DataFrames: DataFrame
    using TexTables, OrderedCollections
    using RDatasets

    export fit
    export fit!
    export coef
    export summary
    export tidy

    export vcovIID, vcovCluster, vcovHetero, vcovBoot

    ## Types ##
    abstract type Model end
    abstract type LinearModel <: Model  end
    abstract type GeneralisedLinearModel <: Model end

    abstract type Fit <: Model end
    abstract type LinearModelFit <: Fit end

    abstract type vcov end

    ## Modules ##
    include("model.jl")
    include("ols.jl")
    include("tsls.jl")
    include("probit.jl")
    include("vcov.jl")
    include("vcovCluster.jl")
    include("vcovHetero.jl")
    include("vcovBoot.jl")
    include("utils.jl")
end # module

max_stars_repo_name: EdJeeOnGitHub/jeeves

max_issues_repo_name: EdJeeOnGitHub/jeeves

max_forks_repo_name: EdJeeOnGitHub/jeeves


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

using DataFrames, RData, CSV, XLSX, LinearAlgebra, Statistics, Intervals

dir = "X:/VIVES/1-Personal/Florian/git/OECD_ICIO/src/"
dir_raw = "C:/Users/u0148308/data/raw/" # location of raw data

# other scripts
include(dir * "functions.jl") # Script with functions to import and transform raw data

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

# OECD rev. 2021
source = "OECD"
revision = "2021"
year = 2015 # specified year
N = 69 # number of countries (originally MEX/CHN twice so 71)
S = 45 # number of industries

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

# Notes:
#   - Z: NS×NS, intermediate demand matrix
#   - Y: NS×N, final demand matrix
#   - W: NS×1, value added vector (notice: in OECD its a row vector)
#   - X: NS×1, output vector (notice: in OECD its a row vector)

Z, Y, W, X = import_data(dir_raw, source, revision, year)

Z, Y, W, X = remove_reorder(Z, Y, W, X, N, S)

any(Z .< 0.0)
any(Y .< 0.0) # possible to have negative values due to inventory!
any(W .< 0.0)
any(X .< 0.0)

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

# other basic matrices

A = Z ./ repeat(X', N*S) # NS×NS, input coefficients
A .= ifelse.(isnan.(A), 0.0, A) # in case X = 0 (otherwise we cannot compute Leontief inverse since NaN in A)

B = inv(I - A) # NS×NS, Leontief inverse

count(broadcast(in, B*[sum(Y[i,:]) for i in 1:N*S] ./ X, 0.95..1.05) .== 0) # gives number of values lying outside of interval

# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Crosscheck computations with TiVA database: https://stats.oecd.org/Index.aspx?datasetcode=TIVA_2021_C1#
# reporter:
#   - country: AUT (n = 2)
#   - industry: D16 (wood and wood products) (s = 8)
#   - row number: (n-1)*S+s = 1*S+8 = 53
# partner:
#   - country: DEU (n = 13)
#   - industry: D41T43 (construction) (s = 25)
#   - column number: (n-1)*S+s = 12*S+25 = 565
#   - column number: (n-1)*S+s = 12*S+8 = 548 (for imports, i.e. D16 as origin industry)
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# 4.1 PROD: Production (gross output), USD million
X

# 4.2 VALU: Value added, USD million
W

# 4.3 PROD_VASH: Value added as a share of Gross Output, by industry, percentage
V = W ./ X # NS×1, (row vector in OECD)

# 4.4 EXGR | EXGR_INT | EXGR_FNL: Gross exports, by industry and by partner country, USD million (f.o.b.)
#   - to obtain total => [sum(EXGR[i, :]) for i in 1:N]
EXGR, EXGR_INT, EXGR_FNL = exports(Z, Y, N, S, "country") # N×N
EXGR, EXGR_INT, EXGR_FNL = exports(Z, Y, N, S, "industry") # NS×N

# 4.5. IMGR | IMGR_INT | IMGR_FNL: Gross imports, by industry and by partner country, USD million (f.o.b.)
#   - to obtain total => [sum(IMGR[:, i]) for i in 1:N]
IMGR, IMGR_INT, IMGR_FNL = imports(Z, Y, N, S, "country") # N×N
IMGR_INT, IMGR_FNL = imports(Z, Y, N, S, "industry") # N×NS, N×N

# 4.6. BALGR | BALGR_INT | BALGR_FNL: Gross trade balance, by partner country, USD million (f.o.b.)
EXGR, EXGR_INT, EXGR_FNL = exports(Z, Y, N, S, "country") # N×N
IMGR, IMGR_INT, IMGR_FNL = imports(Z, Y, N, S, "country") # N×N

#   - to obtain total => [sum(BALGR[i, :]) for i in 1:N], [sum(BALGR[:, i]) for i in 1:N]
#   - look at one row: negative/positive => import more than export (trade deficit) / export more than import (trade surplus)
#   - look at one column: positive/negative => import more than export (trade deficit) / export more than import (trade surplus)
BALGR = EXGR .- IMGR' # N×N
BALGR_INT = EXGR_INT .- IMGR_INT' # N×N
BALGR_FNL = EXGR_FNL .- IMGR_FNL' # N×N

# 4.7. EXGRpSH: Gross exports, partner shares, by industry, percentage
EXGR, EXGR_INT, EXGR_FNL = exports(Z, Y, N, S, "industry") # NS×N

#   - total exports by industry
EXGR_TOTAL = [sum(EXGR[i, :]) for i in 1:N*S]

EXGRpSH_INT = EXGR_INT ./ repeat(EXGR_TOTAL, 1, N) .* 100 # what to do with NaN?
EXGRpSH_FNL = EXGR_FNL ./ repeat(EXGR_TOTAL, 1, N) .* 100
EXGRpSH = EXGR ./ repeat(EXGR_TOTAL, 1, N) .* 100

all([sum(EXGRpSH[i, :]) for i in 1:N*S] .≈ 100)

# 4.8. IMGRpSH: Gross imports, partner shares %, by industry, percentage (DONT UNDERSTAND!)
IMGR, IMGR_INT, IMGR_FNL = imports(Z, Y, N, S, "country") # N×N

#   - total imports
IMGR_TOTAL = [sum(IMGR[:, i]) for i in 1:N]

IMGRpSH_INT = IMGR_INT ./ repeat(IMGR_TOTAL', N) .* 100 # what to do with NaN?
IMGRpSH_FNL = IMGR_FNL ./ repeat(IMGR_TOTAL', N) .* 100
IMGRpSH = IMGR ./ repeat(IMGR_TOTAL', N) .* 100

all([sum(IMGRpSH[:, i]) for i in 1:N] .≈ 100)


text:
# Single embedding training file
using Knet, ArgParse, JLD
include("model3.jl")
include("preprocess.jl")
#include("generator.jl")

function test(param, state, data; perp=false)
    totloss = 0.0
    numofbatch = 0
    for sequence in data
        val = loss(param, deepcopy(state), sequence) #val = loss(param, state, sequence)
        totloss += (perp? exp(val) : val)
        numofbatch += 1
    end
    return totloss / numofbatch
end

function update!(param, state, sequence; lr=1.0, gclip=0.0)
    gloss = lossgradient(param, deepcopy(state), sequence) #    gloss = lossgradient(param, state, sequence)
    gscale = lr
    if gclip > 0
        gnorm = sqrt(mapreduce(sumabs2, +, 0, gloss))
        if gnorm > gclip
            gscale *= gclip / gnorm
        end
    end

    for k=1:length(param)
        axpy!(-gscale, gloss[k], param[k])
    end

    isa(state, Vector{Any}) || error("State should not be Boxed.")
    # TODO: is that really needed?
    for i=1:length(state)
        state[i] = AutoGrad.getval(state[i])
    end
end

function train!(param, state, data, o)
    for sequence in data
        # Only open for gradient check
        if o[:gcheck] > 0
            gradcheck(loss, param, copy(state), sequence; gcheck=o[:gcheck])
        end
        update!(param, state, copy(sequence); lr=o[:lr], gclip=o[:gclip]) #update!(param, state, sequence; lr=o[:lr], gclip=o[:gclip])
    end
end

function main(args=ARGS)
    s = ArgParseSettings()
    s.exc_handler = ArgParse.debug_handler
    @add_arg_table s begin
        ("--trainfile"; default="readme_data.txt" ;help="Training file")
        ("--devfile"; help="Dev file")
        ("--testfile"; help="Test file")
        ("--loadfile"; help="Initialize model from file")
        ("--savefile"; help="Save final model to file")
        ("--atype"; default=(gpu()>=0 ? "KnetArray{Float32}" : "Array{Float32}"); help="array type: Array for cpu, KnetArray for gpu")
        ("--layerconfig"; arg_type=Int; nargs='+'; default=[16]; help="Sizes of the one or more LSTM layers")
        ("--embedding"; arg_type=Int; default=200; help="embedding vector")
        ("--batchsize"; arg_type=Int; default=1; help="Minibatchsize")
        ("--epochs"; arg_type=Int; default=3; help="Number of epochs for training.")
        ("--winit"; arg_type=Float64; default=0.3; help="Initial weights set to winit*randn().")
        ("--decay"; arg_type=Float64; default=0.9; help="Learning rate decay.")
        ("--lr"; arg_type=Float64; default=4.0; help="Initial learning rate.")
        ("--single_embedding"; default=true)
        ("--gcheck"; arg_type=Int; default=0; help="Check N random gradients.")
        ("--gclip"; arg_type=Float64; default=0.0; help="Value to clip the gradient norm at.")
    end
    isa(args, AbstractString) && (args=split(args))
    o = parse_args(args, s; as_symbols=true)
    o[:atype] = eval(parse(o[:atype]))

    for (k,v) in o
        println("$k => $v")
    end
    println("Embedding config is single embedding\n")

    # Data preperation
    tdata = Data(o[:trainfile]; batchsize=o[:batchsize])
    vocabsize = length(tdata.word_to_index)

    # Devdata preperation
    dfile = (o[:devfile] != nothing ? o[:devfile] : o[:trainfile])
    ddata = Data(dfile; batchsize=o[:batchsize], word_to_index=tdata.word_to_index)

    # Model initialization
    param = initparams(o[:atype], o[:layerconfig], o[:embedding], vocabsize, o[:winit])
    state = initstates(o[:atype], o[:layerconfig], o[:batchsize])

    # inital loss value
    initial_loss = test(param, state, tdata; perp=true) #initial_loss = test(param, deepcopy(state), tdata; perp=true)
    println("Initial loss is $initial_loss")
    devloss = test(param, state, ddata; perp=true) #devloss = test(param, deepcopy(state), ddata; perp=true)
    println("Initial dev loss is $devloss")
    devlast = devbest = devloss

    # training started
    for epoch=1:o[:epochs]
        @time train!(param, state, tdata, o) # train!(param, deepcopy(state), tdata, o)
        devloss = test(param, state, ddata; perp=true)
        println("Dev loss for epoch $epoch : $devloss")

        if (epoch % 5) == 0
            trainloss = test(param, state, tdata; perp=true)
            println("Train loss for epoch $epoch : $trainloss")
        end
        # check whether model becomes better
        if devloss < devbest
            devbest = devloss
            if o[:savefile] != nothing
                saveparam = map(p->Array{Float32}, param)
                save(o[:savefile], "model", saveparam, "vocab", tdata.word_to_index, "config", o)
            end
        end

        if devloss > devlast
            o[:lr] *= o[:decay]
            info("New learning rate: $(o[:lr])")
        end
        devlast = devloss
        flush(STDOUT)
    end
end
main(ARGS)


text:
# ============================================================
# Kinetic Theory Methods
# Multi-dispatch for 1) Galerkin and 2) collocation methods
# with different particle velocity settings
# ============================================================

"""
Calculate conservative moments from distribution function
* single component: 1D1F1V, 1D2F1V, 2D1F2V, 2D2F2V
* mixture: 1D4F1V

"""
function uq_moments_conserve(
    f::AbstractArray{<:AbstractFloat,2},
    u::AbstractArray{<:AbstractFloat,1},
    ω::AbstractArray{<:AbstractFloat,1},
) # 1D1F1V

    w = zeros(eltype(f), 3, axes(f, 2))
    for j in axes(w, 2)
        w[:, j] .= KitBase.moments_conserve(f[:, j], u, ω)
    end

    return w

end

#--- 1D2F1V ---#
function uq_moments_conserve(
    h::AbstractArray{<:AbstractFloat,2},
    b::AbstractArray{<:AbstractFloat,2},
    u::AbstractArray{<:AbstractFloat,1},
    ω::AbstractArray{<:AbstractFloat,1},
)

    w = zeros(eltype(f), 3, axes(h, 2))
    for j in axes(w, 2)
        w[:, j] .= KitBase.moments_conserve(h[:, j], b[:, j], u, ω)
    end

    return w

end

#--- 2D1F2V ---#
function uq_moments_conserve(
    f::AbstractArray{<:AbstractFloat,3},
    u::AbstractArray{<:AbstractFloat,2},
    v::AbstractArray{<:AbstractFloat,2},
    ω::AbstractArray{<:AbstractFloat,2},
)

    w = zeros(eltype(f), 4, axes(h, 3))
    for j in axes(w, 2)
        w[:, j] .= KitBase.moments_conserve(f[:, :, j], u, v, ω)
    end

    return w

end

#--- 2D2F2V ---#
function uq_moments_conserve(
    h::AbstractArray{<:AbstractFloat,3},
    b::AbstractArray{<:AbstractFloat,3},
    u::AbstractArray{<:AbstractFloat,2},
    v::AbstractArray{<:AbstractFloat,2},
    ω::AbstractArray{<:AbstractFloat,2},
)

    w = zeros(eltype(f), 4, axes(h, 3))
    for j in axes(w, 2)
        w[:, j] .= KitBase.moments_conserve(h[:, :, j], b[:, :, j], u, v, ω)
    end

    return w

end

"""Multi-component gas"""

#--- 1D4F1V ---#
function uq_moments_conserve(
    h0::AbstractArray{Float64,3},
    h1::AbstractArray{Float64,3},
    h2::AbstractArray{Float64,3},
    h3::AbstractArray{Float64,3},
    u::AbstractArray{Float64,2},
    ω::AbstractArray{Float64,2},
)

    w = zeros(eltype(h0), 5, size(h0, 2), size(h0, 3))
    for j in axes(w, 2)
        w[:, j, :] .= KitBase.mixture_moments_conserve(
            h0[:, j, :],
            h1[:, j, :],
            h2[:, j, :],
            h3[:, j, :],
            u,
            ω,
        )
    end

    return w

end

"""
Calculate equilibrium distribution

"""
function uq_maxwellian(
    uspace::T1,
    prim::T2,
    uq::T3,
) where {T1<:AbstractArray{<:AbstractFloat,1},T2<:AbstractArray{<:Real,2},T3<:AbstractUQ} # 1D1F1V

    if size(prim, 2) == uq.nm + 1 && uq.nm + 1 != uq.nq

        primRan = chaos_ran(prim, 2, uq)

        MRan = zeros(axes(uspace, 1), axes(primRan, 2))
        for j in axes(MRan, 2)
            MRan[:, j] .= KitBase.maxwellian(uspace, primRan[:, j])
        end

        M = ran_chaos(MRan, 2, uq)

        return M

    elseif size(prim, 2) == uq.nq

        M = zeros(axes(uspace, 1), axes(prim, 2))
        for j in axes(M, 2)
            M[:, j] .= KitBase.maxwellian(uspace, prim[:, j])
        end

        return M

    else

        throw("inconsistent random domain size")

    end

end

#--- 2D1F2V ---#
function uq_maxwellian(
    u::T1,
    v::T1,
    prim::T2,
    uq::T3,
) where {T1<:AbstractArray{<:AbstractFloat,2},T2<:AbstractArray{<:Real,2},T3<:AbstractUQ}

    if size(prim, 2) == uq.nm + 1 && uq.nm + 1 != uq.nq

        primRan = chaos_ran(prim, 2, uq)

        MRan = zeros((axes(u)..., axes(primRan, 2)))
        for k in axes(MRan, 3)
            MRan[:, :, k] .= KitBase.maxwellian(u, v, primRan[:, k])
        end

        M = ran_chaos(MRan, 3, uq)

        return M

    elseif size(prim, 2) == uq.nq

        M = zeros((axes(u)..., axes(prim, 2)))
        for k in axes(M, 3)
            M[:, :, k] .= KitBase.maxwellian(u, v, prim[:, k])
        end

        return M

    else

        throw("inconsistent random domain size")

    end

end

#--- 2D2F2V ---#
function uq_maxwellian(
    u::T1,
    v::T1,
    prim::T2,
    uq::T3,
    inK,
) where {T1<:AbstractArray{<:AbstractFloat,2},T2<:AbstractArray{<:Real,2},T3<:AbstractUQ}

    if size(prim, 2) == uq.nm + 1 && uq.nm + 1 != uq.nq

        primRan = chaos_ran(prim, 2, uq)

        HRan = zeros((axes(u)..., axes(primRan, 2)))
        BRan = similar(HRan)
        for k in axes(HRan, 3)
            HRan[:, :, k] .= KitBase.maxwellian(u, v, primRan[:, k])
            BRan[:, :, k] .= HRan[:, :, k] .* inK ./ (2.0 * primRan[end, k])
        end

        H = ran_chaos(HRan, 3, uq)
        B = ran_chaos(BRan, 3, uq)

        return H, B

    elseif size(prim, 2) == uq.nq

        H = zeros((axes(u)..., axes(prim, 2)))
        B = similar(H)
        for k in axes(H, 3)
            H[:, :, k] .= KitBase.maxwellian(u, v, prim[:, k])
            B[:, :, k] .= H[:, :, k] .* inK ./ (2.0 * prim[end, k])
        end

        return H, B

    else

        throw("inconsistent random domain size in settings and solutions")

    end

end

#--- 3D1F3V ---#
function uq_maxwellian(
    u::T1,
    v::T1,
    w::T1,
    prim::T2,
    uq::T3,
) where {T1<:AbstractArray{<:AbstractFloat,3},T2<:AbstractArray{<:Real,2},T3<:AbstractUQ}

    if size(prim, 2) == uq.nm + 1 && uq.nm + 1 != uq.nq

        primRan = chaos_ran(prim, 2, uq)

        MRan = zeros((axes(u)..., axes(primRan, 2)))
        for k in axes(MRan, 4)
            MRan[:, :, :, k] .= KitBase.maxwellian(u, v, w, primRan[:, k])
        end

        M = ran_chaos(MRan, 4, uq)

        return M

    elseif size(prim, 2) == uq.nq

        M = zeros((axes(u)..., axes(prim, 2)))
        for k in axes(M, 4)
            M[:, :, :, k] .= KitBase.maxwellian(u, v, w, prim[:, k])
        end

        return M

    else

        throw("inconsistent random domain size in settings and solutions")

    end

end

"""Multi-component substances"""

#--- 1D4F1V ---#
function uq_maxwellian(
    uspace::T1,
    prim::T2,
    uq::T3,
) where {T1<:AbstractArray{<:AbstractFloat,2},T2<:AbstractArray{<:Real,3},T3<:AbstractUQ}

    if size(prim, 2) == uq.nm + 1 && uq.nm + 1 != uq.nq
        primRan = chaos_ran(prim, 2, uq)

        Mv = zeros(uq.nq, 2)
        Mw = similar(Mv)
        for k in axes(Mv, 2)
            for j in axes(Mv, 1)
                Mv[j, k] = primRan[3, j, k]^2 + 0.5 / primRan[end, j, k]
                Mw[j, k] = primRan[4, j, k]^2 + 0.5 / primRan[end, j, k]
            end
        end

        H0Ran = zeros(axes(uspace, 1), 1:uq.nq, axes(prim, 3))
        H1Ran = similar(H0Ran)
        H2Ran = similar(H0Ran)
        H3Ran = similar(H0Ran)
        for k in axes(H0Ran, 3)
            for j in axes(H0Ran, 2)
                H0Ran[:, j, k] .= KitBase.maxwellian(
                    uspace[:, k],
                    primRan[1, j, k],
                    primRan[2, j, k],
                    primRan[5, j, k],
                )
                H1Ran[:, j, k] .= primRan[3, j, k] .* H0Ran[:, j, k]
                H2Ran[:, j, k] .= primRan[4, j, k] .* H0Ran[:, j, k]
                H3Ran[:, j, k] .= (Mv[j, k] + Mw[j, k]) .* H0Ran[:, j, k]
            end
        end

        H0 = ran_chaos(H0Ran, 2, uq)
        H1 = ran_chaos(H1Ran, 2, uq)
        H2 = ran_chaos(H2Ran, 2, uq)
        H3 = ran_chaos(H3Ran, 2, uq)

        return H0, H1, H2, H3

    elseif size(prim, 2) == uq.nq

        Mv = zeros(axes(prim, 2), 2)
        Mw = similar(Mv)
        for k in axes(Mv, 2)
            for j in axes(Mv, 1)
                Mv[j, k] = prim[3, j, k]^2 + 0.5 / prim[end, j, k]
                Mw[j, k] = prim[4, j, k]^2 + 0.5 / prim[end, j, k]
            end
        end

        H0 = zeros(axes(uspace, 1), axes(prim, 2), axes(prim, 3))
        H1 = similar(H0)
        H2 = similar(H0)
        H3 = similar(H0)
        for k in axes(H0, 3)
            for j in axes(H0, 2)
                H0[:, j, k] .= KitBase.maxwellian(
                    uspace[:, k],
                    prim[1, j, k],
                    prim[2, j, k],
                    prim[5, j, k],
                )
                H1[:, j, k] .= prim[3, j, k] .* H0[:, j, k]
                H2[:, j, k] .= prim[4, j, k] .* H0[:, j, k]
                H3[:, j, k] .= (Mv[j, k] + Mw[j, k]) .* H0[:, j, k]
            end
        end

        return H0, H1, H2, H3

    else

        throw("inconsistent random domain size in settings and solutions")

    end

end

#--- 2D3F2V ---#
function uq_maxwellian(
    u::T1,
    v::T1,
    prim::T2,
    uq::T3,
) where {T1<:AbstractArray{<:AbstractFloat,3},T2<:AbstractArray{<:Real,3},T3<:AbstractUQ}

    if size(prim, 2) == uq.nm + 1 && uq.nm + 1 != uq.nq # galerkin

        primRan = chaos_ran(prim, 2, uq)

        H0Ran = zeros(axes(u, 1), axes(v, 2), 1:uq.nq, axes(prim, 3))
        H1Ran = similar(H0Ran)
        H2Ran = similar(H0Ran)
        for l in axes(H0Ran, 4)
            for k in axes(H0Ran, 3)
                H0Ran[:, :, k, l] .=
                    KitBase.maxwellian(u[:, :, l], v[:, :, l], primRan[:, k, l])
                H1Ran[:, :, k, l] .= H0Ran[:, :, k, l] .* primRan[4, k, l]
                H2Ran[:, :, k, l] .=
                    H0Ran[:, :, k, l] .* (primRan[4, k, l]^2 + 1.0 / (2.0 * prim[5, k, l]))
            end
        end

        H0 = ran_chaos(H0Ran, 3, uq)
        H1 = ran_chaos(H1Ran, 3, uq)
        H2 = ran_chaos(H2Ran, 3, uq)

        return H0, H1, H2

    elseif size(prim, 2) == uq.nq # collocation

        H0 = zeros(axes(u, 1), axes(u, 2), axes(prim, 2), axes(prim, 3))
        H1 = similar(H0)
        H2 = similar(H0)
        for l in axes(H0, 4)
            for k in axes(H0, 3)
                H0[:, :, k, l] .= KitBase.maxwellian(u[:, :, l], v[:, :, l], prim[:, k, l])
                H1[:, :, k, l] .= H0[:, :, k, l] .* prim[4, k, l]
                H2[:, :, k, l] .=
                    H0[:, :, k, l] .* (prim[4, k, l]^2 + 1.0 / (2.0 * prim[5, k, l]))
            end
        end

        return H0, H1, H2

    else

        throw("inconsistent random domain size")

    end

end

"""
Transform reduced distribution functions

"""
function uq_energy_distribution(h::AbstractMatrix{T}, prim, K, uq) where {T}
    if size(h, 2) == uq.nm + 1 && uq.nm + 1 != uq.nq
        hRan = chaos_ran(h, 2, uq)
        bRan = zero(hRan)
        for j in axes(hRan, 2)
            @. bRan[:, j] = hRan[:, j] * K / (2.0 * prim[end, j])
        end
        b = ran_chaos(bRan, 2, uq)
    else
        b = zero(h)
        for j in axes(h, 2)
            @. b[:, j] = h[:, j] * K / (2.0 * prim[end, j])
        end
    end

    return b
end

function uq_energy_distribution(h::AbstractArray{T,3}, prim, K, uq) where {T}
    if size(h, 2) == uq.nm + 1 && uq.nm + 1 != uq.nq
        hRan = chaos_ran(h, 2, uq)
        bRan = zero(hRan)
        for j in axes(hRan, 3)
            @. bRan[:, :, j] = hRan[:, :, j] * K / (2.0 * prim[end, j])
        end
        b = ran_chaos(bRan, 3, uq)
    else
        b = zero(h)
        for j in axes(h, 3)
            @. b[:, :, j] = h[:, :, j] * K / (2.0 * prim[end, j])
        end
    end

    return b
end

"""
Calculate primitive -> conservative variables

"""
function uq_prim_conserve(
    prim::AbstractArray{<:AbstractFloat,2},
    gamma::Real,
    uq::AbstractUQ,
) # single component

    if size(prim, 2) == uq.nm + 1 && uq.nm + 1 != uq.nq

        primRan = chaos_ran(prim, 2, uq)

        wRan = similar(primRan)
        for j in axes(wRan, 2)
            wRan[:, j] .= KitBase.prim_conserve(primRan[:, j], gamma)
        end

        wChaos = zeros(axes(prim))
        for i in axes(wChaos, 1)
            wChaos[i, :] .= ran_chaos(wRan[i, :], uq)
        end

        return wChaos

    elseif size(prim, 2) == uq.nq

        wRan = similar(prim)
        for j in axes(wRan, 2)
            wRan[:, j] .= KitBase.prim_conserve(prim[:, j], gamma)
        end

        return wRan

    else

        throw("inconsistent random domain size in settings and solutions")

    end

end

#--- multiple component ---#
function uq_prim_conserve(
    prim::AbstractArray{<:AbstractFloat,3},
    gamma::Real,
    uq::AbstractUQ,
)

    if size(prim, 2) == uq.nm + 1 && uq.nm + 1 != uq.nq

        primRan = chaos_ran(prim, 2, uq)

        wRan = similar(primRan)
        for k in axes(wRan, 3)
            for j in axes(wRan, 2)
                wRan[:, j, k] .= KitBase.prim_conserve(primRan[:, j, k], gamma)
            end
        end

        wChaos = zeros(axes(prim))
        for k in axes(wChaos, 3)
            for i in axes(wChaos, 1)
                wChaos[i, :, k] .= ran_chaos(wRan[i, :, k], uq)
            end
        end

        return wChaos

    elseif size(prim, 2) == uq.nq

        wRan = similar(prim)
        for k in axes(wRan, 3)
            for j in axes(wRan, 2)
                wRan[:, j, k] .= KitBase.prim_conserve(prim[:, j, k], gamma)
            end
        end

        return wRan

    else

        throw("inconsistent random domain size in settings and solutions")

    end

end

"""
Calculate conservative -> primitive variables

"""
function uq_conserve_prim(w::AbstractArray{<:AbstractFloat,2}, gamma::Real, uq::AbstractUQ) # single component

    if size(w, 2) == uq.nm + 1 && uq.nm + 1 != uq.nq

        wRan = chaos_ran(w, 2, uq)

        primRan = similar(wRan)
        for j in axes(primRan, 2)
            primRan[:, j] .= KitBase.conserve_prim(wRan[:, j], gamma)
        end

        primChaos = zeros(axes(w))
        for i in axes(primChaos, 1)
            primChaos[i, :] .= ran_chaos(primRan[i, :], uq)
        end

        return primChaos

    elseif size(w, 2) == uq.nq

        primRan = similar(w)
        for j in axes(primRan, 2)
            primRan[:, j] .= KitBase.conserve_prim(w[:, j], gamma)
        end

        return primRan

    else

        throw("inconsistent random domain size in settings and solutions")

    end

end

#--- multiple component ---#
function uq_conserve_prim(w::AbstractArray{<:AbstractFloat,3}, gamma::Real, uq::AbstractUQ)

    if size(w, 2) == uq.nm + 1 && uq.nm + 1 != uq.nq

        wRan = chaos_ran(w, 2, uq)

        primRan = similar(wRan)
        for k in axes(primRan, 3)
            for j in axes(primRan, 2)
                primRan[:, j, k] .= KitBase.conserve_prim(wRan[:, j, k], gamma)
            end
        end

        primChaos = zeros(axes(w))
        for k in axes(primChaos, 3)
            for i in axes(primChaos, 1)
                primChaos[i, :, k] .= ran_chaos(primRan[i, :, k], uq)
            end
        end

        return primChaos

    elseif size(w, 2) == uq.nq

        primRan = similar(w)
        for k in axes(primRan, 3)
            for j in axes(primRan, 2)
                primRan[:, j, k] .= KitBase.conserve_prim(w[:, j, k], gamma)
            end
        end

        return primRan

    else

        throw("inconsistent random domain size in settings and solutions")

    end

end

function uq_conserve_prim!(sol::AbstractSolution, γ::Real, uq::AbstractUQ)

    for i in eachindex(sol.w)
        sol.prim[i] .= uq_conserve_prim(sol.w[i], γ, uq)
    end

end

function uq_prim_conserve!(sol::AbstractSolution, γ::Real, uq::AbstractUQ)

    for i in eachindex(sol.prim)
        sol.w[i] .= uq_prim_conserve(sol.prim[i], γ, uq)
    end

end

"""
Calculate speed of sound

"""
function uq_sound_speed(prim::AbstractArray{<:AbstractFloat,2}, gamma::Real, uq::AbstractUQ) # single component

    if size(prim, 2) == uq.nm + 1 && uq.nm + 1 != uq.nq

        primRan = chaos_ran(prim, 2, uq)

        sosRan = zeros(uq.nq)
        for j in eachindex(sosRan)
            sosRan[j] = KitBase.sound_speed(primRan[end, j], gamma)
        end

        return ran_chaos(sosRan, uq)

    elseif size(prim, 2) == uq.nq

        sosRan = zeros(axes(prim, 2))
        for j in eachindex(sosRan)
            sosRan[j] = KitBase.sound_speed(prim[end, j], gamma)
        end

        return sosRan

    else

        throw("inconsistent random domain size in settings and solutions")

    end

end

#--- multiple component ---#
function uq_sound_speed(prim::AbstractArray{<:AbstractFloat,3}, gamma::Real, uq::AbstractUQ)

    if size(prim, 2) == uq.nm + 1 && uq.nm + 1 != uq.nq

        primRan = chaos_ran(prim, 2, uq)

        sosRan = zeros(uq.nq)
        for j in eachindex(sosRan)
            sosRan[j] = max(
                KitBase.sound_speed(primRan[end, j, 1], gamma),
                KitBase.sound_speed(primRan[end, j, 2], gamma),
            )
        end

        return ran_chaos(sosRan, uq)

    elseif size(prim, 2) == uq.nq

        sosRan = zeros(axes(prim, 2))
        for j in eachindex(sosRan)
            sosRan[j] = max(
                KitBase.sound_speed(prim[end, j, 1], gamma),
                KitBase.sound_speed(prim[end, j, 2], gamma),
            )
        end

        return sosRan

    else

        throw("inconsistent random domain size in settings and solutions")

    end

end

"""
Calculate collision time

"""
function uq_vhs_collision_time(
    prim::AbstractArray{<:Real,2},
    muRef::Real,
    omega::Real,
    uq::AbstractUQ,
) # deterministic viscosity

    if size(prim, 2) == uq.nm + 1 && uq.nm + 1 != uq.nq

        primRan = chaos_ran(prim, 2, uq)

        tauRan = zeros(uq.nq)
        for i in eachindex(tauRan)
            tauRan[i] = KitBase.vhs_collision_time(primRan[:, i], muRef, omega)
        end

        return ran_chaos(tauRan, uq)

    elseif size(prim, 2) == uq.nq

        tau = zeros(uq.nq)
        for i in eachindex(tau)
            tau[i] = KitBase.vhs_collision_time(prim[:, i], muRef, omega)
        end

        return tau

    else

        throw("inconsistent random domain size in settings and solutions")

    end

end

#--- stochastic viscosity ---#
function uq_vhs_collision_time(
    prim::Array{<:Real,2},
    muRef::Array{<:Real,1},
    omega::Real,
    uq::AbstractUQ,
)

    if size(prim, 2) == uq.nm + 1 && uq.nm + 1 != uq.nq

        primRan = chaos_ran(prim, 2, uq)
        muRan = chaos_ran(muRef, uq)

        tauRan = zeros(uq.nq)
        for i in eachindex(tauRan)
            tauRan[i] = KitBase.vhs_collision_time(primRan[:, i], muRan[i], omega)
        end

        return ran_chaos(tauRan, uq)

    elseif size(prim, 2) == uq.nq

        tau = zeros(uq.nq)
        for i in eachindex(tau)
            tau[i] = KitBase.vhs_collision_time(prim[:, i], muRef[i], omega)
        end

        return tau

    else

        throw("inconsistent random domain size in settings and solutions")

    end

end

function uq_vhs_collision_time(
    sol::AbstractSolution,
    muRef::Union{Real,Array{<:AbstractFloat,1}},
    omega::Real,
    uq::AbstractUQ,
)
    if ndims(sol.prim) == 1
        return [
            uq_vhs_collision_time(sol.prim[i], muRef, omega, uq) for
            i in eachindex(sol.prim)
        ]
    elseif ndims(sol.prim) == 2
        return [
            uq_vhs_collision_time(sol.prim[i, j], muRef, omega, uq) for
            i in axes(sol.prim, 1), j in axes(sol.prim, 2)
        ]
    end
end

"""
Calculate mixed collision time in AAP model

"""
function uq_aap_hs_collision_time(
    P::Array{<:AbstractFloat,3},
    mi::Real,
    ni::Real,
    me::Real,
    ne::Real,
    kn::Real,
    uq::AbstractUQ,
)

    if size(P, 2) == uq.nm + 1 && uq.nm + 1 != uq.nq

        prim = deepcopy(P[:, 1, :])
        τ = KitBase.aap_hs_collision_time(prim, mi, ni, me, ne, kn)

        return τ

    elseif size(P, 2) == uq.nq

        prim = deepcopy(P[:, end÷2+1, :])
        τ = KitBase.aap_hs_collision_time(prim, mi, ni, me, ne, kn)

        return τ

    else

        throw("inconsistent random domain size in settings and solutions")

    end

end

"""
Calculate pseudo primitive variables in AAP model

"""
function uq_aap_hs_prim(
    prim::Array{<:Real,3},
    tau::Array{<:Real,1},
    mi::Real,
    ni::Real,
    me::Real,
    ne::Real,
    kn::Real,
    uq::AbstractUQ,
)

    if size(prim, 2) == uq.nm + 1 && uq.nm + 1 != uq.nq

        primRan = chaos_ran(prim, 2, uq)

        # hard-sphere molecule
        mixPrimRan = similar(primRan)
        for j in axes(mixPrimRan, 2)
            mixPrimRan[:, j, :] .=
                KitBase.aap_hs_prim(primRan[:, j, :], tau, mi, ni, me, ne, kn)
        end

        mixPrimChaos = similar(prim)
        for k in axes(mixPrimChaos, 3)
            for i in axes(mixPrimChaos, 1)
                mixPrimChaos[i, :, k] .= ran_chaos(mixPrimRan[i, :, k], uq)
            end
        end

        return mixPrimChaos

    elseif size(prim, 2) == uq.nq

        mixPrimRan = similar(prim)
        for j in axes(mixPrimRan, 2)
            mixPrimRan[:, j, :] .=
                KitBase.aap_hs_prim(prim[:, j, :], tau, mi, ni, me, ne, kn)
        end

        return mixPrimRan

    else

        throw("inconsistent random domain size in settings and solutions")

    end

end


text:
@enum Interpolation Linear Nearest

struct Scaling{F, R}
    # a function to scale a value by, e.g. log10, sqrt etc
    scaling_function::F
    # If nothing, no scaling applied!
    range::R
end

Scaling() = Scaling(identity, nothing)

const NoScaling = Scaling{typeof(identity), Nothing}

struct Sampler{N, V} <: AbstractArray{RGBAf, 1}
    # the colors to sample from!
    colors::AbstractArray{T, N} where T
    # or an array of values, which are used to index into colors via interpolation!
    values::V
    # additional alpha that gets multiplied
    alpha::Float64
    # interpolation method for sampling
    interpolation::Interpolation
    scaling::Scaling
end

Base.size(sampler::Sampler) = size(sampler.values)

"""
    interpolated_getindex(cmap::AbstractArray, value::AbstractFloat, norm = (0.0, 1.0))

Like getindex, but accepts values between 0..1 and interpolates those to the full range.
You can use `norm`, to change the range of 0..1 to whatever you want.
"""
function interpolated_getindex(cmap::AbstractArray, value::Number, norm::VecTypes)
    cmin, cmax = norm
    cmin == cmax && error("Can't interpolate in a range where cmin == cmax. This can happen, for example, if a colorrange is set automatically but there's only one unique value present.")
    i01 = clamp((value - cmin) / (cmax - cmin), 0.0, 1.0)
    return interpolated_getindex(cmap, i01)
end

"""
    interpolated_getindex(cmap::AbstractArray, value::AbstractFloat)

Like getindex, but accepts values between 0..1 for `value` and interpolates those to the full range of `cmap`.
"""
function interpolated_getindex(cmap::AbstractArray{T}, i01::AbstractFloat) where T
    if !isfinite(i01)
        error("Looking up a non-finite or NaN value in a colormap is undefined.")
    end
    i1len = (i01 * (length(cmap) - 1)) + 1
    down = floor(Int, i1len)
    up = ceil(Int, i1len)
    down == up && return cmap[down]
    interp_val = i1len - down
    downc, upc = cmap[down], cmap[up]
    return convert(T, (downc * (1.0 - interp_val)) + (upc * interp_val))
end

function nearest_getindex(cmap::AbstractArray, value::AbstractFloat)
    idx = round(Int, i01 * (length(cmap) - 1)) + 1
    return cmap[idx]
end

"""
    apply_scaling(value::Number, scaling::Scaling)

Scales a number to the range 0..1.
"""
function apply_scaling(value::Number, scaling::NoScaling)::Float64
    return Float64(value)
end

function apply_scaling(value::Number, scaling::Scaling)::Float64
    value_scaled = scaling.scaling_function(value)
    if scaling.range === nothing
        return value_scaled
    end
    cmin, cmax = scaling.range
    clamped = clamp((value_scaled - cmin) / (cmax - cmin), 0.0, 1.0)
    if isfinite(clamped)
        return clamped
    else
        return 0.0
    end
end

function Base.getindex(sampler::Sampler, i)::RGBAf
    value = sampler.values[i]
    scaled = apply_scaling(value, sampler.scaling)
    c = if sampler.interpolation == Linear
        interpolated_getindex(sampler.colors, scaled)
    else
        nearest_getindex(sampler.colors, scaled)
    end
    return RGBAf(color(c), alpha(c) * sampler.alpha)
end

function Base.getindex(sampler::Sampler{2, <: AbstractVector{Vec2f}}, i)::RGBAf
    uv = sampler.values[i]
    colors = sampler.colors
    # indexing confirming to OpenGL uv indexing
    wsize = reverse(size(colors))
    wh = wsize .- 1
    x, y = round.(Int, Tuple(uv) .* wh) .+ 1
    c = colors[size(colors, 1) - (y - 1), x]
    return RGBAf(color(c), alpha(c) * sampler.alpha)
end

function sampler(cmap::Union{Symbol, String}, n::Int = 20;
                 scaling=Scaling(), alpha=1.0, interpolation=Linear)
    return sampler(cmap, LinRange(0, 1, n); scaling = scaling, alpha = alpha, interpolation = interpolation)
end

function sampler(cmap::Union{Symbol, String}, values::AbstractVector{<: AbstractFloat};
                 scaling=Scaling(), alpha=1.0, interpolation=Linear)

    cs = PlotUtils.get_colorscheme(cmap)

    colors = getindex.(Ref(cs), values)

    return Sampler(colors, values, alpha, interpolation, scaling)
end

function sampler(cmap::Vector{<: Colorant}, values::AbstractVector{<: AbstractFloat};
                 scaling=Scaling(), alpha=1.0, interpolation=Linear)
    return Sampler(RGBAf.(cmap), values, alpha, interpolation, scaling)
end

function sampler(cmap::AbstractVector, values, crange;
                 alpha=1.0, interpolation=Linear)
    return Sampler(to_color.(cmap), values, alpha, interpolation, Scaling(identity, crange))
end
# uv texture sampler
function sampler(cmap::Matrix{<: Colorant}, uv::AbstractVector{Vec2f};
                 alpha=1.0, interpolation=Linear)
    return Sampler(cmap, uv, alpha, interpolation, Scaling())
end


text:
using ObjectDetector, FileIO
using BenchmarkTools

yolomod = YOLO.v3_COCO(silent=false, w=416, h=416)

batch = emptybatch(yolomod)
img = load(joinpath(dirname(dirname(pathof(ObjectDetector))),"test","images","dog-cycle-car_nonsquare.png"))
img = load(joinpath(dirname(dirname(pathof(ObjectDetector))),"test","images","dog-cycle-car.png"))
batch[:,:,:,1], padding = prepareImage(img, yolomod)
save(joinpath(@__DIR__,"result.png"), cpu(batch[:,:,:,1]))

@time res = yolomod(batch, detectThresh=0.2, overlapThresh=0.8) # Run the model on the length-1 batch

imgBoxes = drawBoxes(img, yolomod, padding, res)
save(joinpath(@__DIR__,"result.png"), imgBoxes)

imgBoxes = drawBoxes(collect(img'), yolomod, padding, res, transpose=false)
save(joinpath(@__DIR__,"result_transposed.png"), imgBoxes)

ObjectDetector.benchmark()

using ImageCore, ObjectDetector, FileIO
img = ones(Gray, 200, 100)

yolomod = YOLO.v3_tiny_COCO(w=416, h=416, silent=true)
batch = emptybatch(yolomod)
batch[:,:,:,1], padding = prepareImage(img, yolomod)

res = collect([ padding[1] padding[2] 1.0-padding[3] 1.0-padding[4] 0.0 0.0;]') #note the transpose!

imgboxes = drawBoxes(img, yolomod, padding, res)
save(joinpath(@__DIR__, "test.png"), imgboxes)
all(imgboxes[1,:] .== Gray(0))
all(imgboxes[:,1] .== Gray(0))
all(imgboxes[end,:] .== Gray(0))
all(imgboxes[:,end] .== Gray(0))
imgboxes


text:
#-----------------------------------------------------------------------# Weight
@recipe function f(wt::Weight; nobs=50)
    xlabel --> "Number of Observations"
    ylabel --> "Weight Value"
    label --> name(wt)
    ylim --> (0, 1)
    linewidth --> 2
    [wt(i) for i in 1:nobs]
end

#-----------------------------------------------------------------------# Fallback
@recipe function f(o::OnlineStat)
    legend --> false
    axis --> false
    grid --> false
    ylim --> (0, 1)
    xlim --> (0, 1)
    annotations --> [(.5, .75, name(o) * ":"), (.5, .4, string(value(o)))]
    zeros(0)
end

#-----------------------------------------------------------------------# residual plot
@recipe function f(o::OnlineStat{VectorOb}, x::AbstractMatrix, y::AbstractVector)
    ylab --> "Residual"
    xlab --> "Observation Index"
    legend --> false
    @series begin
        linetype --> :scatter
        ŷ = predict(o, x)
        eachindex(y), y - ŷ
    end
end

@recipe function f(o::OnlineStat{XY})
    ylab --> "beta_j"
    xlab --> "j"
    seriestype --> :scatter
    coef(o)
end

@recipe function f(o::Series, x::AbstractMatrix, y::AbstractVector)
    for stat in o.stats
        @series begin stat end
    end
end

#-----------------------------------------------------------------------# StatCollection
@recipe function f(s::OnlineStatsBase.StatCollection)
    if :layout in keys(plotattributes)
        for stat in s.stats
            @series begin stat end
        end
    else  # hack to ensure series aren't sent to wrong subplots
        layout --> length(s.stats)
        for i in eachindex(s.stats)
            @series begin
                subplot --> i
                s.stats[i]
            end
        end
    end
end

#-----------------------------------------------------------------------# GroupBy
@recipe function f(o::GroupBy{T, <:Union{KHist, Hist}}) where {T}
    sort!(o.value)
    link --> :all
    for (k, v) in pairs(o.value)
        @series begin
            label --> k
            v
        end
    end
end
@recipe function f(o::GroupBy)
    sort!(o.value)
    link --> :all
    collect(keys(o.value)), value.(collect(values(o.value)))
end

#-----------------------------------------------------------------------# StatLag
@recipe function f(o::StatLag)
    layout --> length(o.value)
    for i in 1:length(o.value)
        @series begin
            subplot --> i
            label --> "Current - $(i-1)"
            o.value[i]
        end
    end
end

#-----------------------------------------------------------------------# AutoCov
@recipe function f(o::AutoCov)
    xlabel --> "Lag"
    ylabel --> "Autocorrelation"
    ylim --> (0, 1)
    seriestype --> :scatter
    autocor(o)
end

#-----------------------------------------------------------------------# CovMatrix
@recipe function f(o::CovMatrix; corr = false)
    seriestype --> :heatmap
    corr ? cor(o) : cov(o)
end

#-----------------------------------------------------------------------# Histograms
@recipe function f(o::HistogramStat; normed=true)
    e, c = edges(o), counts(o)
    inds = findfirst(x -> x > 0, c):findlast(x -> x > 0, c)
    closed = o.left ? :left : :right
    normed --> normed
    Histogram(e[vcat(inds, inds[end] + 1)], c[inds], closed)
end

@recipe function f(o::KHist; normed=true)
    x, y = value(o)
    y2 = normed ? y ./ area(o) : y
    # xlim --> extrema(o.ex)
    fillto --> 0
    alpha --> .5
    x, y2
end

#-----------------------------------------------------------------------# CountMap
@recipe function f(o::CountMap, kys = keys(o); sortby = :keys)
    seriestype --> :bar
    kys = collect(kys)
    vls = [o.value[ky] for ky in kys]
    sortby in [:keys, :values] || @warn("sortby = :$sortby not recognized")
    sp = sortby == :keys ? sortperm(kys) : sortperm(vls)
    x, y = string.(kys[sp]), vls[sp]
    hover --> ["($xi, $yi)" for (xi,yi) in zip(x, y)], :quiet
    label --> "Count"
    x, y
end

#-----------------------------------------------------------------------# Partition
@recipe f(o::AbstractPartition, fun=value) = o.parts, fun

@recipe function f(parts::Vector{Part{T, O}}, fun) where {T, O}
    sort!(parts)
    y = map(part -> fun(part.stat), parts)
    x = midpoint.(parts)
    # if parts[1].a isa Number
    #     xlim --> (parts[1].a, parts[end].b)
    # end
    if y[1] isa Number
        seriestype --> :step
        label --> name(parts[1].stat, false, false)
        x, y
    elseif y[1] isa VectorOb
        label --> name(parts[1].stat, false, false)
        y2 = plotshape(y)
        x2 = eltype(x) == Char ? string.(x) : x  # Plots can't handle Char
        if length(y[1]) == 2
            # seriestype --> :step
            fillto --> y2[:, 1]
            alpha --> .4
            linewidth --> 0
            x2, y2[:, 2]
        else
            x2, y2
        end
    elseif y[1] isa AbstractDict  # CountMap
        kys = []
        for item in y, ky in keys(item)
            ky ∉ kys && push!(kys, ky)
        end
        sort!(kys)
        @series begin
            label --> reshape(kys, (1, length(kys)))
            ylim --> (0, 1)
            linewidth --> .5
            seriestype --> :bar
            if parts[1].a isa Number
                bar_widths --> [p.b - p.a for p in parts]
            end
            y = plotshape(map(x -> reverse(cumsum(probs(x.stat, reverse(kys)))), parts))
            x, y
        end
    else
        error("No plot recipe exists for this kind of partition")
    end
end

plotshape(v::Vector) = v
plotshape(v::Vector{<:VectorOb}) = [v[i][j] for i in eachindex(v), j in eachindex(v[1])]

#---------------------------------------------------------------# [Indexed]Partition Hist
@recipe f(o::IndexedPartition{T,O}) where {T, O<:HistogramStat} = o.parts
@recipe f(o::Partition{T,O}) where {T, O<:HistogramStat} = o.parts

@recipe function f(parts::Vector{Part{T, O}}) where {T, O<:HistogramStat}
    sort!(parts)
    x = []
    y = []
    fillz = Int[]
    for part in parts
        edg = edges(part.stat)
        cnts = counts(part.stat)
        for i in eachindex(cnts)
            if cnts[i] > 0
                # rectangle
                push!(x, part.a); push!(y, edg[i])
                push!(x, part.a); push!(y, edg[i + 1])
                push!(x, part.b); push!(y, edg[i + 1])
                push!(x, part.b); push!(y, edg[i])
                push!(x, NaN); push!(y, NaN);
                # fill color
                push!(fillz, cnts[i])
            end
        end
    end
    @series begin
        seriestype := :shape
        linewidth --> 0
        linealpha --> 0
        legend --> false
        fillz := fillz
        x, y
    end
end

#-----------------------------------------------------------------------# NBClassifier
@recipe function f(o::NBClassifier)
    kys = collect(keys(o))
    layout --> nvars(o) + 1
    for j in 1:nvars(o)
        stats = o[j]
        for (i, s) in enumerate(stats)
            @series begin
                title --> "Var $j"
                legend --> false
                subplot --> j
                s
            end
        end
    end
    @series begin
        subplot --> nvars(o) + 1
        label --> reshape(kys, 1, length(kys))
        framestyle := :none
        zeros(0, length(kys))
    end
end


text:
# It seems like the individual flashes aren't bright enough to navigate.
# However, you might have a better option: the flashes seem to be
# synchronizing!
#
# In the example above, the first time all octopuses flash simultaneously is
# step 195:
#
# After step 193:
#
# 5877777777
# 8877777777
# 7777777777
# 7777777777
# 7777777777
# 7777777777
# 7777777777
# 7777777777
# 7777777777
# 7777777777
#
# After step 194:
# 6988888888
# 9988888888
# 8888888888
# 8888888888
# 8888888888
# 8888888888
# 8888888888
# 8888888888
# 8888888888
# 8888888888
#
# After step 195:
# 0000000000
# 0000000000
# 0000000000
# 0000000000
# 0000000000
# 0000000000
# 0000000000
# 0000000000
# 0000000000
# 0000000000
#
# If you can calculate the exact moments when the octopuses will all flash
# simultaneously, you should be able to navigate through the cavern. What is
# the first step during which all octopuses flash?

# load grid
grid = hcat(map(eachline("input.txt")) do l
  parse.(Int8, split(l, ""))
end...)

# simulate until they all flash
for i in Iterators.countfrom()
  # increment the energy of all octopuses
  grid .+= 1

  # find flashes
  all_flashes = falses(10, 10)
  flashes = grid .> 9
  while any(flashes)
    # increment the energy of the octopuses around flashing octopuses
    @views begin
      grid[1:9, 1:9] .+= flashes[2:10, 2:10]
      grid[1:10, 1:9] .+= flashes[1:10, 2:10]
      grid[2:10, 1:9] .+= flashes[1:9, 2:10]
      grid[2:10, 1:10] .+= flashes[1:9, 1:10]
      grid[2:10, 2:10] .+= flashes[1:9, 1:9]
      grid[1:10, 2:10] .+= flashes[1:10, 1:9]
      grid[1:9, 2:10] .+= flashes[2:10, 1:9]
      grid[1:9, 1:10] .+= flashes[2:10, 1:10]
    end

    # update all_flashes and see if there are any cascading flashes
    all_flashes = all_flashes .|| flashes
    flashes = grid .> 9 .&& .!all_flashes
  end

  # if they've all flashed, we're done
  if all(all_flashes)
    println("All flash at: ", i)
    break
  end

  # set any flashing octopuses back to zero
  grid[all_flashes] .= 0
end


text:
"Alias for `cond`"
|ₚ(x, y) = cond(x, y)

"""
$(SIGNATURES)

Convenience function to condition random variable `x` with a predicate `p`

`condp(x, p) = cond(x, p(x))`

```jldoctest
a = condp(poisson(0.5), iseven)
b = a |ₚ (x -> x < 10)
```

"""
condp(x, p) = cond(x, lift(p)(x))

# "Alias for `condp`"
# |ₚ(x, p) = condp(x, p)


text:
using GridArrays

using Test, LinearAlgebra, DomainSets, Plots, StaticArrays

using GridArrays: MaskedGrid, IndexSubGrid, randompoint,
    component, components, productgrid

using GridArrays.Test: test_generic_grid, test_interval_grid,
    grid_iterator1, grid_iterator2

function delimit(s::AbstractString)
    println()
    println("## ",s)
end

types = (Float64,BigFloat)

include("test_grids.jl")
include("test_gauss.jl")
include("test_periodiccartesianindices.jl")
include("test_boundingbox.jl")
include("test_gauss.jl")
include("test_subgrids.jl")
include("test_randomgrid.jl")
include("test_plots.jl")


text:
using VortexDistributions

# simple test field
Nx = 400; Ny = 400
Lx = 200; Ly = Lx
x = LinRange(-Lx/2,Lx/2, Nx+1)[1:end-1]; y = LinRange(-Ly/2,Ly/2, Ny+1)[1:end-1];
psi0 = one.(x*y') |> complex

psi = Sphere(copy(psi0),x,y)

xp = 1.133
yp = 2.787
vp = PointVortex(xp,yp,1)
sp = ScalarVortex(vp)

vortex!(psi, sp)

vfound = findvortices(psi)

dx = Δ(psi.x)

@test isapprox(vfound[1].xv, xp, rtol=dx/4)
@test isapprox(vfound[1].yv, yp, rtol=dx/4)
@test vfound[1].qv == 1


text:
"""
    Dirac(x)

A *Dirac distribution* is parameterized by its only value `x`, and takes its value with probability 1.

```math
P(X = \\hat{x}) = \\begin{cases}
1 & \\quad \\text{for } \\hat{x} = x, \\\\
0 & \\quad \\text{for } \\hat{x} \\neq x.
\\end{cases}
```

```julia
Dirac(2.5)   # Dirac distribution with value x = 2.5
```

External links:

* [Dirac measure on Wikipedia](http://en.wikipedia.org/wiki/Dirac_measure)
"""
struct Dirac{T} <: DiscreteUnivariateDistribution
    value::T
end

Base.eltype(::Type{Dirac{T}}) where {T} = T

insupport(d::Dirac, x::Real) = x == d.value
minimum(d::Dirac) = d.value
maximum(d::Dirac) = d.value
support(d::Dirac) = (d.value,)

#### Properties
mean(d::Dirac) = d.value
var(d::Dirac{T}) where {T} = zero(T)

mode(d::Dirac) = d.value

entropy(d::Dirac{T}) where {T} = zero(T)

#### Evaluation

pdf(d::Dirac, x::Real) = insupport(d, x) ? 1.0 : 0.0
logpdf(d::Dirac, x::Real) = insupport(d, x) ? 0.0 : -Inf

cdf(d::Dirac, x::Real) = x < d.value ? 0.0 : 1.0
cdf(d::Dirac, x::Integer) = x < d.value ? 0.0 : 1.0

quantile(d::Dirac{T}, p::Real) where {T} = 0 <= p <= 1 ? d.value : T(NaN)

mgf(d::Dirac, t) = exp(t * d.value)
cf(d::Dirac, t) = cis(t * d.value)

#### Sampling

rand(rng::AbstractRNG, d::Dirac) = d.value


text:
module Braiders

#using ..Crypto

import SynchronicBallot
using SynchronicBallot: SocketConfig
using DemeNet: DemeSpec, Notary, Cypher, Signer, Deme, Contract, Certificate, ID, DemeID, DHsym, DHasym
import PeaceVote.braid!
import ..Types: Braid ### The Braid could be defined here

using Pkg.TOML

function Braid(metadata::Vector{UInt8},ballot::Array{UInt8,2})

    ids = ID[]

    N = size(ballot,2)

    for i in 1:N
        id = ID(ballot[:,i],base=16)
        push!(ids,id)
    end

    braid = Braid(nothing,nothing,ids)

    return braid
end

function validate(braid::Braid,id::ID,deme::Deme)
    return id in braid.ids
end

function signbraid(metadata::Vector{UInt8},ballot::Array{UInt8,2},validate::Function,signer::Signer)

    braid = Braid(metadata,ballot)

    @assert validate(braid)

    cert = Certificate(braid,signer)

    sdict = Dict(cert.signature)
    io = IOBuffer()
    TOML.print(io,sdict)
    return take!(io)
end

function Mixer(port,cypher::Cypher,notary::Notary,signer::Signer)
    mixergate = SocketConfig(nothing,DHasym(cypher,notary,signer),cypher.secureio)
    mixermember = SocketConfig(nothing,DHasym(cypher,notary,signer),cypher.secureio)

    return SynchronicBallot.Mixer(port,mixergate,mixermember) # should be renamed to Mixer
end

Mixer(port,deme::Deme,signer::Signer) = Mixer(port,deme.cypher,deme.notary,signer)

struct BraiderConfig{T}
    port::T # braiderport
    ballotport::T # mixerport
    N::UInt8
    M::UInt8
    gateid::ID # braiderid
    mixerid::DemeID
end

struct Braider
    server
    voters::Set{ID}
end

import Base.take!

### Perhaps I also need to put metadata with the ballots channel.
### Because otherwise we may loose a track (a simple bugfix)
function take!(braider::Braider)
    metadata,ballot,signaturesbin = take!(braider.server.ballots)
    braid = Braid(UInt8[],ballot)
    signatures = Dict{String,Any}[TOML.parse(String(i)) for i in signaturesbin]
    contract = Contract(braid,signatures)
    return contract
end

function Braider(braider::BraiderConfig,deme::Deme,mixerdeme::Deme,signer::Signer)

    voters = Set{ID}()

    gatemixer = SocketConfig(braider.mixerid.id,DHasym(mixerdeme.cypher,mixerdeme.notary),mixerdeme.cypher.secureio)
    gatemember =  SocketConfig(voters,DHsym(deme.cypher,deme.notary,signer),deme.cypher.secureio)

    server = SynchronicBallot.GateKeeper(braider.port,braider.ballotport,braider.N,braider.M,gatemixer,gatemember,()->UInt8[]) #() -> Vector{UInt8}(""))
    Braider(server,voters)
end

# New in this context just overloads generated function since it is not necessary.
function Braider(braider::BraiderConfig,deme::Deme,signer::Signer)
    mixeruuid = braider.mixerid.uuid

    mixerdemespec = DemeSpec(mixeruuid)
    mixerdeme = Deme(mixerdemespec)

    Braider(braider,deme,mixerdeme,signer)
end

function Braider(deme::Deme,signer::Signer)
    systemconfig = SystemConfig(deme) ### This is where the config file is verified
    braider = systemconfig.braider
    Braider(braider,deme,signer)
end

function braid!(config::BraiderConfig,deme::Deme,mixerdeme::Deme,voter::Signer,signer::Signer)

    membergate = SocketConfig(config.gateid,DHsym(deme.cypher,deme.notary,signer),deme.cypher.secureio)
    membermixer = SocketConfig(config.mixerid.id,DHasym(mixerdeme.cypher,mixerdeme.notary),deme.cypher.secureio)

    ### I need to take out a validate function, to check the braid formed by the server
    # (2) The hash of the full ledger
    # (3) Some other metadata. Whether the same port was used by everyone to form the braid as well the same mixer and server.

    idbytes = Vector{UInt8}(voter.id,base=16,length=config.M)

    SynchronicBallot.vote(config.port,membergate,membermixer,idbytes,(m,b)->signbraid(m,b,braid->validate(braid,voter.id,deme),signer))
end

function braid!(config::BraiderConfig,deme::Deme,voter::Signer,signer::Signer)
    mixeruuid = config.mixerid.uuid

    mixerdemespec = DemeSpec(mixeruuid)
    mixerdeme = Deme(mixerdemespec)

    braid!(config,deme,mixerdeme,voter,signer)
end

export Mixer, Braider, braid!

end # module


text:
using Test
last_runtimes = Float64[]

#=

	Computing the Fibonacci sequence
	===

	Below, you will find a function `fib` which computes the first `N` (default: 100.000) numbers of the Fibonacci sequence

	> Your mission, should you decide to accept it, is to make `fib` (roughly) 10× faster. 🚀

	There are at least two mistakes in the given code. (One much more relevant than the other).
	If you want, you can find some tips in the hint boxes below.

	Useful functions for `Vector`s might be:
	- `sizehint!()` gives Julia the hint how large the vector might become
	- `zeros(Int64, N)`  creates a vector of zeros with element type `Int64`
	- `Vector{Int64}(undef, N)` creates an uninitialized vector with element type `Int64`

	Check the documentation (e.g. `?sizehint!`) for their function.
=#

function fib(N = 100_000; x₀ = 1, x₁ = 1)
	# modify this code to improve the performance

	x = Int64[]
	append!(x, (x₀, x₁))
	for i in 3:N
		push!(x, x[i-1] + x[i-2])
	end
	return x
end

@test fib()[end] == 2754320626097736315  # this line tests if `fib` works correctly

# this is how you call the function:
x = fib()

#=

How to benchmark a function:

  To test the performance, you can either use `@time`\`@elapsed` or,
  the accurate `@btime`\`@belapsed` from the package `BenchmarkTools`

=#

@time fib()

begin
	# run function again and show hisotry
	current_runtime = @elapsed fib()
	push!(last_runtimes, current_runtime)

	println("Current runtime: ", 1000 * current_runtime, " ms")
	println("Last runtimes: ")
	for t in reverse(last_runtimes)
		println("  ", 1000 * t, " ms")
	end
end

println()

#=

	If you want some hints, continue reading:

	Hint #1"
    ---

	Adding new element to an existing `Vector` can take up much time.
	It is better to allocate enough memory beforehand.

    You could use `sizehint!(x, N)`, which tells the program which size to allocate for that array.

	Or you could allocate directly a vector of size `N`, but then you need to modify the for loop accordingly, since we then don't have to `push!` element into the vector, but just to access the `i`th position.

	Hint #2
	---

	The main problem is this: The vector `x` has type `typeof(x) == Vector{Any}`!
    This wastes a lot of resources.

    To avoid this, you could at least create it as
	`x = Int64[]`.

    Or in combination with Hint #1, you could use
	`x = zeros(Int64, N)` and change the `for` loop accordingly.

=#

max_stars_repo_name: SteffenPL/Julia-for-mathematicians

max_issues_repo_name: SteffenPL/Julia-for-mathematicians

max_forks_repo_name: SteffenPL/Julia-for-mathematicians


text:
function  is_prime(n::Int64)
	if(n<=1)  return false
	end
	if(n==2)  return true
	end
     if( n%2==0)  return false
     end
	i::Int64=3
	while(i*i<=n)
		if( n%i==0)  return false
		end
		i+=2
	end

	return true;

end
"""
# test
for i=1:100
	if( prime(i))
		print(i,"  ")
	end
end
"""

function prime(N::Int64,a::Array{Int64,1})
	  b=fill(true,N);
	  b[1]=false;
	  for i=2:N
	  	for j=i*i:i:N
	  		if(b[i])
	  			b[j]=false
	  		end
	  	end
	  end

	  for i=1:N
	  	if (b[i])
	  		push!(a,i)
	  	end
	  end
end

"""
# test
a=Int64[]
prime(100,a)

print(a)

"""

function prime_table(n::Int)
	b=fill(true,N);
	b[1]=false;
	for i=2:N
	 for j=i*i:i:N
	 	if(b[i])
	 		b[j]=false
	  	end
	  end
	end
	return b
end


text:
## The cells in this file can only be run after running those in drycbl_run_plot.jl

## Transfer the state from the high to the low resolution and plot it.
MicroHH.transfer_state!(f1, f2, g1, g2)
node1[] = s_bot1 .- mean(s_bot1)
node2[] = s_bot2 .- mean(s_bot2)

## Run the model for another 900 sec.
m.timeloop[1].end_time += 1800.
m.timeloop[1].iend_time += 1800 * MicroHH.ifactor
in_progress = true
while in_progress
    global in_progress = step_model!(m)
    node11[] = s_bot1 .- mean(s_bot1)
    node21[] = s_bot2 .- mean(s_bot2)
    node12[] = s1[:, 1, :] .- mean(s1, dims=(1, 2))[:, 1, :]
    node22[] = s2[:, 1, :] .- mean(s2, dims=(1, 2))[:, 1, :]
end


text:
fs2()

temp_config = joinpath(F.PATHS[:folder], "config.md")
write(temp_config, "@def author = \"Stefan Zweig\"\n")
temp_index = joinpath(F.PATHS[:folder], "index.md")
write(temp_index, "blah blah")
temp_index2 = joinpath(F.PATHS[:folder], "index.html")
write(temp_index2, "blah blih")
temp_blah = joinpath(F.PATHS[:folder], "blah.md")
write(temp_blah, "blah blah")
temp_html = joinpath(F.PATHS[:folder], "temp.html")
write(temp_html, "some html")
temp_rnd = joinpath(F.PATHS[:folder], "temp.rnd")
write(temp_rnd, "some random")
temp_css = joinpath(F.PATHS[:css], "temp.css")
write(temp_css, "some css")
temp_js = joinpath(F.PATHS[:libs], "foo.js")
write(temp_js, "some js")

F.process_config()

@testset "Prep outdir" begin
    F.prepare_output_dir()
    @test isdir(F.PATHS[:site])
    temp_out = joinpath(F.PATHS[:site], "tmp.html")
    write(temp_out, "This is a test page.\n")
    # clear is false => file should remain
    F.prepare_output_dir(false)
    @test isfile(temp_out)
    # clear is true => file should go
    F.prepare_output_dir(true)
    @test !isfile(temp_out)
end

@testset "Scan dir" begin
    println("🐝 Testing file tracking...:")
    # it also tests add_if_new_file and last
    md_files       = Dict{Pair{String, String}, Float64}()
    html_files     = empty(md_files)
    other_files    = empty(md_files)
    infra_files    = empty(md_files)
    literate_files = empty(md_files)
    watched_files  = [other_files, infra_files, md_files, html_files, literate_files]
    F.scan_input_dir!(other_files, infra_files, md_files, html_files, literate_files, true)
    @test haskey(md_files, F.PATHS[:folder]=>"blah.md")
    @test md_files[F.PATHS[:folder]=>"blah.md"] == mtime(temp_blah) == stat(temp_blah).mtime
    @test html_files[F.PATHS[:folder]=>"temp.html"] == mtime(temp_html)
    @test other_files[F.PATHS[:folder]=>"temp.rnd"] == mtime(temp_rnd)
end

@testset "Config+write" begin
    F.process_config()
    @test F.GLOBAL_VARS["author"].first == "Stefan Zweig"
    rm(temp_config)
    @test_logs (:warn, "I didn't find a config file. Ignoring.") F.process_config()
    # testing write
    head = "head"
    pg_foot = "\npage_foot"
    foot = "foot {{if hasmath}} {{fill author}}{{end}}"

    out_file = F.form_output_path(F.PATHS[:site], "index.html", :html)
    F.write_page(F.PATHS[:folder], "index.md", head, pg_foot, foot, out_file)

    @test isfile(out_file)
    @test read(out_file, String) == "head\n
\n
blah blah
\n\n\npage_foot\n
\nfoot  Stefan Zweig"
end

temp_config = joinpath(F.PATHS[:folder], "config.md")
write(temp_config, "@def author = \"Stefan Zweig\"\n")
rm(temp_index2)

@testset "Part convert" begin
    write(joinpath(F.PATHS[:layout], "head.html"), raw"""
        <!doctype html>

                <meta charset="UTF-8">
                <link rel="stylesheet" href="/css/main.css">

        """)
    write(joinpath(F.PATHS[:layout], "page_foot.html"), raw"""

                        © All rights reserved.

""")
    write(joinpath(F.PATHS[:layout], "foot.html"), raw"""

        """)

    clear = true

    watched_files = F.fd_setup(; clear=clear)

    F.fd_fullpass(watched_files; clear=clear)

    @test issubset(["css", "libs", "index.html"], readdir(F.PATHS[:site]))
    @test issubset(["temp", "temp.rnd"], readdir(F.PATHS[:site]))
    @test issubset(["index.html"], readdir(joinpath(F.PATHS[:site], "temp")))
    @test issubset(["index.html"], readdir(joinpath(F.PATHS[:site], "blah")))

    @test all(split(read(joinpath(F.PATHS[:site], "index.html"), String)) .== split("<!doctype html>\n\n\t\n\t\t<meta charset=\"UTF-8\">\n\t\t<link rel=\"stylesheet\" href=\"/css/main.css\">\n\t\n\n
\n
blah blah
\n\n
\n\t\t
\n\t\t\t\t© All rights reserved.\n\t\t
\n
\n
\n    \n"))
end

@testset "Err procfile" begin
    write(temp_index, "blah blah { blih etc")
    println("🐝 Testing error message...:")
    @test_throws F.OCBlockError F.process_file_err(:md, F.PATHS[:folder] => "index.md"; clear=false)
end


text:
module PkgVersionHelper

export upcheck

using Pkg

# from https://stackoverflow.com/questions/62667741/julia-is-there-a-way-to-find-the-latest-possible-version-number-of-a-package
function max_ver(pkgname::AbstractString)
    path = joinpath(DEPOT_PATH[1], "registries", "General",
                    first(pkgname, 1), pkgname, "Versions.toml")
    maximum(VersionNumber.(keys(Pkg.TOML.parse(read(path, String)))))
end

"""
    upcheck(; indirect_deps=false)

Checks which package is not up to date. It wil return a
`Dict{String, Tuple{VersionNumber, VersionNumber}}` with these key-value paiars

    `PkgName => (installed_version, latest_version)`

Keyword argument `indirect_deps` if `true` includes the indirect dependencies.
"""
function upcheck(; indirect_deps = false)
    deps = Pkg.dependencies()
    installs = Dict{String, Tuple{VersionNumber, VersionNumber}}()
    for (uuid, dep) in deps
        # is it direct dependency?
        dep.is_direct_dep | indirect_deps || continue

        # is it a standard library
        Pkg.Types.is_stdlib(uuid) && continue

        try
            if max_ver(dep.name) <= dep.version
                continue
            end
            installs[dep.name] = (dep.version, max_ver(dep.name))
        catch e
            println("$(dep.name) may not be a registered package")
        end
    end
    return installs
end

end