text:
using MathProgBase, Ipopt, StatsBase, PyPlot, LaTeXStrings, Base.Iterators.cycle, Clp
PyPlot.matplotlib[:rcParams]["figure.autolayout"] = "True"
srand(0)
PyPlot.matplotlib[:rc]("font", family="serif", serif="Times", size=16)
global_width = 3
type RFTL
    x
    eta
    projection
end
function update(rftl, gradient)
    rftl.x = rftl.projection(rftl.x + rftl.eta * gradient)
end
function get_vector(rftl)
    return rftl.x
end
function play_sound()
    if is_apple()
        run(`say "your program has finished"`)
    end
end
function my_show()
    if is_apple()
        show()
    end
end
function gibbs_sampler(A, b, x)
    n = size(A)[2]
    for iter in 1:n * 3
        j = rand(1:n)
        x_copy = x
        x_copy[j] = 0
        max_coord = minimum((b - A * x_copy) ./ A[:, j])
        max_coord = min(max_coord, 1)
        x[j] = rand() * max_coord
    end
    return x
end
function meta_frank_wolfe(dim_lam, X_list, projection, f, gradient, eta, T, K)
    x0 = zeros(dim_lam)
    rftl = [RFTL(x0, eta, projection) for k in 1:K]
    reward = []
    for t in 1:T
        X = X_list[t]
        v = [get_vector(rftl[k]) for k in 1:K]
        x = zeros(K + 1, dim_lam)
        for k in 1:K
            update(rftl[k], gradient(x[k, :], X))
            x[k + 1, :] = x[k, :] + v[k] / K
        end
        push!(reward, f(x[K + 1, :], X))
    end
    reward = cumsum(reward)
    return reward
end
function frank_wolfe(dim_lam, X, A, b, f, gradient, K)
    x = zeros(dim_lam)
    for k in 1:K
        grad_x = gradient(x, X)
        sol = linprog(-grad_x, A, '<', b, 0.0, 1.0, ClpSolver())
        if sol.status == :Optimal
            v = sol.sol
        else
            error("No solution was found.")
        end
        x += v / K
    end
    return x
end
function frank_wolfe2(dim_lam, A, b, gradient, K)
    x = zeros(dim_lam)
    for k in 1:K
        grad_x = gradient(x)
        sol = linprog(-grad_x, A, '<', b, 0.0, 1.0, ClpSolver())
        if sol.status == :Optimal
            v = sol.sol
        else
            error("No solution was found.")
        end
        x += v / K
    end
    return x
end
function online_gradient_ascent(dim_lam, X_list, projection, f, gradient, eta, T)
    x0 = zeros(dim_lam)
    rftl = RFTL(x0, eta, projection)
    reward = []
    for t in 1:T
        X = X_list[t]
        v = get_vector(rftl)
        push!(reward, f(v, X))
        update(rftl, gradient(v, X))
    end
    reward = cumsum(reward)
    return reward
end
function online_random(dim_lam, X_list, f, linear_projection, T)
    reward = []
    for t in 1:T
        X = X_list[t]
        v = rand(dim_lam)
        v = linear_projection(v)
        push!(reward, f(v, X))
    end
    reward = cumsum(reward)
    return reward
end
function online_random_gibbs(dim_lam, X_list, f, A, b, T)
    reward = zeros(T)
    gibbs_v = zeros(size(A)[2])
    for t in 1:T
        X = X_list[t]
        for j in 1:100
            gibbs_v = gibbs_sampler(A, b, gibbs_v)
            reward[t] = max(reward[t], f(gibbs_v, X))
        end
    end
    reward = cumsum(reward)
    return reward
end
function generate_projection_function(A, b)
    n = size(A)[2]
    function projection(x0)
        sol = quadprog(-x0, eye(n), A, '<', b, 0.0, 1.0, IpoptSolver(print_level=0))
        if sol.status == :Optimal
            return sol.sol
        end
        error("No solution was found.")
    end
    return projection
end

function offline(X_list, f, gradient, A, b, K)
    T = length(X_list)
    reward = zeros(T)
    gradient_cumul = x->sum([gradient(x, params) for params in X_list])
    v = frank_wolfe2(size(A)[2], A, b, gradient_cumul, K)
    for i in 1:T
        reward[i] = f(v, X_list[i])
    end
    return cumsum(reward)
end

function experiment_coverage(size_U, n_sets, T, K)
    function generate_a_param()
        B = []
        for i in 1:n_sets
            lo = round(Int, size_U * 0.2)
            hi = round(Int, size_U * 0.8)
            set_size = rand(lo:hi)
            push!(B, IntSet(sample(1:size_U, set_size, replace=false)))
        end
        w = rand(size_U) * 100
        return (B, w)
    end
    X_list = [generate_a_param() for i in 1:T]
    function f_cover(x, params)
        B = params[1]
        w = params[2]
        res = 0
        for u in 1:size_U
            prod = 1
            for i in 1:n_sets
                if u in B[i]
                    prod *= (1 - x[i])
                end
            end
            res += w[u] * (1 - prod)
        end
        return res
    end
    function concave_f(x)
        res = 0
        for u in 1:size_U
            res += w[u] * min(1, sum([x[i] for i in xrange(n_sets) if u in B[i]]))
        end
        return res
    end
    function partial(x, i, B, w)
        res = 0
        for u in 1:size_U
            prod = 1
            for j in 1:n_sets
                if j != i && u in B[j]
                    prod *= (1 - x[j])
                end
            end
            res += w[u] * prod
        end
        return res
    end
    function gradient_cover(x, params)
        B = params[1]
        w = params[2]
        res = zeros(n_sets)
        for i in 1:n_sets
            res[i] = partial(x, i, B, w)
        end
        return res
    end
    function stochastic_gradient_cover(x, params)
        B = params[1]
        w = params[2]
        function func(S)
            if isempty(S)
                return 0
            end
            all_sets = [B[i] for i in S]
            union_set = reduce(union, all_sets)
            return sum([w[i] for i in union_set])
        end
        function gradient_i(i)
            R = [j for j in 1:n_sets if j != i && rand() < x[j]]
            R_i = vcat(R, [i])
            return func(R_i) - func(R)
        end
        grad = zeros(n_sets)
        for i in 1:n_sets
            grad[i] = gradient_i(i)
        end
        return grad
    end
    function gradient_concave(x, params)
        B = params[1]
        w = params[2]
        function partial_concave(i)
            res = 0
            for u in 1:size_U
                if !(u in B[i])
                    continue
                end
                if sum([x[j] for j in 1:n_sets if u in B[j]]) <= 1
                    res += w[u]
                end
            end
            return res
        end
        grad = zeros(n_sets)
        for i in 1:n_sets
            grad[i] = partial_concave(i)
        end
        return grad
    end
    n_constraints = 2
    A = rand(n_constraints, n_sets)
    b = 1.0
    linear_projection = generate_projection_function(A, b)
    reward_offline = offline(X_list, f_cover, gradient_cover, A, b, K)
    println(reward_offline)
    function generate_figure_regret()
        reward_list = []
        push!(reward_list, reward_offline - meta_frank_wolfe(n_sets, X_list,
                                                    linear_projection, f_cover, gradient_cover, 0.1, T, K))
        push!(reward_list, reward_offline - meta_frank_wolfe(n_sets, X_list,
                                                    linear_projection, f_cover, gradient_cover, 0.01, T, K))
        push!(reward_list, reward_offline - online_gradient_ascent(n_sets, X_list,
                                                linear_projection, f_cover, gradient_cover, 0.1, T))
        push!(reward_list, reward_offline - online_gradient_ascent(n_sets, X_list,
                                                linear_projection, f_cover, gradient_cover, 0.01, T))
        push!(reward_list, reward_offline - online_gradient_ascent(n_sets, X_list,
                                                linear_projection, f_cover, gradient_concave, 0.1, T))
        push!(reward_list, reward_offline - online_gradient_ascent(n_sets, X_list,
                                                linear_projection, f_cover, gradient_concave, 0.01, T))
        push!(reward_list, reward_offline - online_random_gibbs(n_sets, X_list, f_cover, A, b, T))
        labels = [L"FW ($\eta=0.1$)", L"FW ($\eta=0.01$)", L"GA ($\eta=0.1$)", L"GA ($\eta=0.01$)", L"SurrGA ($\eta=0.1$)", L"SurrGA ($\eta=0.01$)", "OfflApprMax"]
        marker = cycle((".", ",", "+", "_", "o", "x", "*"))
        linecycler = cycle(("-", "--", "-.", ":"))
        println(reward_list)
        for zipped in zip(reward_list, labels, marker, linecycler)
            reward, label, marker_iter, line_iter = zipped
            plot(1:T, reward, label=label, linestyle=line_iter, linewidth=global_width)
        end
        ticklabel_format(style="sci", axis="y", scilimits=(0, 0))
        legend(loc="best")
        xlabel("Iteration index")
        ylabel(L"$(1-1/e)$-Regret")
        grid()
        # savefig("coverage_regret_T2.pdf")
        my_show()
    end
    function generate_figure_regret_stochastic()
        reward_list = []
        push!(reward_list, reward_offline - meta_frank_wolfe(n_sets, X_list,
                                                    linear_projection, f_cover, stochastic_gradient_cover, 0.1, T, K))
        push!(reward_list, reward_offline - meta_frank_wolfe(n_sets, X_list,
                                                    linear_projection, f_cover, stochastic_gradient_cover, 0.01, T, K))
        push!(reward_list, reward_offline - online_gradient_ascent(n_sets, X_list,
                                                linear_projection, f_cover, stochastic_gradient_cover, 0.1, T))
        push!(reward_list, reward_offline - online_gradient_ascent(n_sets, X_list,
                                                linear_projection, f_cover, stochastic_gradient_cover, 0.01, T))
        push!(reward_list, reward_offline - online_gradient_ascent(n_sets, X_list,
                                                linear_projection, f_cover, gradient_concave, 0.1, T))
        push!(reward_list, reward_offline - online_gradient_ascent(n_sets, X_list,
                                                linear_projection, f_cover, gradient_concave, 0.01, T))
        push!(reward_list, reward_offline - online_random_gibbs(n_sets, X_list, f_cover, A, b, T))
        labels = [L"FW ($\eta=0.1$)", L"FW ($\eta=0.01$)", L"GA ($\eta=0.1$)", L"GA ($\eta=0.01$)", L"SurrGA ($\eta=0.1$)", L"SurrGA ($\eta=0.01$)", "OfflApprMax"]
        marker = cycle((".", ",", "+", "_", "o", "x", "*"))
        linecycler = cycle(("-", "--", "-.", ":"))
        println(reward_list)
        for zipped in zip(reward_list, labels, marker, linecycler)
            reward, label, marker_iter, line_iter = zipped
            plot(1:T, reward, label=label, linestyle=line_iter, linewidth=global_width)
        end
        ticklabel_format(style="sci", axis="y", scilimits=(0, 0))
        legend(loc="best")
        xlabel("Iteration index")
        ylabel(L"$(1-1/e)$-Regret")
        grid()
        # savefig("stoch_coverage_regret_T2.pdf")
        my_show()
    end
    # generate_figure_regret()
    generate_figure_regret_stochastic()
end
function experiment_nqp(dim_x, T, K)
    function generate_a_param(dim_x, alpha=100)
        return rand(dim_x, dim_x) * (-alpha)
    end
    n_constraints = 2
    A = rand(n_constraints, dim_x)
    b = 1.0
    linear_projection = generate_projection_function(A, b)
    function f_nqp(x, H)
        return ((x / 2 - 1)' * H * x)
    end
    function gradient_nqp(x, H)
        return (H + H') / 2 * x - H' * ones(dim_x)
    end
    function generate_figure_regret()
        all_reward_lists = []
        for experiment_iter in 1:1
            X_list = [generate_a_param(dim_x) for i in 1:T]
            reward_offline = offline(X_list, f_nqp, gradient_nqp, A, b, K)
            reward_list = []
            push!(reward_list, reward_offline - meta_frank_wolfe(dim_x, X_list,
                                                        linear_projection, f_nqp, gradient_nqp, 0.001, T, K))
            push!(reward_list, reward_offline - meta_frank_wolfe(dim_x, X_list,
                                                        linear_projection, f_nqp, gradient_nqp, 0.0008, T, K))
            push!(reward_list, reward_offline - online_gradient_ascent(dim_x, X_list,
                                                    linear_projection, f_nqp, gradient_nqp, 0.001, T))
            push!(reward_list, reward_offline - online_gradient_ascent(dim_x, X_list,
                                                    linear_projection, f_nqp, gradient_nqp, 0.0008, T))
            push!(all_reward_lists, reward_list)
        end
        all_reward_lists = mean(all_reward_lists)
        labels = [L"FW ($\eta=0.1$)", L"FW ($\eta=0.01$)", L"GA ($\eta=0.1$)", L"GA ($\eta=0.01$)"]
        marker = cycle((".", ",", "+", "_", "o", "x", "*"))
        linecycler = cycle(("-", "--", "-.", ":"))
        println(all_reward_lists)
        for zipped in zip(all_reward_lists, labels, marker, linecycler)
            reward, label, marker_iter, line_iter = zipped
            plot(1:T, reward, label=label,  linestyle=line_iter, linewidth=global_width)
        end
        ticklabel_format(style="sci", axis="y", scilimits=(0, 0))
        legend(loc="best")
        xlabel("Iteration index")
        ylabel(L"$(1-1/e)$-Regret")
        grid()
        play_sound()
        my_show()
    end
    function generate_figure_regret_dim(alpha)
        dim_x_list = 5:5:50
        all_reward_lists = []
        reward_list = zeros(4, size(dim_x_list)[1])
        for experiment_iter in 1:5
            for (idx, dim_x) in enumerate(dim_x_list)
                tic()
                println("parameter\t", idx, ", ", dim_x)
                X_list = [generate_a_param(dim_x, alpha) for i in 1:T]
                A = rand(n_constraints, dim_x)
                b = 1.0
                linear_projection = generate_projection_function(A, b)
                X_list_sum = sum(X_list)
                v = frank_wolfe(dim_x, X_list_sum, A, b, f_nqp, gradient_nqp, K)
                reward_offline = f_nqp(v, X_list_sum)
                reward_list[:, idx] =  reward_offline - [meta_frank_wolfe(dim_x, X_list, linear_projection, f_nqp, gradient_nqp, 0.001, T, K)[end],
                    meta_frank_wolfe(dim_x, X_list, linear_projection, f_nqp, gradient_nqp, 0.0008, T, K)[end],
                    online_gradient_ascent(dim_x, X_list, linear_projection, f_nqp, gradient_nqp, 0.001, T)[end],
                    online_gradient_ascent(dim_x, X_list, linear_projection, f_nqp, gradient_nqp, 0.0008, T)[end]]
                toc()
            end
            push!(all_reward_lists, reward_list)
        end
        reward_list = mean(all_reward_lists)
        reward_list = [reward_list[i, :] for i in 1:size(reward_list)[1]]
        labels = [L"FW ($\eta=0.1$)", L"FW ($\eta=0.01$)", L"GA ($\eta=0.1$)", L"GA ($\eta=0.01$)"]
        marker = cycle((".", ",", "+", "_", "o", "x", "*"))
        linecycler = cycle(("-", "--", "-.", ":"))
        print(reward_list)
        for zipped in zip(reward_list, labels, marker, linecycler)
            reward, label, marker_iter, line_iter = zipped
            plot(dim_x_list, reward, label=label, linestyle=line_iter, linewidth=global_width)
        end
        legend(loc="best")
        ticklabel_format(style="sci", axis="y", scilimits=(0, 0))
        xlabel(L"Dimension of $\mathbf{x}$")
        ylabel(L"$(1-1/e)$-Regret")
        grid()
        play_sound()
        my_show()
    end
    generate_figure_regret()
    # generate_figure_regret_dim(1)
end
function experiment_design(dim_lam, dim_x_vectors, T, K)
    function generate_a_param(dim_lam, dim_x_vectors)
        mat = randn(dim_x_vectors, dim_lam)
        return mat
    end
    X_list = [generate_a_param(dim_lam, dim_x_vectors) for i in 1:T]
    n_constraints = 2
    A = rand(n_constraints, dim_lam)
    b = 1.0
    linear_projection = generate_projection_function(A, b)
    offset_constant = 1
    multp = 1
    function f(lam, X)
        lam += offset_constant
        lam *= multp
        M = 0
        for i in 1:dim_lam
            M += lam[i] * X[:, i] * X[:, i]'
        end
        return log(det(M))
    end
    function gradient(lam, X)
        lam += offset_constant
        lam *= multp
        M = 0
        for i in 1:dim_lam
            M += lam[i] * X[:, i] * X[:, i]'
        end
        grad = zeros(lam)
        for i in 1:size(lam)[1]
            grad[i] = X[:, i]' * inv(M) * X[:, i]
        end
        return grad
    end
    function gradient_orth(lam, X)
        return 1 ./ (lam + offset_constant)
    end
    function f_a(lam, X)
        lam += offset_constant
        M = 0
        for i in 1:dim_lam
            M += lam[i] * X[:, i] * X[:, i]'
        end
        return -trace(inv(M))
    end
    function gradient_a(lam, X)
        lam += offset_constant
        M = 0
        for i in 1:dim_lam
            M += lam[i] * X[:, i] * X[:, i]'
        end
        inv_M = inv(M)
        grad = zeros(lam)
        for i in 1:size(lam)[1]
            grad[i] = norm(inv_M * X[:, i])^2
        end
        return grad
    end
    function generate_figure_regret(f, gradient)
        reward_list = []
        reward_offline = offline(X_list, f, gradient, A, b, K)
        println(reward_offline)
        push!(reward_list, reward_offline - meta_frank_wolfe(dim_lam, X_list,
                                                     linear_projection, f, gradient, 0.8, T, K))
        push!(reward_list, reward_offline - meta_frank_wolfe(dim_lam, X_list,
                                                     linear_projection, f, gradient, 1, T, K))
        push!(reward_list, reward_offline - online_gradient_ascent(dim_lam, X_list,
                                                  linear_projection, f, gradient, 0.8, T))
        push!(reward_list, reward_offline - online_gradient_ascent(dim_lam, X_list,
                                                  linear_projection, f, gradient, 1, T))
        labels = [L"FW ($\eta=200$)", L"FW ($\eta=100$)", L"GA ($\eta=200$)", L"GA ($\eta=100$)", "OfflApprMax"]
        marker = cycle((".", ",", "+", "_", "o", "x", "*"))
        linecycler = cycle(("-", "--", "-.", ":"))
        print(reward_list)
        for zipped in zip(reward_list, labels, marker, linecycler)
            reward, label, marker_iter, line_iter = zipped
            plot(1:T, reward, label=label,  linestyle=line_iter, linewidth=global_width)
        end
        ticklabel_format(style="sci", axis="y", scilimits=(0, 0))
        legend(loc="best")
        xlabel("Iteration index")
        ylabel(L"$(1-1/e)$-Regret")
        grid()
        play_sound()
        my_show()
    end
    function generate_figure_regret_K()
        param_list = 1:20
        function f_cum(x)
            return sum([f(x, params) for params in X_list])
        end
        function gradient_cum(x)
            return sum([gradient(x, params) for params in X_list])
        end
        v = frank_wolfe2(dim_lam, A, b, gradient_cum, K)
        reward_offline = f_cum(v)
        ga1 = online_gradient_ascent(dim_lam, X_list, linear_projection, f, gradient, 0.8, T)[end]
        ga2 = online_gradient_ascent(dim_lam, X_list, linear_projection, f, gradient, 1, T)[end]
        reward_list = zeros(4, size(param_list)[1])
        for (idx, K) in enumerate(param_list)
            idx_th_col = reward_offline - [meta_frank_wolfe(dim_lam, X_list, linear_projection, f, gradient, 0.8, T, K)[end],
                meta_frank_wolfe(dim_lam, X_list, linear_projection, f, gradient, 1, T, K)[end], ga1, ga2]
            reward_list[:, idx] = idx_th_col
        end
        reward_list = [reward_list[i, :] for i in 1:size(reward_list)[1]]
        labels = [L"FW ($\eta=200$)", L"FW ($\eta=100$)", L"GA ($\eta=200$)", L"GA ($\eta=100$)"]
        marker = cycle((".", ",", "+", "_", "o", "x", "*"))
        linecycler = cycle(("-", "--", "-.", ":"))
        print(reward_list)
        for zipped in zip(reward_list, labels, marker, linecycler)
            reward, label, marker_iter, line_iter = zipped
            plot(param_list, reward, label=label, linestyle=line_iter, linewidth=global_width)
        end
        legend(loc="best")
        ticklabel_format(style="sci", axis="y", scilimits=(0, 0))
        xlabel(L"$K$")
        ylabel(L"$(1-1/e)$-Regret")
        grid()
        # savefig("design_regret_K3.pdf")
        play_sound()
        my_show()
    end
    generate_figure_regret(f, gradient)
end
experiment_nqp(10, 100, 100)
# experiment_design(20, 20, 50, 50)
# experiment_coverage(20, 60, 100, 100)


text:
using GeometricIntegrators.Config
using GeometricIntegrators.Integrators
using GeometricIntegrators.Solvers
using GeometricIntegrators.Tableaus
using GeometricIntegrators.TestProblems.LotkaVolterra
using GeometricIntegrators.Utils
using Test

using GeometricIntegrators.TestProblems.LotkaVolterra: Δt, nt

set_config(:nls_atol, 8eps())
set_config(:nls_rtol, 2eps())

idae = lotka_volterra_2d_idae()
pdae = lotka_volterra_2d_pdae()
vdae = lotka_volterra_2d_vdae()

int  = IntegratorFIRK(lotka_volterra_2d_ode(), getTableauGLRK(8), Δt)
sol  = integrate(int, nt)
refx = sol.q[:,end]

@testset "$(rpad("VPARK integrators",80))" begin

    dint = Integrator(idae, getTableauSymplecticProjection(:pglrk1ps, getCoefficientsGLRK(1), getCoefficientsGLRK(1)), Δt)
    dsol = integrate(dint, nt)

    # println(rel_err(dsol.q, refx))
    @test rel_err(dsol.q, refx) < 1E-6

    dint = Integrator(idae, getTableauSymplecticProjection(:pglrk2ps, getCoefficientsGLRK(2), getCoefficientsGLRK(2)), Δt)
    dsol = integrate(dint, nt)

    # println(rel_err(dsol.q, refx))
    @test rel_err(dsol.q, refx) < 1E-11

    dint = Integrator(idae, getTableauGLRKpSymplectic(1), Δt)
    dsol = integrate(dint, nt)

    # println(rel_err(dsol.q, refx))
    @test rel_err(dsol.q, refx) < 1E-6

    dint = Integrator(idae, getTableauGLRKpSymplectic(2), Δt)
    dsol = integrate(dint, nt)

    # println(rel_err(dsol.q, refx))
    @test rel_err(dsol.q, refx) < 1E-11

    dint = Integrator(idae, getTableauLobIIIAIIIB2pSymplectic(), Δt)
    dsol = integrate(dint, nt)

    # println(rel_err(dsol.q, refx))
    @test rel_err(dsol.q, refx) < 2E-6

    dint = Integrator(idae, getTableauLobIIIAIIIB3pSymplectic(), Δt)
    dsol = integrate(dint, nt)

    # println(rel_err(dsol.q, refx))
    @test rel_err(dsol.q, refx) < 8E-5

    dint = Integrator(idae, getTableauLobIIIAIIIB4pSymplectic(), Δt)
    dsol = integrate(dint, nt)

    # println(rel_err(dsol.q, refx))
    @test rel_err(dsol.q, refx) < 8E-12

end

@testset "$(rpad("VSPARK integrators",80))" begin
    ### VSPARK Integrators ###

    dint = IntegratorVSPARK(idae, getTableauSPARKGLRK(1), Δt)
    # dsol = integrate(dint, nt)
    # TODO
    # println(rel_err(dsol.q, refx))
    # @test rel_err(dsol.q, refx) < 1E-6

    ### VSPARKprimary Integrators ###

    dint = Integrator(idae, getTableauVSPARKGLRKpMidpoint(1), Δt)
    dsol = integrate(dint, nt)

    # println(rel_err(dsol.q, refx))
    @test rel_err(dsol.q, refx) < 1E-6

    dint = Integrator(idae, getTableauVSPARKGLRKpMidpoint(2), Δt)
    dsol = integrate(dint, nt)

    # println(rel_err(dsol.q, refx))
    @test rel_err(dsol.q, refx) < 1E-11

    dint = Integrator(idae, getTableauVSPARKGLRKpSymplectic(1), Δt)
    dsol = integrate(dint, nt)

    # println(rel_err(dsol.q, refx))
    @test rel_err(dsol.q, refx) < 1E-6

    dint = Integrator(idae, getTableauVSPARKGLRKpSymplectic(2), Δt)
    dsol = integrate(dint, nt)

    # println(rel_err(dsol.q, refx))
    @test rel_err(dsol.q, refx) < 1E-11

    dint = Integrator(idae, getTableauVSPARKGLRKpSymmetric(1), Δt)
    dsol = integrate(dint, nt)

    # println(rel_err(dsol.q, refx))
    @test rel_err(dsol.q, refx) < 1E-6

    dint = Integrator(idae, getTableauVSPARKGLRKpSymmetric(2), Δt)
    dsol = integrate(dint, nt)

    # println(rel_err(dsol.q, refx))
    @test rel_err(dsol.q, refx) < 1E-11

    dint = IntegratorVSPARKprimary(idae, getTableauVSPARKLobIIIAIIIB2pSymmetric(), Δt)
    dsol = integrate(dint, nt)

    # println(rel_err(dsol.q, refx))
    @test rel_err(dsol.q, refx) < 2E-6

    dint = IntegratorVSPARKprimary(idae, getTableauVSPARKLobIIIAIIIB3pSymmetric(), Δt)
    dsol = integrate(dint, nt)

    # println(rel_err(dsol.q, refx))
    @test rel_err(dsol.q, refx) < 8E-5

    dint = IntegratorVSPARKprimary(idae, getTableauVSPARKLobIIIAIIIB4pSymmetric(), Δt)
    dsol = integrate(dint, nt)

    # println(rel_err(dsol.q, refx))
    @test rel_err(dsol.q, refx) < 8E-12

    ### VSPARKsecondary Integrators ###

    dint = Integrator(vdae, getTableauVSPARKLobIIIAB(2), Δt)
    dsol = integrate(dint, nt)

    # println(rel_err(dsol.q, refx))
    @test rel_err(dsol.q, refx) < 4E-6

    dint = Integrator(vdae, getTableauVSPARKLobIIIAB(3), Δt)
    dsol = integrate(dint, nt)

    # println(rel_err(dsol.q, refx))
    @test rel_err(dsol.q, refx) < 2E-11

    dint = Integrator(vdae, getTableauVSPARKLobIIIAB(4), Δt)
    dsol = integrate(dint, nt)

    # println(rel_err(dsol.q, refx))
    @test rel_err(dsol.q, refx) < 1E-15

    dint = Integrator(vdae, getTableauVSPARKLobIIIC(2), Δt)
    dsol = integrate(dint, nt)

    # println(rel_err(dsol.q, refx))
    @test rel_err(dsol.q, refx) < 4E-6

    dint = Integrator(vdae, getTableauVSPARKLobIIIC(3), Δt)
    dsol = integrate(dint, nt)

    # println(rel_err(dsol.q, refx))
    @test rel_err(dsol.q, refx) < 2E-11

    dint = Integrator(vdae, getTableauVSPARKLobIIIC(4), Δt)
    dsol = integrate(dint, nt)

    # println(rel_err(dsol.q, refx))
    @test rel_err(dsol.q, refx) < 1E-15

    dint = Integrator(vdae, getTableauVSPARKLobIIID(2), Δt)
    dsol = integrate(dint, nt)

    # println(rel_err(dsol.q, refx))
    @test rel_err(dsol.q, refx) < 4E-6

    dint = Integrator(vdae, getTableauVSPARKLobIIID(3), Δt)
    dsol = integrate(dint, nt)

    # println(rel_err(dsol.q, refx))
    @test rel_err(dsol.q, refx) < 2E-11

    dint = Integrator(vdae, getTableauVSPARKLobIIID(4), Δt)
    dsol = integrate(dint, nt)

    # println(rel_err(dsol.q, refx))
    @test rel_err(dsol.q, refx) < 2E-15

    dint = Integrator(vdae, getTableauVSPARKLobIIIE(2), Δt)
    dsol = integrate(dint, nt)

    # println(rel_err(dsol.q, refx))
    @test rel_err(dsol.q, refx) < 1E-6

    dint = Integrator(vdae, getTableauVSPARKLobIIIE(3), Δt)
    dsol = integrate(dint, nt)

    # println(rel_err(dsol.q, refx))
    @test rel_err(dsol.q, refx) < 1E-11

    dint = Integrator(vdae, getTableauVSPARKLobIIIE(4), Δt)
    dsol = integrate(dint, nt)

    # println(rel_err(dsol.q, refx))
    @test rel_err(dsol.q, refx) < 1E-15

    dint = Integrator(vdae, getTableauVSPARKGLRKLobIIIAB(1), Δt)
    dsol = integrate(dint, nt)

    # println(rel_err(dsol.q, refx))
    @test rel_err(dsol.q, refx) < 4E-6

    dint = Integrator(vdae, getTableauVSPARKGLRKLobIIIAB(2), Δt)
    dsol = integrate(dint, nt)

    # println(rel_err(dsol.q, refx))
    @test rel_err(dsol.q, refx) < 1E-11

    dint = Integrator(vdae, getTableauVSPARKGLRKLobIIIAB(3), Δt)
    dsol = integrate(dint, nt)

    # println(rel_err(dsol.q, refx))
    @test rel_err(dsol.q, refx) < 1E-15

    dint = Integrator(vdae, getTableauVSPARKGLRKLobIIIC(1), Δt)
    dsol = integrate(dint, nt)

    # println(rel_err(dsol.q, refx))
    @test rel_err(dsol.q, refx) < 4E-6

    dint = Integrator(vdae, getTableauVSPARKGLRKLobIIIC(2), Δt)
    dsol = integrate(dint, nt)

    # println(rel_err(dsol.q, refx))
    @test rel_err(dsol.q, refx) < 1E-11

    dint = Integrator(vdae, getTableauVSPARKGLRKLobIIIC(3), Δt)
    dsol = integrate(dint, nt)

    # println(rel_err(dsol.q, refx))
    @test rel_err(dsol.q, refx) < 1E-15

    dint = Integrator(vdae, getTableauVSPARKGLRKLobIIID(1), Δt)
    dsol = integrate(dint, nt)

    # println(rel_err(dsol.q, refx))
    @test rel_err(dsol.q, refx) < 4E-6

    dint = Integrator(vdae, getTableauVSPARKGLRKLobIIID(2), Δt)
    dsol = integrate(dint, nt)

    # println(rel_err(dsol.q, refx))
    @test rel_err(dsol.q, refx) < 1E-11

    dint = Integrator(vdae, getTableauVSPARKGLRKLobIIID(3), Δt)
    dsol = integrate(dint, nt)

    # println(rel_err(dsol.q, refx))
    @test rel_err(dsol.q, refx) < 1E-15

    dint = Integrator(vdae, getTableauVSPARKGLRKLobIIIE(1), Δt)
    dsol = integrate(dint, nt)

    # println(rel_err(dsol.q, refx))
    @test rel_err(dsol.q, refx) < 1E-6

    dint = Integrator(vdae, getTableauVSPARKGLRKLobIIIE(2), Δt)
    dsol = integrate(dint, nt)

    # println(rel_err(dsol.q, refx))
    @test rel_err(dsol.q, refx) < 1E-11

    dint = Integrator(vdae, getTableauVSPARKGLRKLobIIIE(3), Δt)
    dsol = integrate(dint, nt)

    # println(rel_err(dsol.q, refx))
    @test rel_err(dsol.q, refx) < 1E-15
end

@testset "$(rpad("HPARK integrators",80))" begin

    dint = Integrator(pdae, getTableauHPARK(:hpark_glrk1, getCoefficientsGLRK(1), getCoefficientsGLRK(1)), Δt)
    dsol = integrate(dint, nt)

    # println(rel_err(dsol.q, refx))
    @test rel_err(dsol.q, refx) < 2E-6

    dint = Integrator(pdae, getTableauHPARK(:hpark_glrk2, getCoefficientsGLRK(2), getCoefficientsGLRK(2)), Δt)
    dsol = integrate(dint, nt)

    # println(rel_err(dsol.q, refx))
    @test rel_err(dsol.q, refx) < 8E-7

    dint = Integrator(pdae, getTableauHPARKGLRK(1), Δt)
    dsol = integrate(dint, nt)

    # println(rel_err(dsol.q, refx))
    @test rel_err(dsol.q, refx) < 2E-6

    dint = Integrator(pdae, getTableauHPARKGLRK(2), Δt)
    dsol = integrate(dint, nt)

    # println(rel_err(dsol.q, refx))
    @test rel_err(dsol.q, refx) < 8E-7

    dint = Integrator(pdae, getTableauHPARKLobIIIAIIIB2(), Δt)
    # dsol = integrate(dint, nt)
    # TODO
    # println(rel_err(dsol.q, refx))
    # @test rel_err(dsol.q, refx) < 2E-2

    dint = Integrator(pdae, getTableauHPARKLobIIIAIIIB3(), Δt)
    # dsol = integrate(dint, nt)
    # TODO
    # println(rel_err(dsol.q, refx))
    # @test rel_err(dsol.q, refx) < 8E-2

    dint = Integrator(pdae, getTableauHPARKLobIIIAIIIB4(), Δt)
    # dsol = integrate(dint, nt)
    # TODO
    # println(rel_err(dsol.q, refx))
    # @test rel_err(dsol.q, refx) < 8E-2
end

@testset "$(rpad("HSPARK integrators",80))" begin
    ### HSPARK Integrators ###

    dint = IntegratorHSPARK(pdae, getTableauSPARKGLRK(1), Δt)
    dsol = integrate(dint, nt)

    # println(rel_err(dsol.q, refx))
    @test rel_err(dsol.q, refx) < 1E-6

    dint = IntegratorHSPARK(pdae, getTableauSPARKGLRK(2), Δt)
    dsol = integrate(dint, nt)

    # println(rel_err(dsol.q, refx))
    @test rel_err(dsol.q, refx) < 1E-11

    ### HSPARKprimary Integrators ###

    dint = Integrator(pdae, getTableauHSPARKGLRKpSymmetric(1), Δt)
    dsol = integrate(dint, nt)

    # println(rel_err(dsol.q, refx))
    @test rel_err(dsol.q, refx) < 4E-6

    dint = Integrator(pdae, getTableauHSPARKGLRKpSymmetric(2), Δt)
    dsol = integrate(dint, nt)

    # println(rel_err(dsol.q, refx))
    @test rel_err(dsol.q, refx) < 4E-6

    dint = Integrator(pdae, getTableauHSPARKLobIIIAIIIB2pSymmetric(), Δt)
    dsol = integrate(dint, nt)

    # println(rel_err(dsol.q, refx))
    @test rel_err(dsol.q, refx) < 4E-6

    dint = Integrator(pdae, getTableauHSPARKLobIIIAIIIB3pSymmetric(), Δt)
    dsol = integrate(dint, nt)

    # println(rel_err(dsol.q, refx))
    @test rel_err(dsol.q, refx) < 4E-6

    dint = Integrator(pdae, getTableauHSPARKLobIIIAIIIB4pSymmetric(), Δt)
    dsol = integrate(dint, nt)

    # println(rel_err(dsol.q, refx))
    @test rel_err(dsol.q, refx) < 4E-6
end

ext: jl


text:
IntPoint(xy::NTuple{2,T}) where T<:Integer = IntPoint(xy[1],xy[2])

IntPoint2Tuple(a::IntPoint)=(a.X,a.Y)
IntPoint2Vector(i::IntPoint)=[i.X,i.Y]

function tupleextrema(x::Vector{NTuple{2,T}}) where T
    ex1 = extrema([xi[1] for xi in x])
    ex2 = extrema([xi[2] for xi in x])
    ( (ex1[1],ex2[1]), (ex1[2],ex2[2]) )
end

maxCint_one_two(T)=(typemax(Cint),T(1),T(2))
function _scale_j2c(x::Vector{Z},minx::Z,width::Z) where Z<:NTuple{2,T} where T
	n=length(x)
	ix = hcat( [[z...] for z in x]... ) # 2×length(x) Array
	maxi,o,t=maxCint_one_two(T)
	sx = maxi*(2*(ix.-minx)./width .- 1) # broadcast works for (Array,NTuple)
	cx = Array{Cint}(undef,size(sx))
	tb=sx.>maxi; ts=sx.<-maxi;
	cx[tb].=maxi; cx[ts].=-maxi
	representable = .~(tb.|ts)
	cx[representable] = round.(Cint,sx[representable])
	NTuple{2,Cint}[Tuple(cx[:,i]) for i=1:n]
end
function _scale_c2j(x::Vector{Z},minx::Y,width::Y) where {Z<:NTuple{2,I},Y<:NTuple{2,T}} where {I,T<:AbstractFloat}
	maxi,o,t=maxCint_one_two(T)
	[ convert.(T, (xi./maxi./t .+ (o/t) ).*width .+ minx) for xi in x ]
end
function _scale_c2j(x::Vector{Z},minx::Y,width::Y) where {Z<:NTuple{2,I},Y<:NTuple{2,T}} where {I,T<:Integer}
	maxi,o,t=maxCint_one_two(T)
	[   round.(T, (xi./maxi./t .+ (o/t) ).*width .+ minx) for xi in x ]
end

function min_width(x::Vector{NTuple{2,T}}) where T
    minx,maxx=tupleextrema(x)
    return (minx,(maxx.-minx))
end

to_intpoints(x::Vector{NTuple{2,T}}) where T = IntPoint.(_scale_j2c(x,min_width(x)...))
to_intpoints(x::Vector{Z},minx::Z,width::Z) where Z<:NTuple{2,T} where T = IntPoint.(_scale_j2c(x,minx,width))
function to_intpoints(x::Vector{Vector{NTuple{2,T}}}) where T
	minx,width=min_width(vcat(x...)) # determine width and minimum from *all* points
	[to_intpoints(xi,minx,width) for xi in x]
end

function from_intpoints(x::Vector{NTuple{2,T}},i::Vector{Vector{IntPoint}}) where T
    minx,width=min_width(x)
    [_scale_c2j(k,minx,width) for k in [IntPoint2Tuple.(j) for j in i]]
end
function from_intpoints(x::Vector{NTuple{2,T}},i::Vector{IntPoint}) where T
    minx,width=min_width(x)
    [_scale_c2j(k,minx,width) for k in IntPoint2Tuple.(i)]
end

function vector2ntuple(N::Integer,T::DataType,v::Vector{Vector{R}}) where R<:Number
    NTuple{N,T}[convert.(T,Tuple(x)) for x in v]
end
function vector2ntuple(N::Integer,T::DataType,v::Vector{Vector{Vector{R}}}) where R<:Number
    Vector{NTuple{N,T}}[vector2ntuple(N,T,x) for x in v]
end
function ntuple2vector(t::Vector{NTuple{N,T}}) where {N,T}
    Vector{T}[[x...] for x in t]
end
function ntuple2vector(t::Vector{Vector{NTuple{N,T}}}) where {N,T}
    Vector{Vector{T}}[ntuple2vector(x) for x in t]
end


text:
"""
    solve(s::Solver, m::DRO; kwargs...)

Solves given DRO model `m` using solver `s`.
"""
function solve(s::S, m::DRO; kwargs...) where {S<:Solver}
    reset_stats(s, m)

    @assert S <: Philpott && typeof(m.d) <: Ltwo "Philpott solver is supported only for the DRO model with l2 norm"

    p, t, bytes, gctime, memallocs = @timed optimal(s, m.d, m; kwargs...)
    add_stats(t, bytes)

    return p
end

"""
    solve(s::Solver, m::Simplex; kwargs...)

Solves given Simplex model `m` using solver `s`.
"""
function solve(s::S, m::Simplex; kwargs...) where {S<:Solver}
    reset_stats(s, m)

    @assert !(S <: Philpott) "Philpott solver is supported only for the DRO model with l2 norm"

    p, t, bytes, gctime, memallocs = @timed optimal(s, m; kwargs...)
    add_stats(t, bytes)

    return p
end

max_stars_repo_name: VaclavMacha/Projections

max_issues_repo_name: VaclavMacha/Projections

max_forks_repo_name: VaclavMacha/Projections


text:
using DSGE
using Test
using DataFrames: DataFrame
using Dates

@testset "Miscellaneous data handling functions" begin
    # Previous and next quarter arithmetic
    q = Dates.Date(1913,12,13)
    DSGE.prev_quarter()
    DSGE.next_quarter()
    @test DSGE.prev_quarter(q) == Dates.Date(1913,09,30)
    @test DSGE.next_quarter(q) == Dates.Date(1914,03,31)

    start_date = Dates.Date(2000,01,01)
    end_date   = Dates.Date(2010,01,01)
    DSGE.get_quarter_ends(start_date,end_date)
    DSGE.subtract_quarters(end_date, start_date)

    df1 = DataFrame(a = [1])
    global df2 = DataFrame()
    df1, df2 = DSGE.reconcile_column_names(df1, df2)
    @test names(df1) == names(df2)
    @test size(df2, 1) == 0

    quartertodate("11q3")
    quartertodate("1997q4")
    quartertodate("1985-Q1")
    @test_throws Meta.ParseError quartertodate("2005q9")
    @test_throws Meta.ParseError quartertodate("12345")

    global df = DataFrame(date = ["1913-12-23", "1992-11-14", "2002-01-01", "2014-12-19"],
                          x = [1, 2, missing, 4])
    DSGE.format_dates!(:date, df)
end

# Test hpfilter, ensuring missing data are handled correctly
missingfront = [missing; missing]
missingback  = [missing]
y  = [0.; 1.; 0.; 1.]
yn = [missingfront; y; missingback]
yf, yt   = hpfilter(y, 1600)
yfn, ytn = hpfilter(yn, 1600)

@testset "HP Filter and missing data handling" begin
    @test isequal(missingfront, yfn[1:length(missingfront)])
    @test isequal(missingfront, ytn[1:length(missingfront)])
    @test isequal(missingback, yfn[(end-length(missingback)+1):end])
    @test isequal(missingback, ytn[(end-length(missingback)+1):end])
end

nothing


text:
"""
Abstract supertype that extends `ManualNeighborhoodRule` for neighborhood-based
dispersal rules that update surounding cells based on the values of the
current cell, as if dispersing outwards from the current cell.

The result should be identical to [`InwardsDispersal`](@ref) but may be more
efficient than when a small number of cells are occupied. It is less efficient
when a large proportion of the grid is occupied.
"""
abstract type OutwardsDispersal{R,W} <: ManualNeighborhoodRule{R,W} end

@inline applyrule!(data, rule::OutwardsDispersal{R,W}, state, index) where {R,W} = begin
    state == zero(state) && return
    hood = neighborhood(rule)
    sum = mapsetneighbor!(data[W], hood, rule, state, index)
    update_state!(data[W], hood, state, index, sum)
    return
end

"""
    OutwardsBinaryDispersal(neighborhood)
    OutwardsBinaryDispersal(; neighborhood=DispersalKernel{3}())
    OutwardsBinaryDispersal{R,W}(neighborhood)

Cells in the surrounding neighborhood have some propability of
invasion if the current cell is occupied.

Pass grid name `Symbol`s to `R` and `W` type parameters to use specific grids.

$(FIELDDOCTABLE)
"""
@Kernel @Probabilistic struct OutwardsBinaryDispersal{R,W} <: OutwardsDispersal{R,W} end

"""
    OutwardsPopulationDispersal(neighborhood)
    OutwardsPopulationDispersal(; neighborhood=DispersalKernel{3}())
    OutwardsPopulationDispersal{R,W}(neighborhood)

Dispersal reduces the current cell population, increasing the populations of the
cells in the surrounding neighborhood.

Pass grid name `Symbol`s to `R` and `W` type parameters to use specific grids.

$(FIELDDOCTABLE)
"""
@Kernel struct OutwardsPopulationDispersal{R,W} <: OutwardsDispersal{R,W} end

@inline update_state!(grid, hood, state, index, sum) =
    @inbounds return grid[index...] -= sum
@inline update_state!(grid, hood, state::Bool, index, sum) = state

@inline setneighbor!(data::WritableGridData, hood::Neighborhood, rule::OutwardsPopulationDispersal,
                     state::AbstractFloat, hood_index, dest_index) = begin
    @inbounds propagules = state * kernel(hood)[hood_index...]
    @inbounds data[dest_index...] += propagules
    propagules
end

@inline setneighbor!(data::WritableGridData, hood::Neighborhood, rule::OutwardsBinaryDispersal,
                     state::Bool, hood_index, dest_index) = begin
    @inbounds rand() * kernel(hood)[hood_index...] > rule.prob_threshold || return zero(state)
    @inbounds data[dest_index...] |= oneunit(state)
    oneunit(state)
end


text:
immutable SkiLineDynamics{G<:Function}
    p::PowseekerMDP{G}
    pos::Vec2
    dir::Vec2
end

function rk4step(d::SkiLineDynamics, vel::Float64)
    h = d.p.dt
    vec = SVector(vel, 0.0)
    k1 = f(d, vec)
    k2 = f(d, vec+h/2*k1)
    k3 = f(d, vec+h/2*k2)
    k4 = f(d, vec+h*k3)
    return vec + h/6*(k1+2*k2+2*k3+k4)
end

function f(d::SkiLineDynamics, vec::Vec2)
    vel = vec[1]
    dist = vec[2]
    h = d.p.dt
    mfs = d.p.max_flat_speed
    theta = atan(dot(d.p.gradient(d.pos+dist*d.dir), d.dir)) # slope angle
    acc = d.p.force*max(0.0, mfs-vel)/mfs - 9.8*sin(theta) - sign(vel)*9.8*(vel/d.p.terminal_vel)^2
    return SVector(acc, vel*cos(theta))
end


text:
# ```@setup binning
# # all examples are run on an example global.toml and downsized maps.
# ARGS = ["example.toml"]
# ```

configfile = ARGS[1]

## get modules and utility functions
using Plots
using TOML
using Healpix
include("util.jl")

config = TOML.parsefile(configfile)
nside = config["general"]["nside"]
lmax = nside2lmax(nside)

pl = PlanckReferenceCov(joinpath(config["scratch"], "plicref"))

## read binning scheme
binfile = joinpath(@__DIR__, "../", "input", "binused.dat")
P, lb = util_planck_binning(binfile; lmax=lmax);

##
freq1, freq2 = "100", "100"
Wl = util_planck_beam_Wl(freq1, "hm1", freq2, "hm2", :EE, :EE;
    lmax=lmax, beamdir=joinpath(config["scratch"], "beams"))
plot(Wl, yaxis=:log10, label="\$B_{\\ell}\$")

##
run_name = config["general"]["name"]
mapids = [k for k in keys(config["map"])]
spectrapath = joinpath(config["scratch"], "rawspectra")

mapid1 = "P$(freq1)hm1"
mapid2 = "P$(freq2)hm2"
spec = CSV.read(joinpath(spectrapath,"$(run_name)_$(mapid1)x$(mapid2).csv"), DataFrame)
plot(spec.ell, spec.ell.^2 .* spec.EE, label="$(run_name)_$(mapid1)x$(mapid2)",
    xlabel="multipole moment", ylabel="\$\\ell^2 C_{\\ell}^{EE}\$", xlim=(0,2nside))

##
lbref, cbref, err_ref, _ = extract_spec_and_cov(pl, "EE_$(freq1)x$(freq2)")

##
cl = SpectralVector(copy(spec.EE))
this_lmax = max(lastindex(cl), lastindex(Wl))
pixwinT, pixwinP = pixwin(nside; pol=true)
cl[0:this_lmax] ./= Wl[0:this_lmax] .* pixwinP[1:(this_lmax+1)].^2
cl[0:1] .= 0.0
cb = P * parent(cl)
planck_bin_choice = findfirst(lb .≥ lbref[1]):findlast(lb .≤ lbref[end])

nside_cut = 1:length(planck_bin_choice)
plot(lb[planck_bin_choice], (cb[planck_bin_choice] .- cbref[nside_cut]) ./ err_ref[nside_cut],
    ylim=(-2,2), label=label="$(run_name)_$(mapid1)x$(mapid2)")

max_forks_repo_name: sgiardie/PSpipe


text:
predicate_spec =
[
    (:name, AbstractString, nothing),
    (:_type, AbstractString, nothing),
    (:operands, Vector, nothing),
    (:range, Vector, nothing),
    (:scale, Vector, nothing),
    (:item, Any, nothing),
]

primitivefactory(:VegaPredicate, predicate_spec)

max_stars_repo_name: jpthor/vega


text:
using Test
using IntegerIB

#bach chorale data, taken from Nori Jacoby's webpage.
bach_histogram = [0.00001 392 14 582 896 262 106; 12 0.0001 15 5 674 17 182; 13 6 0.0001 44 12 32 3; 239 134 16 0.0001 380 13 173; 2099 25 46 93 0.00001 211 2; 96 206 25	95	172	0.00001 36; 499 1 4 9 24 20 0.00001]
pxy = bach_histogram./sum(bach_histogram)
model = IB(pxy, 0.95, "DIB")
brute_optimize!(model)
h, ixt, iyt, L = calc_metrics(model)
@test size(model.qt_x, 1) == 5
@test round(h, digits = 2) == 2.06
@test round(ixt, digits = 2) == 2.06
@test round(iyt, digits = 2) == 0.78
@test round(L, digits = 2) == 1.32
#toy time-series presenting one predictable pattern "b" -> "a"
toy_ts = ["a", "b", "a", "b", "a", "c", "d", "b", "a", "a", "d", "b", "a", "b", "a", "d", "c", "d", "b", "a", "d", "a", "c", "b", "a", "a", "b", "a", "c", "b", "a"]
model = IB(toy_ts, 500)
IB_optimize!(model)
@test (size(model.qt_x, 1) == 2 || size(model.qt_x, 1) == 3) #we expect the algorithm to cluster at least 2 labels together.


text:
using ArrayLayouts, Test

include("test_layouts.jl")
include("test_muladd.jl")
include("test_ldiv.jl")


text:
#Computes the d-th power residue symbol (a/b)_d in F_q
#with d|q-1

function power_residue(a::PolyElem{T}, b::PolyElem{T}, d::fmpz, q::fmpz) where T<: RingElem
    if mod(q-1,d)!=0
        return error("Must have d|q-1")
    end
    e = div(q-1, d)

    if b==0 || a==0
        return error("Must have a,b!=0")
    end

    if gcd(a,b)!=1
        return 0
    end

    a = mod(a,b)

    #top part constant -> formula
    if degree(a)==0
        return powmod( powmod(fmpz(coeff(a,0)), e, q),  degree(b), q)
    end

    #switch top and bottom part -> recursion (general reciprocity law)
    symbol = mod(  mod( powmod(powmod(fmpz(leading_coefficient(a)), e, q), degree(b), q) *
    powmod(powmod(fmpz(leading_coefficient(b)), e, q), -degree(a), q),  q) *
    power_residue(b,a,d,q),   q)

    if isodd(e) && isodd(degree(a)) && isodd(degree(b))
        symbol = mod( (-1) * symbol, q)
    end

    return symbol
end


text:
# Each character on a computer is assigned a unique code and the preferred
# standard is ASCII (American Standard Code for Information Interchange). For
# example, uppercase A = 65, asterisk (*) = 42, and lowercase k = 107.
#
# A modern encryption method is to take a text file, convert the bytes to
# ASCII, then XOR each byte with a given value, taken from a secret key. The
# advantage with the XOR function is that using the same encryption key on the
# cipher text, restores the plain text; for example, 65 XOR 42 = 107, then 107
# XOR 42 = 65.
#
# For unbreakable encryption, the key is the same length as the plain text
# message, and the key is made up of random bytes. The user would keep the
# encrypted message and the encryption key in different locations, and
# without both "halves", it is impossible to decrypt the message.
#
# Unfortunately, this method is impractical for most users, so the modified
# method is to use a password as a key. If the password is shorter than the
# message, which is likely, the key is repeated cyclically throughout the
# message. The balance for this method is using a sufficiently long password
# key for security, but short enough to be memorable.
#
# Your task has been made easy, as the encryption key consists of three lower
# case characters. Using p059_cipher.txt (right click and 'Save Link/Target
# As...'), a file containing the encrypted ASCII codes, and the knowledge that
# the plain text must contain common English words, decrypt the message and
# find the sum of the ASCII values in the original text.

using ProjectEulerSolutions

# Helper function for finding the most common element in an array.
function most_common(input::Array{Int8})::Int8
    d = Dict{eltype(input), Int}()
    for val in input
        if isa(val, Number) && isnan(val)
            continue
        end
        d[val] = get!(d, val, 0) + 1
    end
    top_key = 0
    most_val = 0
    for k in keys(d)
        if d[k] > most_val
            most_val = d[k]
            top_key = k
        end
    end
    return top_key
end

# Brute force solution, try all keys and keep the one with the most number
# of spaces.
function p059solution_bruteforce(path::String="")::Integer
    if !isfile(path)
        if path != ""
            println("Cannot find path: '", path, "'")
        end
        return -1
    end

    keylen = 3
    all_chars = [parse.(Int8, split(line, ',')) for line in eachline(path)][1]
    orig_len = length(all_chars)

    # Make a multiple of keylen for broadcasting our key
    diff_len = length(all_chars) % keylen
    for i in 1:diff_len
        push!(all_chars, 0)
    end
    all_chars = reshape(all_chars, keylen, :)

    start_keyval = Int8('a')
    end_keyval = Int8('z')
    key = fill(Int8(start_keyval), keylen)

    max_space = 0
    best_key = copy(key)
    for i in start_keyval:end_keyval
        key[1] = i
        for j in start_keyval:end_keyval
            key[2] = j
            for k in start_keyval:end_keyval
                key[3] = k
                x = all_chars .⊻ key
                count_space = sum(x .== Int8(' '))
                if count_space > max_space
                    max_space = count_space
                    best_key = copy(key)
                end
            end
        end
    end
    return sum(reshape(all_chars .⊻ best_key, length(all_chars), :)[1:orig_len])
end

# Here we assume that space is the most common character and xor the most common
# value from the file with ' ' for each section of the key to find the final key
# value.
function p059solution_xor(path::String="")::Integer
    if !isfile(path)
        if path != ""
            println("Cannot find path: '", path, "'")
        end
        return -1
    end

    keylen = 3
    all_chars = [parse.(Int8, split(line, ',')) for line in eachline(path)][1]
    orig_len = length(all_chars)
    # Make a multiple of keylen for broadcasting our key
    diff_len = length(all_chars) % keylen
    for i in 1:diff_len
        push!(all_chars, 0)
    end

    # Most likely char is a space.  XOR to find the key value that gives a space
    # for each key
    all_chars = reshape(all_chars, keylen, :)
    key = fill(Int8('a'), keylen)
    for i in 1:keylen
        key[i] = most_common(all_chars[i,:]) ⊻ Int8(' ')
    end

    return sum(reshape(all_chars .⊻ key, length(all_chars), :)[1:orig_len])
end

p059 = Problems.Problem(Dict("Brute Force" => p059solution_bruteforce,
                             "XOR" => p059solution_xor))

input_path = "data/p059_cipher.txt"

Problems.benchmark(p059, input_path)

max_stars_repo_name: gnujosh/julia-euler-project

max_issues_repo_name: gnujosh/julia-euler-project

max_forks_repo_name: gnujosh/julia-euler-project


text:
const NO_DIRECTION = 0
const NORTH = 1
const EAST = 2
const SOUTH = 3
const WEST = 4

const CARDINAL_POINTS = (NORTH, EAST, SOUTH, WEST)

const REAL = 1
const DEPARTURE = 2
const ARRIVAL = 3


text:
# target tracking -- cylindrical manifold (r, θ)

using Distributed

addprocs(3)

using Caesar
using DifferentialEquations
using ApproxManifoldProducts
using Caesar, RoME
using DataFrames
using TransformUtils

@everywhere using Caesar

@everywhere setForceEvalDirect!(true)

## Plotting libraries

using Colors
using Plots # for ODE
using KernelDensityEstimatePlotting, RoMEPlotting
using AbstractPlotting, Makie
using FileIO
using DelimitedFiles

## parameters for experiment

include(joinpath(dirname(@__FILE__),"GenerateSimData.jl"))

mc = 1; N = 100;

for mc = 1:50, N = [200;]

## Start building the factor graph in polar coordinates
include(joinpath(dirname(@__FILE__),"BuildFactorGraph.jl"))

##

# writeGraphPdf(fg, show=true)

# ensureAllInitialized!(fg)
getSolverParams(fg).drawtree=false
getSolverParams(fg).showtree=false
getSolverParams(fg).N=N
tree, smt, hist = solveTree!(fg)

# savejld(fg, file="fg100_rstd10_sol3.jld2")

##  Prep data for plotting

include(joinpath(dirname(@__FILE__),"PrepPlottingData.jl"))

## Calculate error statistics

mkdir( joinpath(dirname(@__FILE__),"exports","N$(N)_MC$(mc)") )
fname = joinpath(dirname(@__FILE__),"exports","N$(N)_MC$(mc)","meanRanErr.txt")
writedlm(fname, dferr["mean".==dferr[:,3],2], ',' )
fname = joinpath(dirname(@__FILE__),"exports","N$(N)_MC$(mc)","maxRanErr.txt")
writedlm(fname, dferr["max".==dferr[:,3],2], ',')
fname = joinpath(dirname(@__FILE__),"exports","N$(N)_MC$(mc)","meanAngErr.txt")
writedlm(fname, dferrang["mean".==dferrang[:,3],2], ',')
fname = joinpath(dirname(@__FILE__),"exports","N$(N)_MC$(mc)","maxAngErr.txt")
writedlm(fname, dferrang["max".==dferrang[:,3],2], ',')

end

## Generate all plots

include(joinpath(dirname(@__FILE__),"GeneratePlots.jl"))

## Look a factor graph local output -- More debugging

# pll, plc = plotLocalProductCylinder(fg, :t0, show=false, scale=0.1)

0

#


text:
module PackageWithBuildSpecificTestDeps

f(x) = x

end


text:
export PrioritizedSweepingSampleModel, update!, sample

using DataStructures

import StatsBase: sample

"""
    PrioritizedSweepingSampleModel(θ::Float64=1e-4)

See more details at Section (8.4) on Page 168 of the book *Sutton, Richard S., and Andrew G. Barto. Reinforcement learning: An introduction. MIT press, 2018.*
"""
mutable struct PrioritizedSweepingSampleModel <: AbstractSampleBasedModel
    experiences::Dict{Tuple{Any,Any},Tuple{Float64,Bool,Any}}
    PQueue::PriorityQueue{Tuple{Any,Any},Float64}
    predecessors::Dict{Any,Set{Tuple{Any,Any,Float64,Bool}}}
    θ::Float64
    sample_count::Int
    PrioritizedSweepingSampleModel(θ::Float64 = 1e-4) =
        new(
            Dict{Tuple{Any,Any},Tuple{Float64,Bool,Any}}(),
            PriorityQueue{Tuple{Any,Any},Float64}(Base.Order.Reverse),
            Dict{Any,Set{Tuple{Any,Any,Float64,Bool}}}(),
            θ,
            0,
        )
end

function extract_transitions(
    buffer::EpisodeTurnBuffer,
    model::PrioritizedSweepingSampleModel,
)
    if length(buffer) > 0
        state(buffer)[end-1],
        action(buffer)[end-1],
        reward(buffer)[end],
        terminal(buffer)[end],
        state(buffer)[end]
    else
        nothing
    end
end

function update!(m::PrioritizedSweepingSampleModel, transition, P)
    s, a, r, d, s′ = transition
    m.experiences[(s, a)] = (r, d, s′)
    if P >= m.θ
        m.PQueue[(s, a)] = P
    end
    if !haskey(m.predecessors, s′)
        m.predecessors[s′] = Set{Tuple{Any,Any,Float64,Bool}}()
    end
    push!(m.predecessors[s′], (s, a, r, d))
end

function sample(m::PrioritizedSweepingSampleModel)
    if length(m.PQueue) > 0
        s, a = dequeue!(m.PQueue)
        r, d, s′ = m.experiences[(s, a)]
        m.sample_count += 1
        s, a, r, d, s′
    else
        nothing
    end
end


text:
function numericRoot(residFnc::Function, measurement, parameters, x0::Vector{Float64})
  return (nlsolve(   (res, X) -> residFnc(res, measurement, parameters, X), x0 )).zero
end

function freshSamples!(ccwl::CommonConvWrapper, N::Int=1)
  ccwl.measurement = getSample(ccwl.usrfnc!, N)
  nothing
end

function shuffleXAltD(X::Vector{Float64}, Alt::Vector{Float64}, d::Int, p::Vector{Int})
  # n = length(X)
  Y = deepcopy(Alt)
  for i in 1:d
    Y[p[i]] = X[i]
  end
  return Y
end

# Shuffle incoming X into random positions in fr.Y
# shuffled fr.Y will be placed back into fr.X[:,fr.gwp.particleidx] upon fr.gwp.usrfnc(x, res)
function shuffleXAltD!(ccwl::CommonConvWrapper, X::Vector{Float64})
  # populate defaults from existing values
  for i in 1:ccwl.xDim
    ccwl.cpt[Threads.threadid()].Y[i] = ccwl.cpt[Threads.threadid()].X[i, ccwl.cpt[Threads.threadid()].particleidx]
  end
  # populate as many measurment dimensions randomly for calculation
  for i in 1:ccwl.zDim
    ccwl.cpt[Threads.threadid()].Y[ccwl.cpt[Threads.threadid()].p[i]] = X[i]
  end
  nothing
end

function (ccw::CommonConvWrapper)(res::Vector{Float64}, x::Vector{Float64})
  shuffleXAltD!(ccw, x)
  ccw.params[ccw.varidx][:, ccw.cpt[Threads.threadid()].particleidx] = ccw.cpt[Threads.threadid()].Y
  ret = ccw.usrfnc!(res, ccw.cpt[Threads.threadid()].factormetadata, ccw.cpt[Threads.threadid()].particleidx, ccw.measurement, ccw.params[ccw.cpt[Threads.threadid()].activehypo]...) # optmize the view here, re-use the same memory
  return ret
end

function (ccw::CommonConvWrapper)(x::Vector{Float64})
  ccw.params[ccw.varidx][:, ccw.cpt[Threads.threadid()].particleidx] = x #ccw.Y
  ccw.usrfnc!(ccw.cpt[Threads.threadid()].res, ccw.cpt[Threads.threadid()].factormetadata, ccw.cpt[Threads.threadid()].particleidx, ccw.measurement, ccw.params[ccw.cpt[Threads.threadid()].activehypo]...)
end

function numericRootGenericRandomizedFnc!(
            ccwl::CommonConvWrapper{T};
            perturb::Float64=1e-10,
            testshuffle::Bool=false ) where {T <: FunctorPairwiseMinimize}
  #
  fill!(ccwl.cpt[Threads.threadid()].res, 0.0) # 1:frl.xDim
  r = optimize( ccwl, ccwl.cpt[Threads.threadid()].X[:, ccwl.cpt[Threads.threadid()].particleidx] ) # ccw.gg
  # TODO -- clearly lots of optmization to be done here
  ccwl.cpt[Threads.threadid()].Y[:] = r.minimizer
  ccwl.cpt[Threads.threadid()].X[:,ccwl.cpt[Threads.threadid()].particleidx] = ccwl.cpt[Threads.threadid()].Y
  nothing
end

## TODO desperately needs cleaning up and refactoring
# Solve free variable x by root finding residual function fgr.usrfnc(x, res)
# randomly shuffle x dimensions if underconstrained by measurement z dimensions
# small random perturbation used to prevent trivial solver cases, div by 0 etc.
# result stored in fgr.Y
# fr.X must be set to memory ref the param[varidx] being solved, at creation of fr
function numericRootGenericRandomizedFnc!(
            ccwl::CommonConvWrapper{T};
            perturb::Float64=1e-10,
            testshuffle::Bool=false ) where {T <: FunctorPairwise}
  #
  # ststr = "thrid=$(Threads.threadid()), zDim=$(ccwl.zDim), xDim=$(ccwl.xDim)\n"
  # ccall(:jl_, Void, (Any,), ststr)
  if ccwl.zDim < ccwl.xDim && !ccwl.partial || testshuffle
    # less measurement dimensions than variable dimensions -- i.e. shuffle
    shuffle!(ccwl.cpt[Threads.threadid()].p)
    for i in 1:ccwl.xDim
      ccwl.cpt[Threads.threadid()].perturb[1:ccwl.zDim] = perturb*randn(ccwl.zDim)
      ccwl.cpt[Threads.threadid()].X[ccwl.cpt[Threads.threadid()].p[1:ccwl.zDim], ccwl.cpt[Threads.threadid()].particleidx] += ccwl.cpt[Threads.threadid()].perturb
      r = nlsolve(  ccwl,
                    ccwl.cpt[Threads.threadid()].X[ccwl.cpt[Threads.threadid()].p[1:ccwl.zDim], ccwl.cpt[Threads.threadid()].particleidx] # this is x0
                 )
      if r.f_converged
        shuffleXAltD!( ccwl, r.zero )
        break;
      else
        # TODO -- report on this bottleneck, useful for optimization of code
        # @show i, ccwl.p, ccwl.xDim, ccwl.zDim
        temp = ccwl.cpt[Threads.threadid()].p[end]
        ccwl.cpt[Threads.threadid()].p[2:end] = ccwl.cpt[Threads.threadid()].p[1:(end-1)]
        ccwl.cpt[Threads.threadid()].p[1] = temp
        if i == ccwl.xDim
          error("numericRootGenericRandomizedFnc could not converge, i=$(i), ccwl.usrfnc!=$(typeof(ccwl.usrfnc!))")
        end
      end
    end
    #shuffleXAltD!( ccwl, r.zero ) # moved up
  elseif ccwl.zDim >= ccwl.xDim && !ccwl.partial
    # equal or more measurement dimensions than variable dimensions -- i.e. don't shuffle
    ccwl.cpt[Threads.threadid()].perturb[1:ccwl.xDim] = perturb*randn(ccwl.xDim)
    ccwl.cpt[Threads.threadid()].X[1:ccwl.xDim, ccwl.cpt[Threads.threadid()].particleidx] += ccwl.cpt[Threads.threadid()].perturb[1:ccwl.xDim] # moved up

    # str = "nlsolve, thrid_=$(Threads.threadid()), partidx=$(ccwl.cpt[Threads.threadid()].particleidx), X=$(ccwl.cpt[Threads.threadid()].X[1:ccwl.xDim,ccwl.cpt[Threads.threadid()].particleidx])"
    # ccall(:jl_, Void, (Any,), str)
    r = nlsolve( ccwl, ccwl.cpt[Threads.threadid()].X[1:ccwl.xDim,ccwl.cpt[Threads.threadid()].particleidx] )
    # ccall(:jl_, Void, (Any,), "nlsolve.zero=$(r.zero)")
    if sum(isnan.(( r ).zero)) == 0
      ccwl.cpt[Threads.threadid()].Y[1:ccwl.xDim] = ( r ).zero
    else
      # TODO print this output as needed
      # warn("got NaN, ccwl.cpt[Threads.threadid()].particleidx = $(ccwl.cpt[Threads.threadid()].particleidx), r=$(r)")
      str = "ccw.thrid_=$(Threads.threadid()), got NaN, ccwl.cpt[Threads.threadid()].particleidx = $(ccwl.cpt[Threads.threadid()].particleidx), r=$(r)\n"
      ccall(:jl_, Void, (Any,), str)
      ccall(:jl_, Void, (Any,), ccwl.usrfnc!)
      for thatlen in 1:length(ccwl.params)
        str = "thatlen=$thatlen, ccwl.params[thatlen][:, ccwl.cpt[Threads.threadid()].particleidx]=$(ccwl.params[thatlen][:, ccwl.cpt[Threads.threadid()].particleidx])\n"
        ccall(:jl_, Void, (Any,), str)
      end
    end
  elseif ccwl.partial
    # improve memory management in this function
    ccwl.cpt[Threads.threadid()].p[1:length(ccwl.usrfnc!.partial)] = Int[ccwl.usrfnc!.partial...] # TODO -- move this line up and out of inner loop
    r = nlsolve(  ccwl,
                  ccwl.cpt[Threads.threadid()].X[ccwl.cpt[Threads.threadid()].p[1:ccwl.zDim], ccwl.cpt[Threads.threadid()].particleidx] # this is x0
               )
    shuffleXAltD!( ccwl, r.zero )
  else
    ccall(:jl_, Void, (Any,), "ERROR: Unresolved numeric solve case")
    error("Unresolved numeric solve case")
  end
  ccwl.cpt[Threads.threadid()].X[:,ccwl.cpt[Threads.threadid()].particleidx] = ccwl.cpt[Threads.threadid()].Y
  nothing
end

"""
    $(SIGNATURES)

Perform multimodal incremental smoothing and mapping (mm-iSAM) computations over given factor graph `fgl::FactorGraph` on the local computer.  A pdf of the Bayes (Junction) tree will be generated in the working folder with `drawpdf=true`
"""
function batchSolve!(fgl::FactorGraph; drawpdf::Bool=false)
  tree = wipeBuildNewTree!(fgl, drawpdf=drawpdf)
  inferOverTree!(fgl, tree)
  tree
end

#


text:
module TestModes
    using Base.Test
    import FiniteDiff

    @testset "Mode type definitions and constructors" begin
        modes = (
            FiniteDiff.ForwardMode,
            FiniteDiff.BackwardMode,
            FiniteDiff.CentralMode,
            FiniteDiff.ComplexMode,
            FiniteDiff.HessianMode,
        )

        for mode in modes
            m = mode()
            @test isa(m, mode)
            @test isa(m, FiniteDiff.Mode)
        end
    end
end


text:
@test approxeq( differential_entropy([0.0 1.0 2.0 3.0]), 2.77, 1e-2 )

@test approxeq( differential_entropy([0 1 2 3]), 2.77, 1e-2 )
@test approxeq( differential_entropy([0 1 3 2]), 2.77, 1e-2 )
@test approxeq( differential_entropy([0 2 1 3]), 2.77, 1e-2 )
@test approxeq( differential_entropy([0 2 3 1]), 2.77, 1e-2 )
for samples in permutations([0 1 2 3])
  @test approxeq( differential_entropy(samples'), 2.77, 1e-2 )
end

@test approxeq( differential_entropy([0.0 1.0 2.0 2.0]), 2.23, 1e-2 )

@test approxeq( differential_entropy([0 1 2 2]), 2.23, 1e-2 )
@test approxeq( differential_entropy([0 2 1 2]), 2.23, 1e-2 )
@test approxeq( differential_entropy([0 2 2 1]), 2.23, 1e-2 )
for samples in permutations([0 1 2 2])
  @test approxeq( differential_entropy(samples'), 2.23, 1e-2 )
end

@test approxeq( differential_entropy([0 1 1 2]), 1.98, 1e-2 )
for samples in permutations([0 1 1 2])
  @test approxeq( differential_entropy(samples'), 1.98, 1e-2 )
end

# Actual entropy is 0.25*log(0.25) + 0.75*log(0.75)
@test approxeq( differential_entropy([0 1 1 1]), 1.48, 1e-2 )
@test approxeq( differential_entropy([1 2 2 2]), 1.48, 1e-2 )
@test approxeq( differential_entropy([2 3 3 3]), 1.48, 1e-2 )
for samples in permutations([0 1 1 1])
  @test approxeq( differential_entropy(samples'), 1.48, 1e-2 )
end

@test approxeq( differential_entropy([0 1 0 1]), 1.48, 1e-2 )
@test approxeq( differential_entropy([1 0 0 1]), 1.48, 1e-2 )
for samples in permutations([0 0 1 1])
  @test approxeq( differential_entropy(samples'), 1.48, 1e-2 )
end

# Uniform distribution of width alpha should have entropy of log_2(alpha)
@repeatedly begin
  N = 1000
  alpha = 10.0 * rand()
  samples_from_uniform_distr = alpha * rand(1, N)
  expected = log(2, alpha)

  for k in 3:6
    ent = differential_entropy(samples_from_uniform_distr, k)
    #println("estimated = $ent, expected = $expected")
    @test approxeq( ent, expected, 0.15 )
  end
end

#
# Compare to analytically calculated entropies for known distributions, as listed on:
#  http://en.wikipedia.org/wiki/Differential_entropy
#
using Distributions
N = 1000

function diff_entropy_approx_eq(distr, expected, delta = 0.10, N = N)
  samples = rand(distr, N)
  ent = differential_entropy(samples')
  approxeq(ent, expected, delta)
end

# Convert from nats to bits
nats_to_bits(natsvalue) = natsvalue / log(2.0)

@repeatedly begin

  # Uniform distribution
  a = rand(-100.0:100.0)
  b = rand((a+0.01):1e-5:(a+100.0))
  expected_entropy = nats_to_bits( log(b-a) )
  @test diff_entropy_approx_eq( Uniform(a, b), expected_entropy, 0.20 )

  # Normal distribution
  mu = rand(-10.0:1e-6:10.0)
  sigma = rand(0.0:1e-2:100.0)
  expected_entropy = nats_to_bits( log(sigma * sqrt(2*pi*e)) )
  @test diff_entropy_approx_eq( Normal(mu, sigma), expected_entropy, 0.20 )

  # Exponential distribution
  lambda = rand(0.0:1e-6:10.0)
  scale = 1/lambda
  expected_entropy = nats_to_bits( 1 - log(lambda) )
  @test diff_entropy_approx_eq( Exponential(scale), expected_entropy, 0.20 )

  # Rayleigh distribution
  expected_entropy = nats_to_bits( 1 + log(scale / sqrt(2)) + eulergamma / 2 )
  @test diff_entropy_approx_eq( Rayleigh(scale), expected_entropy, 0.15 )

  # Laplace distribution
  expected_entropy = nats_to_bits( 1 + log(2 * scale) )
  @test diff_entropy_approx_eq( Laplace(mu, scale), expected_entropy, 0.20 )

  # Gamma distribution
  alpha = rand(1:20)
  scale = rand(0.1:1e-3:3.0)
  expected_entropy = nats_to_bits( log(scale * gamma(alpha)) + (1 - alpha) * digamma(alpha) + alpha )
  @test diff_entropy_approx_eq( Gamma(alpha, scale), expected_entropy, 0.20 )

end


text:
module TestSpecies

  using GBIF
  using Test

  iver = taxon("Iris versicolor", rank=:SPECIES)
  @test iver.species == Pair("Iris versicolor", 5298019)

  iver_id = taxon(5298019)  # test using taxon id
  @test iver_id.species == iver.species

  @test iver_id == taxon(Pair("Iris versicolor", 5298019))

  i_sp = taxon(iver.genus.first; rank=:GENUS, family=iver.family.first, strict=false)
  @test ismissing(i_sp.species)

  iver_occ = occurrences(iver)
  @test typeof(iver_occ) <: GBIFRecords

  iver_occ_spain = occurrences(iver, Dict{Any,Any}("country" => "ES"))
  @test typeof(iver_occ_spain) <: GBIFRecords

end


text:
module Configurations

import JLD2

using ..Trajectories
using ..Metrics
using ..Force
using ..UpdateMethods

export Configuration

struct Configuration
    traj :: Trajectory
    workheat :: WorkHeat
end

function Configuration(x0::Real, f::AbstractForce, h::Heun, τ::Real)
    t0 = 0.0
    traj = heun_fill(x0, h, (p, t) -> total_force(f, p, t), (t0, τ))
    workheat = get_workheat(traj, f, h.step)
    Configuration(traj, workheat)
end

function save_configuration(conf::Configuration, filename)
    @JLD2.save filename conf
end

function load_configuration(filename)::Configuration
    @JLD2.load filename conf
    return conf
end

end


text:
using NearestNeighbors
using Benchmarks
using Mmap

runtimes = []
runtimesreordered = []

function create_tree(n, reorder=false)
    filename = tempname()
    d = 10
    data = Mmap.mmap(filename, Matrix{Float32}, (d, n))
    data[:] = rand(Float32, d, n)
    if reorder
        reorderbuffer = Mmap.mmap(filename, Matrix{Float32}, (d, n))
        t = injectdata(DataFreeTree(KDTree, data, reorderbuffer = reorderbuffer), reorderbuffer)
    else
        t = injectdata(DataFreeTree(KDTree, data), data)
    end

    return t, data, filename
end

function knnbench(tree, data, n, N)
    ind = rand(1:n, N)
    knn(tree, data[:,ind], 3)[2]
end

function bench()
    runtimes = []
    runtimesreordered = []
    ns = [10, 100, 1_000, 10_000, 100_000, 1_000_000, 10_000_000]
    for n in ns
        t, data, filename = create_tree(n)
        tr, datar, filenamer = create_tree(n, true)

        bm = @benchmark knnbench(t, data, n, 1000)
        push!(runtimes, mean(bm.samples.elapsed_times) / 1e9)

        bmr = @benchmark knnbench(tr, datar, n, 1000)
        push!(runtimesreordered, mean(bmr.samples.elapsed_times) / 1e9)

        rm(filename)
        rm(filenamer)
    end

    println("Speedups through reordering:")
    for i in 1:length(ns)
        println("$(ns[i]): $(runtimes[i] ./ runtimesreordered[i])")
    end
    return
end
bench()


text:
@doc raw"""
    module NMR

provides routines for processing 1D NMR spectra, as well as for simple Hilbert space
spin dynamics simulations.

### Version History

#### v0.6.5
- added `medianBaseline()` for baseline correction
- added an entry for DSS in the HMDB table

"""
module NMR

import Pkg

   __precompile__(true);

   include("DataSet.jl")
   include("AutoPhase.jl")
   include("Scaling.jl")
   include("SimpleGraphics.jl")
   include("SimplePlot.jl")
   include("Contour.jl")
   include("LinearPredict.jl")
   include("MapPlot.jl")
   include("PauliMatrix.jl")
   include("SpinSim.jl")
   include("Peaks.jl")
   include("HMDB.jl")

   include("Bruker.jl")
   include("Varian.jl")

  function __init__()
    println("Module NMR v0.6.5\n(c)mu 2018-2021");
  end

end


text:
using Primes

function find_tiles(layers::Int, target::Int)
    prime_list = primes(Int(floor(sqrt(12 * layers + 5))))

    count = 0
    for n in 1:layers
        if is_prime(6 * n + 1, prime_list) && is_prime(6 * n - 1, prime_list) && is_prime(12 * n + 5, prime_list)
            count += 1
            if count == target
                println(3 * n^2 - 3 * n + 2)
            end
        end

        if is_prime(6 * n + 5, prime_list) && is_prime(6 * n - 1, prime_list) && is_prime(12 * n - 7, prime_list)
            count += 1
            if count == target
                println(3 * n^2 + 3 * n + 1)
            end
        end
    end
end

function is_prime(n::Int, prime_list::Array)
    root = sqrt(n)
    for prime in prime_list
        if prime > root
            return true
        end

        if n % prime == 0
            return false
        end
    end
    return true
end

function main()
    layers = 70000
    target = 2000
    find_tiles(layers, target)
end

main()

max_stars_repo_name: HarrisonGreen/Project-Euler-Solutions

max_issues_repo_name: HarrisonGreen/Project-Euler-Solutions

max_forks_repo_name: HarrisonGreen/Project-Euler-Solutions


text:
# utility/z-test.jl

using Test: @testset

list = [
"downsample.jl"
"eql_root.jl"
"interp1.jl"
"jinc.jl"
"map_many.jl"
"mask.jl"
"max_percent_diff.jl"
"ndgrid.jl"
"rect.jl"
"reverser.jl"
"rmsd100.jl"
]

for file in list
	@testset "$file" begin
		include(file)
	end
end


text:
using SpecialFunctions, StatsFuns

@testset "logmvgamma" begin
    @testset "type behavior" for eltya in (Float32, Float64)
        p = rand(1:50)
        a = rand(eltya)
        # add p since loggamma is only define for positive arguments
        @test typeof(logmvgamma(p, a + p)) == eltya
    end

    @testset "consistent with loggamma" for eltya in (Float32, Float64)
        #  Γ₁(a) = Γ(a), Γ₂(a) = √π Γ(a) Γ(a - 0.5), etc
        a = rand(eltya) + 1 # add one since loggamma is only define for positive arguments
        @test logmvgamma(1, a) ≈ loggamma(a)
        @test logmvgamma(2, a) ≈    eltya(0.5logπ) + loggamma(a) + loggamma(a - eltya(0.5))
        @test logmvgamma(3, a) ≈ eltya((3/2)*logπ) + loggamma(a) + loggamma(a - eltya(0.5)) + loggamma(a - one(a))
    end

    @testset "consistent with itself" for eltya in (Float32, Float64)
        #  Γᵢ(a) = (π^{i-1/2}) Γ(a) Γᵢ₋₁(a - 0.5)
        for p in 1:50
            a = rand(eltya) + p # add p since loggamma is only define for positive arguments
            @test logmvgamma(p, a) ≈ eltya((p/2-1/2)*logπ) + loggamma(a) + logmvgamma(p - 1, a - eltya(0.5))
        end
    end
end


text:
using CUSOLVER
using CUDArt
using Base.Test

m = 15
n = 10
k = 1

@testset for elty in [Float32, Float64, Complex64, Complex128]
    @testset "potrf!" begin
        A    = rand(elty,n,n)
        A    = A*A' #posdef
        d_A  = CudaArray(A)
        d_A  = CUSOLVER.potrf!('U',d_A)
        h_A  = to_host(d_A)
        cA,_ = LAPACK.potrf!('U',A)
        @test h_A ≈ cA
        A    = rand(elty,m,n)
        d_A  = CudaArray(A)
        @test_throws DimensionMismatch CUSOLVER.potrf!('U',d_A)
        A    = zeros(elty,n,n)
        d_A  = CudaArray(A)
        @test_throws Base.LinAlg.SingularException CUSOLVER.potrf!('U',d_A)
    end
    @testset "potrs!" begin
        A     = rand(elty,n,n)
        A     = A*A' #posdef
        d_A   = CudaArray(A)
        d_A   = CUSOLVER.potrf!('U',d_A)
        h_A   = to_host(d_A)
        B     = rand(elty,n,n)
        d_B   = CudaArray(B)
        d_B   = CUSOLVER.potrs!('U',d_A,d_B)
        h_B   = to_host(d_B)
        @test h_B ≈ LAPACK.potrs!('U',h_A,B)
        A    = rand(elty,m,n)
        d_A  = CudaArray(A)
        @test_throws DimensionMismatch CUSOLVER.potrs!('U',d_A,d_B)
        A    = rand(elty,n,n)
        d_A  = CudaArray(A)
        B     = rand(elty,m,m)
        d_B   = CudaArray(B)
        @test_throws DimensionMismatch CUSOLVER.potrs!('U',d_A,d_B)
    end

    @testset "getrf!" begin
        A          = rand(elty,m,n)
        d_A        = CudaArray(A)
        d_A,d_ipiv = CUSOLVER.getrf!(d_A)
        h_A        = to_host(d_A)
        h_ipiv     = to_host(d_ipiv)
        alu        = Base.LinAlg.LU(h_A, convert(Vector{Int},h_ipiv), zero(Int))
        @test A ≈ full(alu)
        A    = zeros(elty,n,n)
        d_A  = CudaArray(A)
        @test_throws Base.LinAlg.SingularException CUSOLVER.getrf!(d_A)
    end

    @testset "getrs!" begin
        A          = rand(elty,n,n)
        d_A        = CudaArray(A)
        d_A,d_ipiv = CUSOLVER.getrf!(d_A)
        B          = rand(elty,n,n)
        d_B        = CudaArray(B)
        d_B        = CUSOLVER.getrs!('N',d_A,d_ipiv,d_B)
        h_B        = to_host(d_B)
        @test h_B  ≈ A\B
        A          = rand(elty,m,n)
        d_A        = CudaArray(A)
        @test_throws DimensionMismatch CUSOLVER.getrs!('N',d_A,d_ipiv,d_B)
        A          = rand(elty,n,n)
        d_A        = CudaArray(A)
        B          = rand(elty,m,n)
        d_B        = CudaArray(B)
        @test_throws DimensionMismatch CUSOLVER.getrs!('N',d_A,d_ipiv,d_B)
    end

    @testset "geqrf!" begin
        A         = rand(elty,m,n)
        d_A       = CudaArray(A)
        d_A,d_tau = CUSOLVER.geqrf!(d_A)
        h_A       = to_host(d_A)
        h_tau     = to_host(d_tau)
        qra       = Base.LinAlg.QR(h_A, h_tau)
        @test A ≈ full(qra)
    end

    @testset "ormqr!" begin
        A         = rand(elty,n,n)
        d_A       = CudaArray(A)
        d_A,d_tau = CUSOLVER.geqrf!(d_A)
        h_A       = to_host(d_A)
        h_tau     = to_host(d_tau)
        qra       = Base.LinAlg.QR(h_A, h_tau)
        B         = rand(elty,n,n)
        d_B       = CudaArray(B)
        d_B       = CUSOLVER.ormqr!('L','N',d_A,d_tau,d_B)
        h_B       = to_host(d_B)
        @test h_B ≈ qra[:Q]*B
    end

    @testset "sytrf!" begin
        A          = rand(elty,n,n)
        A          = A + A' #symmetric
        d_A        = CudaArray(A)
        d_A,d_ipiv = CUSOLVER.sytrf!('U',d_A)
        h_A        = to_host(d_A)
        h_ipiv     = to_host(d_ipiv)
        A, ipiv    = LAPACK.sytrf!('U',A)
        @test ipiv == h_ipiv
        @test A ≈ h_A
        A    = rand(elty,m,n)
        d_A  = CudaArray(A)
        @test_throws DimensionMismatch CUSOLVER.sytrf!('U',d_A)
        A    = zeros(elty,n,n)
        d_A  = CudaArray(A)
        @test_throws Base.LinAlg.SingularException CUSOLVER.sytrf!('U',d_A)
    end

    @testset "gebrd!" begin
        A                             = rand(elty,m,n)
        d_A                           = CudaArray(A)
        d_A, d_D, d_E, d_TAUQ, d_TAUP = CUSOLVER.gebrd!(d_A)
        h_A                           = to_host(d_A)
        h_D                           = to_host(d_D)
        h_E                           = to_host(d_E)
        h_TAUQ                        = to_host(d_TAUQ)
        h_TAUP                        = to_host(d_TAUP)
        A,d,e,q,p                     = LAPACK.gebrd!(A)
        #@test A ≈ h_A
        @test d ≈ h_D
        @test e ≈ h_E
        @test q ≈ h_TAUQ
        @test p ≈ h_TAUP
    end

    @testset "gesvd!" begin
        A              = rand(elty,m,n)
        d_A            = CudaArray(A)
        d_U, d_S, d_Vt = CUSOLVER.gesvd!('A','A',d_A)
        h_S            = to_host(d_S)
        h_U            = to_host(d_U)
        h_Vt           = to_host(d_Vt)
        svda           = svdfact(A,thin=false)
        @test h_U ≈ svda[:U]
        @test h_S ≈ svdvals(A)
        @test h_Vt ≈ svda[:Vt]
    end
end


text:
include("common.jl")

fp = open(RESULT_PATH, "w")
println(fp, "family,instance,solver,formulation,optimal,solve_time,true_primal,model_primal,dual")

for family in FAMILIES, instance in readdir(joinpath(ROOT_DIR, family))
    Q, c = MIPforQP.BoxQP.parse_box_qp(joinpath(ROOT_DIR, family, instance))
    primal_model = MIPforQP.BoxQP.formulate_box_qp(Q, c, primal_factory, epigraph=EPIGRAPH)
    optimize!(primal_model)
    @assert primal_status(primal_model) == MOI.FEASIBLE_POINT
    cutoff = JuMP.objective_value(primal_model)
    for (solver, method) in CONFIGURATIONS
        factory, bridge = FACTORIES_AND_BRIDGES[(solver, method)]
        model = MIPforQP.BoxQP.formulate_box_qp(Q, c, factory(cutoff), epigraph=EPIGRAPH)
        bridge(model)
        JuMP.optimize!(model)
        true_pb = Inf
        model_pb = Inf
        if primal_status(model) == MOI.FEASIBLE_POINT
            model_pb = JuMP.objective_value(model)
            x = [JuMP.variable_by_name(model, "x[$i]") for i in 1:length(c)]
            @assert isa(x, Vector{JuMP.VariableRef})
            x_val = JuMP.value.(x)
            true_pb = -0.5dot(x_val, Q * x_val) - dot(c, x_val)
        end
        db = JuMP.objective_bound(model)
        st = JuMP.solve_time(model)
        opt = JuMP.termination_status(model) == MOI.OPTIMAL
        println(fp, join([family, instance, solver, method, opt, st, true_pb, model_pb, db], ","))
        flush(fp)
    end
end

close(fp)

max_stars_repo_name: joehuchette/compact-mips-for-qp

max_issues_repo_name: joehuchette/compact-mips-for-qp

max_forks_repo_name: joehuchette/compact-mips-for-qp


text:
module JungleHelperLantern

using ..Ahorn, Maple

@mapdef Entity "JungleHelper/Lantern" Lantern(x::Integer, y::Integer, sprite::String="", onlyIfMaddyNotHolding::Bool=false)

const placements = Ahorn.PlacementDict(
    "Lantern (Jungle Helper)" => Ahorn.EntityPlacement(
        Lantern
    )
)

sprite = "JungleHelper/Lantern/LanternEntity/lantern_00"

function Ahorn.selection(entity::Lantern)
    x, y = Ahorn.position(entity)
    return Ahorn.getSpriteRectangle(sprite, x, y + 5)
end

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Lantern, room::Maple.Room) = Ahorn.drawSprite(ctx, sprite, 0, 5)

end

max_stars_repo_name: Jackal-Celeste/JungleHelper

max_issues_repo_name: Jackal-Celeste/JungleHelper

max_forks_repo_name: Jackal-Celeste/JungleHelper


text:
function Polygon{T}(
    domain::Polygon{T},
    parent::ArbField = domain.parent,
) where {T<:Union{arb,fmpq}}
    return Polygon(copy(domain.angles), copy(domain.vertices), parent)
end

function Polygon(
    angles::AbstractVector{T},
    vertices::AbstractVector,
    parent::ArbField,
) where {T}
    @assert length(angles) == length(vertices)
    if T == fmpq || T <: Rational
        S = fmpq
        angles = fmpq.(angles)
    else
        S = arb
        angles = parent.(angles)
    end
    vertices = [SVector{2,arb}(parent.(vertex)) for vertex in vertices]
    return Polygon{S}(angles, vertices, parent)
end

function Base.show(io::IO, domain::Polygon)
    print(
        io,
        "Polygon with $(length(boundaries(domain))) sides and $(domain.parent.prec) bits of precision",
    )
end

boundaries(domain::Polygon) = 1:length(domain.angles)

angle(domain::Polygon{fmpq}, i::Integer) = domain.parent(π) * domain.angles[i]
angle(domain::Polygon{arb}, i::Integer) = domain.angles[i]

angledivπ(domain::Polygon{fmpq}, i::Integer) = domain.angles[i]
angledivπ(domain::Polygon{arb}, i::Integer) = domain.angles[i] / domain.parent(π)

angles(domain::Polygon) = [angle(domain, i) for i in boundaries(domain)]
anglesdivπ(domain::Polygon) = [angledivπ(domain, i) for i in boundaries(domain)]

vertex(domain::Polygon, i::Integer) = domain.vertices[i]
vertices(domain::Polygon) = domain.vertices

function area(domain::Polygon)
    # https://en.wikipedia.org/wiki/Polygon#Area
    A = domain.parent(0)
    vs = vertices(domain)
    for i in eachindex(vs)
        A +=
            vs[i][1] * vs[mod1(i + 1, length(vs))][2] -
            vs[mod1(i + 1, length(vs))][1] * vs[i][2]
    end

    return abs(A) / 2
end
center(domain::Polygon) = sum(vertices(domain)) / length(boundaries(domain))

# TODO: Add this
function Base.in(xy, domain::Polygon)
    vs = vertices(domain)

    isleft(p0, p1, p2) =
        (p1[1] - p0[1]) * (p2[2] - p0[2]) - (p2[1] - p0[1]) * (p1[2] - p0[2])

    wn = 0
    for i in eachindex(vs)
        vi = vs[i]
        vip1 = vs[mod1(i + 1, length(vs))]
        if vi[2] <= xy[2]
            if !(vip1[2] <= xy[2])
                if isleft(vi, vip1, xy) > 0
                    wn += 1
                end
            end
        else
            if vip1[2] <= xy[2]
                if isleft(vi, vip1, xy) < 0
                    wn -= 1
                end
            end
        end
    end

    return !iszero(wn)
end

function boundary_parameterization(t, domain::Polygon, i::Integer)
    i ∈ boundaries(domain) ||
        throw(ArgumentError("attempt to get vertex $i from a $(typeof(domain))"))
    v = vertex(domain, mod1(i, length(boundaries(domain))))
    w = vertex(domain, mod1(i + 1, length(boundaries(domain))))

    return v .+ t .* (w - v)
end

function boundary_points(domain::Polygon, i::Integer, n::Integer; distribution = :chebyshev)
    i ∈ boundaries(domain) ||
        throw(ArgumentError("attempt to get vertex $i from a $(typeof(domain))"))

    v = vertex(domain, mod1(i, length(boundaries(domain))))
    w = vertex(domain, mod1(i + 1, length(boundaries(domain))))

    if distribution == :linear
        f = j -> domain.parent(j // (n + 1))
    elseif distribution == :chebyshev
        f = j -> (cospi(domain.parent((2j - 1) // 2n)) + 1) / 2
    elseif distribution == :exponential
        m = (n + 1) // 2
        f = j -> begin
            d = exp(-domain.parent(m - j))
            if j < m
                return d / 2
            elseif j > m
                return 1 - inv(d) / 2
            else
                return domain.parent(0.5)
            end
        end
    elseif distribution == :root_exponential
        m = (n + 1) // 2
        f =
            j -> begin
                d = exp(-4(sqrt(domain.parent(m)) - sqrt(domain.parent(min(j, n - j + 1)))))
                if j < m
                    return d / 2
                elseif j > m
                    return 1 - d / 2
                else
                    return domain.parent(0.5)
                end
            end
    else
        throw(ArgumentError("no distribution named $distribution"))
    end

    points = Vector{SVector{2,arb}}(undef, n)
    for j = 1:n
        t = f(j)
        points[j] = v .+ t .* (w - v)
    end

    return points, fill(i, n)
end

# TODO: Fix this
function interior_points(domain::Polygon, n::Integer; rng = MersenneTwister(42))
    points = Vector{SVector{2,arb}}(undef, n)
    xmin, xmax = extrema(getindex.(vertices(domain), 1))
    ymin, ymax = extrema(getindex.(vertices(domain), 2))
    for i = 1:n
        pt = SVector(xmin + rand(rng) * (xmax - xmin), ymin + rand(rng) * (ymax - ymin))
        while !(pt ∈ domain)
            pt = SVector(xmin + rand(rng) * (xmax - xmin), ymin + rand(rng) * (ymax - ymin))
        end

        points[i] = pt
    end

    return points
end


text:
# Hack: https://github.com/timholy/CodeTracking.jl/issues/45
#==
function Revise.get_def(method::Method; modified_files=Revise.revision_queue)
    yield()   # magic bug fix for the OSX test failures. TODO: figure out why this works (prob. Julia bug)
    filename = Revise.fixpath(String(method.file))
    if startswith(filename, "REPL[")
        isdefined(Base, :active_repl) || return false
        fi = Revise.add_definitions_from_repl(filename)
        hassig = false
        for (mod, exs) in fi.modexsigs
            for sigs in values(exs)
                hassig |= !isempty(sigs)
            end
        end
        return hassig
    end
    id = Revise.get_tracked_id(method.module; modified_files=modified_files)
    id === nothing && return false
    pkgdata = Revise.pkgdatas[id]
    filename = relpath(filename, pkgdata)
    if Revise.hasfile(pkgdata, filename)
        def = Revise.get_def(method, pkgdata, filename)
        def !== nothing && return true
    end
    # Lookup can fail for macro-defined methods, see https://github.com/JuliaLang/julia/issues/31197
    # We need to find the right file.
    if method.module == Base || method.module == Core || method.module == Core.Compiler
        @warn "skipping $method to avoid parsing too much code"
        CodeTracking.method_info[method.sig] = missing
        return false
    end
    parentfile, included_files = Revise.modulefiles(method.module)
    if parentfile !== nothing
        def = Revise.get_def(method, pkgdata, relpath(parentfile, pkgdata))
        def !== nothing && return true
        for modulefile in included_files
            def = Revise.get_def(method, pkgdata, relpath(modulefile, pkgdata))
            def !== nothing && return true
        end
    end
    # As a last resort, try every file in the package
    for file in Revise.srcfiles(pkgdata)
        def = Revise.get_def(method, pkgdata, file)
        def !== nothing && return true
    end
    @warn "$(method.sig) was not found"
    # So that we don't call it again, store missingness info in CodeTracking
    CodeTracking.method_info[method.sig] = missing
    return false
end
==#


text:
export VPGPolicy

"""
Vanilla Policy Gradient

VPGPolicy(;kwargs)

# Keyword arguments
- `approximator`,
- `baseline`,
- `dist`, distribution function of the action
- `γ`, discount factor
- `α_θ`, step size of policy parameter
- `α_w`, step size of baseline parameter
- `batch_size`,
- `rng`,
- `loss`,
- `baseline_loss`,

if the action space is continuous,
then the env should transform the action value, (such as using tanh),
in order to make sure low ≤ value ≤ high
"""
Base.@kwdef mutable struct VPGPolicy{
    A<:NeuralNetworkApproximator,
    B<:Union{NeuralNetworkApproximator,Nothing},
    S,
    R<:AbstractRNG,
} <: AbstractPolicy
    approximator::A
    baseline::B = nothing
    action_space::S
    dist::Any
    γ::Float32 = 0.99f0 # discount factor
    α_θ = 1.0f0 # step size of policy
    α_w = 1.0f0 # step size of baseline
    batch_size::Int = 1024
    rng::R = Random.GLOBAL_RNG
    loss::Float32 = 0.0f0
    baseline_loss::Float32 = 0.0f0
end

"""
About continuous action space, see
* [Diagonal Gaussian Policies](https://spinningup.openai.com/en/latest/spinningup/rl_intro.html#stochastic-policies
* [Clipped Action Policy Gradient](https://arxiv.org/pdf/1802.07564.pdf)
"""

function (π::VPGPolicy)(env::AbstractEnv)
    to_dev(x) = send_to_device(device(π.approximator), x)

    logits = env |> state |> to_dev |> π.approximator |> send_to_host

    if π.action_space isa AbstractVector
        dist = logits |> softmax |> π.dist
        action = π.action_space[rand(π.rng, dist)]
    elseif π.action_space isa Interval
        dist = π.dist.(logits...)
        action = rand.(π.rng, dist)[1]
    else
        error("not implemented")
    end
    action
end

function (π::VPGPolicy)(env::MultiThreadEnv)
    error("not implemented")
    # TODO: can PG support multi env? PG only get updated at the end of an episode.
end

function RLBase.update!(
    trajectory::ElasticSARTTrajectory,
    policy::VPGPolicy,
    env::AbstractEnv,
    ::PreActStage,
    action,
)
    push!(trajectory[:state], state(env))
    push!(trajectory[:action], action)
end

function RLBase.update!(
    t::ElasticSARTTrajectory,
    ::VPGPolicy,
    ::AbstractEnv,
    ::PreEpisodeStage,
)
    empty!(t)
end

RLBase.update!(::VPGPolicy, ::ElasticSARTTrajectory, ::AbstractEnv, ::PreActStage) = nothing

function RLBase.update!(
    π::VPGPolicy,
    traj::ElasticSARTTrajectory,
    ::AbstractEnv,
    ::PostEpisodeStage,
)
    model = π.approximator
    to_dev(x) = send_to_device(device(model), x)

    states = traj[:state]
    actions = traj[:action] |> Array # need to convert ElasticArray to Array, or code will fail on gpu. `log_prob[CartesianIndex.(A, 1:length(A))`
    gains = traj[:reward] |> x -> discount_rewards(x, π.γ)

    for idx in Iterators.partition(shuffle(1:length(traj[:terminal])), π.batch_size)
        S = select_last_dim(states, idx) |> Array |> to_dev
        A = actions[idx]
        G = gains[idx] |> x -> Flux.unsqueeze(x, 1) |> to_dev
        # gains is a 1 column array, but the output of flux model is 1 row, n_batch columns array. so unsqueeze it.

        if π.baseline isa NeuralNetworkApproximator
            gs = gradient(Flux.params(π.baseline)) do
                δ = G - π.baseline(S)
                loss = mean(δ .^ 2) * π.α_w # mse
                ignore() do
                    π.baseline_loss = loss
                end
                loss
            end
            update!(π.baseline, gs)
        elseif π.baseline isa Nothing
            # Normalization. See
            # (http://rail.eecs.berkeley.edu/deeprlcourse-fa17/f17docs/hw2_final.pdf)
            # (https://web.stanford.edu/class/cs234/assignment3/solution.pdf)
            # normalise should not be used with baseline. or the loss of the policy will be too small.
            δ = G |> x -> normalise(x; dims = 2)
        end

        gs = gradient(Flux.params(model)) do
            if π.action_space isa AbstractVector
                log_prob = S |> model |> logsoftmax
                log_probₐ = log_prob[CartesianIndex.(A, 1:length(A))]
            elseif π.action_space isa Interval
                dist = π.dist.(model(S)...) # TODO: this part does not work on GPU. See: https://github.com/JuliaStats/Distributions.jl/issues/1183 .
                log_probₐ = logpdf.(dist, A)
            end
            loss = -mean(log_probₐ .* δ) * π.α_θ
            ignore() do
                π.loss = loss
            end
            loss
        end
        update!(model, gs)
    end
end


text:
# Special purpose updating for operations which can be done in-place. This function is
# just internal and free-form; it is not a method of `accumulate!` directly as it does
# not adhere to the expected method signature form, i.e. `accumulate!(value, rule, args)`.
# Instead it's `_update!(old, new, extrastuff...)` and is not specific to any particular
# rule.

_update!(x, y) = x + y
_update!(x::Array{T,N}, y::AbstractArray{T,N}) where {T,N} = x .+= y

_update!(x, ::Zero) = x
_update!(::Zero, y) = y
_update!(::Zero, ::Zero) = Zero()

function _update!(x::NamedTuple{Ns}, y::NamedTuple{Ns}) where Ns
    return NamedTuple{Ns}(map(p->_update!(getproperty(x, p), getproperty(y, p)), Ns))
end

function _update!(x::NamedTuple, y, p::Symbol)
    new = NamedTuple{(p,)}((_update!(getproperty(x, p), y),))
    return merge(x, new)
end

function _update!(x::NamedTuple{Ns}, y::NamedTuple{Ns}, p::Symbol) where Ns
    return _update!(x, getproperty(y, p), p)
end

function _checked_rrule(f, args...; kwargs...)
    r = rrule(f, args...; kwargs...)
    r isa Nothing && _throw_checked_rrule_error(f, args...; kwargs...)
    return r
end

@noinline function _throw_checked_rrule_error(f, args...; kwargs...)
    io = IOBuffer()
    print(io, "can't differentiate `", f, '(')
    join(io, map(arg->string("::", typeof(arg)), args), ", ")
    if !isempty(kwargs)
        print(io, ";")
        join(io, map(((k, v),)->string(k, "=", v), kwargs), ", ")
    end
    print(io, ")`; no matching `rrule` is defined")
    throw(ArgumentError(String(take!(io))))
end


text:
include("error.jl")

const GDALVERSION = Ref{VersionNumber}()
const GDAL_DATA = Ref{String}()
const PROJ_LIB = Ref{String}()

function __init__()
    # register custom error handler
    funcptr = @cfunction(gdaljl_errorhandler, Ptr{Cvoid}, (CPLErr, Cint, Cstring))
    cplseterrorhandler(funcptr)

    # get GDAL version number
    versionstring = gdalversioninfo("RELEASE_NAME")
    GDALVERSION[] = VersionNumber(versionstring)

    # set GDAL_DATA location, this overrides setting the environment variable
    GDAL_DATA[] = joinpath(GDAL_jll.artifact_dir, "share", "gdal")
    cplsetconfigoption("GDAL_DATA", GDAL_DATA[])

    # set path to CA certificates
    ca_path = @static if VERSION >= v"1.6"
        ca_roots()
    else
        ca_roots_path()
    end
    if ca_path !== nothing
        cplsetconfigoption("CURL_CA_BUNDLE", ca_path)
    end

    # set PROJ_LIB location, this overrides setting the environment variable
    PROJ_LIB[] = joinpath(PROJ_jll.artifact_dir, "share", "proj")
    osrsetprojsearchpaths([PROJ_LIB[]])

    # register all known configured GDAL drivers
    gdalallregister()
end


text:
module CMAES
##

    using Optim, Distributions, Parameters, PDMats
    using LinearAlgebra

    import Optim: AbstractObjective, ZerothOrderOptimizer, ZerothOrderState, initial_state, update_state!,
                  trace!, assess_convergence, AbstractOptimizerState, update!, value, value!,
                  pick_best_x, pick_best_f, f_abschange, x_abschange

    const MatTypes = Union{Matrix,Diagonal}

    @static if VERSION < v"0.7.0-DEV.3510"
        import Base.LinAlg: copy_oftype, checksquare
    else
        using LinearAlgebra
        import LinearAlgebra: copy_oftype, checksquare
    end

    import Base.(^)
    function (^)(A::Diagonal{T}, p::Real) where T

        n = checksquare(A)
        TT = Base.promote_op(^, T, typeof(p))
        retmat = copy_oftype(A, TT)
        for i in 1:n
            retmat[i, i] = retmat[i, i] ^ p
        end
        return retmat
    end

    abstract type CovMatrix end

    struct CMA{K<:MatTypes, CM <: CovMatrix} <: Optim.ZerothOrderOptimizer
        λ::Int
        σ::Float64
    end

    init_matrix(::Type{Matrix},n) = Matrix(Diagonal(ones(n)))
    init_matrix(::Type{Diagonal},n) = Diagonal(ones(n))

    mutable struct CMAState{T <: AbstractFloat, CM <: CovMatrix} <: Optim.ZerothOrderState
        x::Array{T,1}
        x_previous::Array{T,1} #this is used by default convergence test
        f_x::T
        f_x_previous::T
        #
        w::Vector{T}
        n::Int
        μ_w::T
        c_σ::T
        d_σ::T
        c_c::T
        #
        p_σ::Vector{T}
        p_c::Vector{T}
        cov_mat::CM
        m::Vector{T}
        chi_D::T
        σ::T
        σ_0::T
        λ::Int
        t::Int
        xs::Vector{Vector{T}}
    end

    ∑(x) = sum(x)
    weights(μ) = log(μ + 1/2) .- log.(1:μ) #[(log(μ+1)-log(i)) for i=1:μ] / ∑( log(μ+1)-log(j) for j=1:μ )

    include("FullCovMatrix.jl")
    include("BD_CovMatrix.jl")

    CMA(λ::Int, σ::Real) = CMA{Matrix}(λ,σ)
    function CMA(σ::Real; dimension=0)
         @assert dimension > 0
         CMA{Matrix,BD_CovMatrix}(4+floor(Int,3*log(dimension)), σ)
    end

    function CMA(::Type{T},σ::Real; dimension=0) where T <: CovMatrix
         @assert dimension > 0
         CMA{Matrix,T}(4+floor(Int,3*log(dimension)), σ)
    end

    ## Dict support (view previous commit)

    npara(xinit::Vector) = length(xinit)

    ##

    function initial_state(method::CMA{K,CM}, options, d, xinit) where {K <: MatTypes, CM <: CovMatrix}

        λ = method.λ
        n = npara(xinit)

        μ = floor(Int,λ/2)
        w = weights(μ)
        normalize!(w, 1)

        m = xinit
        n, μ_w, c_σ, d_σ, c_c = init_constants(method,n,λ,w,μ)
        p_σ, p_c = zeros(n), ones(n)
        cov_mat = CM(K, n, μ_w)

        chi_D = √n*(1-1/(4n) + 1/(21n^2))
        xs = [zeros(n) for i=1:μ]
        σ_0 = copy(method.σ)

        CMAState(
            copy(xinit),copy(xinit),Inf,Inf,w,
            n, μ_w, c_σ, d_σ, c_c,
            p_σ, p_c,cov_mat,m,chi_D,method.σ,σ_0,λ,1,xs,
        )
    end

    function trace!(tr, d, state, iteration, method::CMA, options)
        dt = Dict()
        if options.extended_trace
            dt["x"] = state.xs
        end
        g_norm = 0.0
        update!(tr,
        iteration,
        state.f_x,
        g_norm,
        dt,
        options.store_trace,
        options.show_trace,
        options.show_every,
        options.callback)
    end

    function eval_objfun(T, λ, d, x)
        fx = zeros(T,λ)
        for i=1:λ
            value!(d,x[i])
            fx[i] = value(d)
            isnan(fx[i]) && error("Function returned NaN for input $(x[i])")
        end
        fx
    end

    function update_H_σ(::Type{CM},p_σ, c_σ, t, n) where CM <: CovMatrix
        norm(p_σ) < √(1-(1-c_σ)^(2*(t+1)))*(1.4 + 2/(n+1)) ? 1.0 : 0.0
    end

    function update_σ(::Type{CM},σ, c_σ, d_σ, p_σ, chi_D) where CM <: CovMatrix
    σ * exp(c_σ/d_σ*(norm(p_σ)/chi_D -1))
    end

    ⊗(::Type{Matrix{T}},x)   where T = x * x'
    ⊗(::Type{Diagonal{T,Vector{T}}},x) where T = Diagonal(x.^2)

    function init_constants(method::CMA{K},n,λ,w,μ) where K <: MatTypes

        μ_w = 1 / sum(w.^2)
        c_σ =  (μ_w + 2) / (n + μ_w + 5)

        d_σ = 1 + c_σ + 2*max(0, √((μ_w-1)/(n+1)) -1)

        c_c = (4 + μ_w/n)/(n + 4 + 2μ_w/n) #Benchmarking a BI-Population CMA-ES on the BBOB-2009 Function Testbed
#        c_c = 4 / (n + 4)

#        if K <: Diagonal #A Simple Modification in CMA-ES Achieving Linear Time and Space Complexity
#            c_cov = min(1, c_cov*(n+2)/3) # I need to clamp here for low n
#        end

        n, μ_w, c_σ, d_σ, c_c
    end

#(w, μeff, cc, cσ, c1, cμ, dσ) = ([0.637043, 0.28457, 0.0783872], 2.0286114646100617, 0.5714285714285714, 0.5017818438926943, 0.09747248265066792, 0.038593139193450914, 1.5017818438926942)([0.637043, 0.28457, 0.0783872], 2.0286114646100617, 0.5714285714285714, 0.5017818438926943, 0.09747248265066792, 0.038593139193450914, 1.5017818438926942)

    function assess_convergence(state::CMAState, d, options)

        x_converged, f_converged, g_converged, converged, f_increased = Optim.default_convergence_assessment(state, d, options)
        f_increased = false #disable this one, since error can increase with CMAES

        @unpack n, λ, p_c, σ, σ_0,t, cov_mat = state
        C,D = cov_matrix(cov_mat), eigvals(cov_mat)

        MaxIter = false# t > 100 + 50*(n+3)*2/√(λ)

        TolX = t > 1 && all(p_c * σ/σ_0 .< 1e-16) && all(sqrt.(diag(C)) * σ/σ_0  .< 1e-16)
        TolUpSigma = σ/σ_0 > 10^20 * √(maximum(D))
        ConditionCov = abs(maximum(D)) / abs(minimum(D)) > 10^14

        converged = MaxIter || TolX || TolUpSigma || ConditionCov

        x_converged, f_converged, g_converged, converged, f_increased
    end

    f_abschange(d::AbstractObjective, state::CMAState) = f_abschange(state.f_x, state.f_x_previous)
    x_abschange(state::CMAState) = x_abschange(state.x, state.x_previous)

    pick_best_f(f_increased, state::CMAState, d) = state.f_x
    pick_best_x(f_increased, state::CMAState) = state.x

##
end

##

##