text:
@testset "AD tests" begin
    ddim = 3
    dmean = zeros(ddim)
    cov_mat = Matrix{Float64}(I, ddim, ddim)
    cov_vec = ones(ddim)
    cov_num = 1.0
    norm_val_vec = ones(ddim)
    norm_val_mat = ones(ddim, 2)
    alpha = ones(4)
    dir_val_vec = fill(0.25, 4)
    dir_val_mat = fill(0.25, 4, 2)
    beta_mat = rand(MatrixBeta(ddim, ddim, ddim))
    tested = []
    function test_info(name)
        if !(name in tested)
            push!(tested, name)
            @info("Testing: $(name)")
        end
    end

    function filldist_expr(sym, n = 1; dims = (1,))
        x = [gensym(:x) for i in 1:n]
        :(($(x...),) -> filldist($sym($(x...),), $(dims...)))
    end
    function arraydist_expr(sym, n = 1; dims = (1,))
        x = [gensym(:x) for i in 1:n]
        :(($(x...),) -> arraydist(reshape(fill($sym($(x...),), prod($(dims))), $(dims...))))
    end
    function filldist_spec(dist::DistSpec; disttype = :uni, n = 2, d = 1)
        @assert disttype in (:uni, :multi)
        disttype == :uni && dist.x isa Vector && return
        disttype == :multi && dist.x isa Matrix && return
        # Broken
        (dist.name in (:VonMises, :TriangularDist)) && return
        (dist.name == :Weibull && n isa Tuple) && return
        # Tests are failing for matrix covariance vectorized MvNormal
        # Note the 2 MvNormal cases in broken_multi_cont_dists
        if dist.name in (:MvNormal, :TuringMvNormal, :MvLogNormal, :TuringMvNormal) &&
            any(x -> isa(x, Matrix), dist.θ)
            return
        end
        if disttype == :multi
            name = filldist_expr(dist.name, length(dist.θ); dims = (n,))
            x = hcat(fill(dist.x, n)...,)
            x = d == 1 ? x : fill(x, d)
        else
            if n isa Tuple
                name = filldist_expr(dist.name, length(dist.θ); dims = n)
                x = fill(dist.x, n...,)
                x = d == 1 ? x : fill(x, d)
            else # n number
                name = filldist_expr(dist.name, length(dist.θ); dims = (n,))
                x = fill(dist.x, n)
                x = d == 1 ? x : hcat(fill(x, d)...,)
            end
        end
        return DistSpec(name, dist.θ, x)
    end
    function arraydist_spec(dist::DistSpec; disttype = :uni, n = 2, d = 1)
        @assert disttype in (:uni, :multi)
        (disttype == :uni) && (dist.x isa Vector) && return
        (disttype == :multi) && (dist.x isa Matrix) && return
        # Broken
        (dist.name in (:VonMises, :TriangularDist)) && return
        (dist.name == :NormalInverseGaussian && n isa Tuple) && return
        if dist.name in (:MvNormal, :TuringMvNormal, :MvLogNormal, :TuringMvNormal) &&
            any(x -> isa(x, Matrix), dist.θ)
            return
        end
        if disttype == :multi
            name = arraydist_expr(dist.name, length(dist.θ); dims = (n,))
            x = hcat(fill(dist.x, n)...,)
            x = d == 1 ? x : fill(x, d)
        else
            if n isa Tuple
                name = arraydist_expr(dist.name, length(dist.θ); dims = n)
                x = fill(dist.x, n...,)
                x = d == 1 ? x : fill(x, d)
            else
                name = arraydist_expr(dist.name, length(dist.θ); dims = (n,))
                x = fill(dist.x, n)
                x = d == 1 ? x : hcat(fill(x, d)...,)
            end
        end
        return DistSpec(name, dist.θ, x)
    end

    uni_disc_dists = [
        DistSpec(:Bernoulli, (0.45,), 1),
        DistSpec(:Bernoulli, (0.45,), [1, 1]),
        DistSpec(:Bernoulli, (0.45,), 0),
        DistSpec(:Bernoulli, (0.45,), [0, 0]),

        DistSpec(:((a, b) -> BetaBinomial(10, a, b)), (2.0, 1.0), 5),
        DistSpec(:((a, b) -> BetaBinomial(10, a, b)), (2.0, 1.0), [5, 5]),

        DistSpec(:(p -> Binomial(10, p)), (0.5,), 5),
        DistSpec(:(p -> Binomial(10, p)), (0.5,), [5, 5]),

        DistSpec(:(p -> Categorical(p / sum(p))), ([0.45, 0.55],), 1),
        DistSpec(:(p -> Categorical(p / sum(p))), ([0.45, 0.55],), [1, 1]),

        DistSpec(:Geometric, (0.45,), 3),
        DistSpec(:Geometric, (0.45,), [3, 3]),

        DistSpec(:NegativeBinomial, (3.5, 0.5), 1),
        DistSpec(:NegativeBinomial, (3.5, 0.5), [1, 1]),

        DistSpec(:Poisson, (0.5,), 1),
        DistSpec(:Poisson, (0.5,), [1, 1]),

        DistSpec(:Skellam, (1.0, 2.0), -2),
        DistSpec(:Skellam, (1.0, 2.0), [-2, -2]),

        DistSpec(:PoissonBinomial, ([0.5, 0.5],), 0),
        DistSpec(:PoissonBinomial, ([0.5, 0.5],), [0, 0]),

        DistSpec(:TuringPoissonBinomial, ([0.5, 0.5],), 0),
        DistSpec(:TuringPoissonBinomial, ([0.5, 0.5],), [0, 0]),
    ]

    uni_cont_dists = [
        DistSpec(:Arcsine, (), 0.5),
        DistSpec(:Arcsine, (1.0,), 0.5),
        DistSpec(:Arcsine, (0.0, 2.0), 0.5),
        DistSpec(:Arcsine, (), [0.5]),
        #DistSpec(:Arcsine, (1.0,), [0.5]),
        #DistSpec(:Arcsine, (0.0, 2.0), [0.5]),

        DistSpec(:Beta, (), 0.5),
        DistSpec(:Beta, (1.0,), 0.5),
        DistSpec(:Beta, (1.0, 2.0), 0.5),
        #DistSpec(:Beta, (), [0.5]),
        #DistSpec(:Beta, (1.0,), [0.5]),
        #DistSpec(:Beta, (1.0, 2.0), [0.5]),

        DistSpec(:BetaPrime, (), 0.5),
        DistSpec(:BetaPrime, (1.0,), 0.5),
        DistSpec(:BetaPrime, (1.0, 2.0), 0.5),
        #DistSpec(:BetaPrime, (), [0.5]),
        #DistSpec(:BetaPrime, (1.0,), [0.5]),
        #DistSpec(:BetaPrime, (1.0, 2.0), [0.5]),

        DistSpec(:Biweight, (), 0.5),
        DistSpec(:Biweight, (1.0,), 0.5),
        DistSpec(:Biweight, (1.0, 2.0), 0.5),
        #DistSpec(:Biweight, (), [0.5]),
        #DistSpec(:Biweight, (1.0,), [0.5]),
        #DistSpec(:Biweight, (1.0, 2.0), [0.5]),

        DistSpec(:Cauchy, (), 0.5),
        DistSpec(:Cauchy, (1.0,), 0.5),
        DistSpec(:Cauchy, (1.0, 2.0), 0.5),
        #DistSpec(:Cauchy, (), [0.5]),
        #DistSpec(:Cauchy, (1.0,), [0.5]),
        #DistSpec(:Cauchy, (1.0, 2.0), [0.5]),

        DistSpec(:Chi, (1.0,), 0.5),
        #DistSpec(:Chi, (1.0,), [0.5]),

        DistSpec(:Chisq, (1.0,), 0.5),
        #DistSpec(:Chisq, (1.0,), [0.5]),

        DistSpec(:Cosine, (1.0, 1.0), 0.5),
        #DistSpec(:Cosine, (1.0, 1.0), [0.5]),

        DistSpec(:Epanechnikov, (1.0, 1.0), 0.5),
        #DistSpec(:Epanechnikov, (1.0, 1.0), [0.5]),

        DistSpec(:((s)->Erlang(1, s)), (1.0,), 0.5), # First arg is integer
        #DistSpec(:((s)->Erlang(1, s)), (1.0,), [0.5]),

        DistSpec(:Exponential, (1.0,), 0.5),
        #DistSpec(:Exponential, (1.0,), [0.5]),

        DistSpec(:FDist, (1.0, 1.0), 0.5),
        #DistSpec(:FDist, (1.0, 1.0), [0.5]),

        DistSpec(:Frechet, (), 0.5),
        DistSpec(:Frechet, (1.0,), 0.5),
        DistSpec(:Frechet, (1.0, 2.0), 0.5),
        #DistSpec(:Frechet, (), [0.5]),
        #DistSpec(:Frechet, (1.0,), [0.5]),
        #DistSpec(:Frechet, (1.0, 2.0), [0.5]),

        DistSpec(:Gamma, (), 0.5),
        DistSpec(:Gamma, (1.0,), 0.5),
        DistSpec(:Gamma, (1.0, 2.0), 0.5),
        #DistSpec(:Gamma, (), [0.5]),
        #DistSpec(:Gamma, (1.0,), [0.5]),
        #DistSpec(:Gamma, (1.0, 2.0), [0.5]),

        DistSpec(:GeneralizedExtremeValue, (1.0, 1.0, 1.0), 0.5),
        #DistSpec(:GeneralizedExtremeValue, (1.0, 1.0, 1.0), [0.5]),

        DistSpec(:GeneralizedPareto, (), 0.5),
        DistSpec(:GeneralizedPareto, (1.0, 2.0), 0.5),
        DistSpec(:GeneralizedPareto, (0.0, 2.0, 3.0), 0.5),
        #DistSpec(:GeneralizedPareto, (), [0.5]),
        #DistSpec(:GeneralizedPareto, (1.0, 2.0), [0.5]),
        #DistSpec(:GeneralizedPareto, (0.0, 2.0, 3.0), [0.5]),

        DistSpec(:Gumbel, (), 0.5),
        DistSpec(:Gumbel, (1.0,), 0.5),
        DistSpec(:Gumbel, (1.0, 2.0), 0.5),
        #DistSpec(:Gumbel, (), [0.5]),
        #DistSpec(:Gumbel, (1.0,), [0.5]),
        #DistSpec(:Gumbel, (1.0, 2.0), [0.5]),

        DistSpec(:InverseGamma, (), 0.5),
        DistSpec(:InverseGamma, (1.0,), 0.5),
        DistSpec(:InverseGamma, (1.0, 2.0), 0.5),
        #DistSpec(:InverseGamma, (), [0.5]),
        #DistSpec(:InverseGamma, (1.0,), [0.5]),
        #DistSpec(:InverseGamma, (1.0, 2.0), [0.5]),

        DistSpec(:InverseGaussian, (), 0.5),
        DistSpec(:InverseGaussian, (1.0,), 0.5),
        DistSpec(:InverseGaussian, (1.0, 2.0), 0.5),
        #DistSpec(:InverseGaussian, (), [0.5]),
        #DistSpec(:InverseGaussian, (1.0,), [0.5]),
        #DistSpec(:InverseGaussian, (1.0, 2.0), [0.5]),

        DistSpec(:Kolmogorov, (), 0.5),
        #DistSpec(:Kolmogorov, (), [0.5]),

        DistSpec(:Laplace, (), 0.5),
        DistSpec(:Laplace, (1.0,), 0.5),
        DistSpec(:Laplace, (1.0, 2.0), 0.5),
        #DistSpec(:Laplace, (), [0.5]),
        #DistSpec(:Laplace, (1.0,), [0.5]),
        #DistSpec(:Laplace, (1.0, 2.0), [0.5]),

        DistSpec(:Levy, (), 0.5),
        DistSpec(:Levy, (0.0,), 0.5),
        DistSpec(:Levy, (0.0, 2.0), 0.5),
        #DistSpec(:Levy, (), [0.5]),
        #DistSpec(:Levy, (0.0,), [0.5]),
        #DistSpec(:Levy, (0.0, 2.0), [0.5]),

        DistSpec(:((a, b) -> LocationScale(a, b, Normal())), (1.0, 2.0), 0.5),
        #DistSpec(:((a, b) -> LocationScale(a, b, Normal())), (1.0, 2.0), [0.5]),

        DistSpec(:Logistic, (), 0.5),
        DistSpec(:Logistic, (1.0,), 0.5),
        DistSpec(:Logistic, (1.0, 2.0), 0.5),
        #DistSpec(:Logistic, (), [0.5]),
        #DistSpec(:Logistic, (1.0,), [0.5]),
        #DistSpec(:Logistic, (1.0, 2.0), [0.5]),

        DistSpec(:LogitNormal, (), 0.5),
        DistSpec(:LogitNormal, (1.0,), 0.5),
        DistSpec(:LogitNormal, (1.0, 2.0), 0.5),
        #DistSpec(:LogitNormal, (), [0.5]),
        #DistSpec(:LogitNormal, (1.0,), [0.5]),
        #DistSpec(:LogitNormal, (1.0, 2.0), [0.5]),

        DistSpec(:LogNormal, (), 0.5),
        DistSpec(:LogNormal, (1.0,), 0.5),
        DistSpec(:LogNormal, (1.0, 2.0), 0.5),
        #DistSpec(:LogNormal, (), [0.5]),
        #DistSpec(:LogNormal, (1.0,), [0.5]),
        #DistSpec(:LogNormal, (1.0, 2.0), [0.5]),

        DistSpec(:Normal, (), 0.5),
        DistSpec(:Normal, (1.0,), 0.5),
        DistSpec(:Normal, (1.0, 2.0), 0.5),
        #DistSpec(:Normal, (), [0.5]),
        #DistSpec(:Normal, (1.0,), [0.5]),
        #DistSpec(:Normal, (1.0, 2.0), [0.5]),

        DistSpec(:NormalCanon, (1.0, 2.0), 0.5),
        #DistSpec(:NormalCanon, (1.0, 2.0), [0.5]),

        DistSpec(:NormalInverseGaussian, (1.0, 2.0, 1.0, 1.0), 0.5),
        #DistSpec(:NormalInverseGaussian, (1.0, 2.0, 1.0, 1.0), [0.5]),

        DistSpec(:Pareto, (), 1.5),
        DistSpec(:Pareto, (1.0,), 1.5),
        DistSpec(:Pareto, (1.0, 1.0), 1.5),
        #DistSpec(:Pareto, (), [1.5]),
        #DistSpec(:Pareto, (1.0,), [1.5]),
        #DistSpec(:Pareto, (1.0, 1.0), [1.5]),

        DistSpec(:PGeneralizedGaussian, (), 0.5),
        DistSpec(:PGeneralizedGaussian, (1.0, 1.0, 1.0), 0.5),
        #DistSpec(:PGeneralizedGaussian, (), [0.5]),
        #DistSpec(:PGeneralizedGaussian, (1.0, 1.0, 1.0), [0.5]),

        DistSpec(:Rayleigh, (), 0.5),
        DistSpec(:Rayleigh, (1.0,), 0.5),
        #DistSpec(:Rayleigh, (), [0.5]),
        #DistSpec(:Rayleigh, (1.0,), [0.5]),

        DistSpec(:Semicircle, (1.0,), 0.5),
        #DistSpec(:Semicircle, (1.0,), [0.5]),

        DistSpec(:SymTriangularDist, (), 0.5),
        DistSpec(:SymTriangularDist, (1.0,), 0.5),
        DistSpec(:SymTriangularDist, (1.0, 2.0), 0.5),
        #DistSpec(:SymTriangularDist, (), [0.5]),
        #DistSpec(:SymTriangularDist, (1.0,), [0.5]),
        #DistSpec(:SymTriangularDist, (1.0, 2.0), [0.5]),

        DistSpec(:TDist, (1.0,), 0.5),
        #DistSpec(:TDist, (1.0,), [0.5]),

        DistSpec(:TriangularDist, (1.0, 2.0), 1.5),
        DistSpec(:TriangularDist, (1.0, 3.0, 2.0), 1.5),
        #DistSpec(:TriangularDist, (1.0, 2.0), [1.5]),
        #DistSpec(:TriangularDist, (1.0, 3.0, 2.0), [1.5]),

        DistSpec(:Triweight, (1.0, 1.0), 1.0),
        #DistSpec(:Triweight, (1.0, 1.0), [1.0]),

        DistSpec(:((mu, sigma, l, u) -> truncated(Normal(mu, sigma), l, u)), (0.0, 1.0, 1.0, 2.0), 1.5),
        #DistSpec(:((mu, sigma, l, u) -> truncated(Normal(mu, sigma), l, u)), (0.0, 1.0, 1.0, 2.0), [1.5]),

        DistSpec(:Uniform, (), 0.5),
        DistSpec(:Uniform, (0.0, 1.0), 0.5),
        #DistSpec(:Uniform, (), [0.5]),
        #DistSpec(:Uniform, (0.0, 1.0), [0.5]),

        DistSpec(:TuringUniform, (), 0.5),
        DistSpec(:TuringUniform, (0.0, 1.0), 0.5),
        #DistSpec(:TuringUniform, (), [0.5]),
        #DistSpec(:TuringUniform, (0.0, 1.0), [0.5]),

        DistSpec(:VonMises, (), 1.0),
        #DistSpec(:VonMises, (), [1.0]),

        DistSpec(:Weibull, (), 1.0),
        DistSpec(:Weibull, (1.0,), 1.0),
        DistSpec(:Weibull, (1.0, 1.0), 1.0),
        #DistSpec(:Weibull, (), [1.0]),
        #DistSpec(:Weibull, (1.0,), [1.0]),
        #DistSpec(:Weibull, (1.0, 1.0), [1.0]),
    ]
    broken_uni_cont_dists = [
        # Zygote
        DistSpec(:Chernoff, (), 0.5),
        # Broken in Distributions even without autodiff
        DistSpec(:(()->KSDist(1)), (), 0.5),
        DistSpec(:(()->KSOneSided(1)), (), 0.5),
        DistSpec(:StudentizedRange, (1.0, 2.0), 0.5),
        # Dispatch error caused by ccall
        DistSpec(:NoncentralBeta, (1.0, 2.0, 1.0), 0.5),
        DistSpec(:NoncentralChisq, (1.0, 2.0), 0.5),
        DistSpec(:NoncentralF, (1, 2, 1), 0.5),
        DistSpec(:NoncentralT, (1, 2), 0.5),
        # Stackoverflow caused by SpecialFunctions.besselix
        DistSpec(:VonMises, (1.0,), 1.0),
        DistSpec(:VonMises, (1, 1), 1),
    ]

    multi_disc_dists = [
        # Vector x
        DistSpec(:((p) -> Multinomial(2, p / sum(p))), (fill(0.5, 2),), [2, 0]),
        # Matrix x
        #DistSpec(:((p) -> Multinomial(2, p / sum(p))), (fill(0.5, 2),), [2 2; 0 0]);
    ]
    xmulti_disc_dists = [
        multi_disc_dists;
        # Vector x
        filter(!isnothing, filldist_spec.(uni_disc_dists; n = 2, d = 1));
        filter(!isnothing, arraydist_spec.(uni_disc_dists; n = 2, d = 1));
        # Matrix x
        filter(!isnothing, filldist_spec.(uni_disc_dists; n = 2, d = 2));
        filter(!isnothing, arraydist_spec.(uni_disc_dists; n = 2, d = 2));
    ]

    multi_cont_dists = [
        # Vector x
        DistSpec(:MvNormal, (dmean, cov_mat), norm_val_vec),
        DistSpec(:MvNormal, (dmean, cov_vec), norm_val_vec),
        DistSpec(:MvNormal, (dmean, Diagonal(cov_vec)), norm_val_vec),
        DistSpec(:MvNormal, (dmean, cov_num), norm_val_vec),
        DistSpec(:((m, v) -> MvNormal(m, v*I)), (dmean, cov_num), norm_val_vec),
        DistSpec(:MvNormal, (cov_mat,), norm_val_vec),
        DistSpec(:MvNormal, (cov_vec,), norm_val_vec),
        DistSpec(:MvNormal, (Diagonal(cov_vec),), norm_val_vec),
        DistSpec(:(cov_num -> MvNormal($ddim, cov_num)), (cov_num,), norm_val_vec),
        DistSpec(:TuringMvNormal, (dmean, cov_mat), norm_val_vec),
        DistSpec(:TuringMvNormal, (dmean, cov_vec), norm_val_vec),
        DistSpec(:TuringMvNormal, (dmean, Diagonal(cov_vec)), norm_val_vec),
        DistSpec(:TuringMvNormal, (dmean, cov_num), norm_val_vec),
        DistSpec(:((m, v) -> TuringMvNormal(m, v*I)), (dmean, cov_num), norm_val_vec),
        DistSpec(:TuringMvNormal, (cov_mat,), norm_val_vec),
        DistSpec(:TuringMvNormal, (cov_vec,), norm_val_vec),
        DistSpec(:TuringMvNormal, (Diagonal(cov_vec),), norm_val_vec),
        DistSpec(:(cov_num -> TuringMvNormal($ddim, cov_num)), (cov_num,), norm_val_vec),
        DistSpec(:MvLogNormal, (dmean, cov_mat), norm_val_vec),
        DistSpec(:MvLogNormal, (dmean, cov_vec), norm_val_vec),
        DistSpec(:MvLogNormal, (dmean, Diagonal(cov_vec)), norm_val_vec),
        DistSpec(:MvLogNormal, (dmean, cov_num), norm_val_vec),
        DistSpec(:MvLogNormal, (cov_mat,), norm_val_vec),
        DistSpec(:MvLogNormal, (cov_vec,), norm_val_vec),
        DistSpec(:MvLogNormal, (Diagonal(cov_vec),), norm_val_vec),
        DistSpec(:(cov_num -> MvLogNormal($ddim, cov_num)), (cov_num,), norm_val_vec),

        DistSpec(:Dirichlet, (alpha,), dir_val_vec),

        # Matrix case
        DistSpec(:MvNormal, (dmean, cov_vec), norm_val_mat),
        DistSpec(:MvNormal, (dmean, Diagonal(cov_vec)), norm_val_mat),
        DistSpec(:MvNormal, (dmean, cov_num), norm_val_mat),
        DistSpec(:((m, v) -> MvNormal(m, v*I)), (dmean, cov_num), norm_val_mat),
        DistSpec(:MvNormal, (cov_vec,), norm_val_mat),
        DistSpec(:MvNormal, (Diagonal(cov_vec),), norm_val_mat),
        DistSpec(:(cov_num -> MvNormal($ddim, cov_num)), (cov_num,), norm_val_mat),
        DistSpec(:MvLogNormal, (dmean, cov_vec), norm_val_mat),
        DistSpec(:MvLogNormal, (dmean, Diagonal(cov_vec)), norm_val_mat),
        DistSpec(:MvLogNormal, (dmean, cov_num), norm_val_mat),
        DistSpec(:MvLogNormal, (cov_vec,), norm_val_mat),
        DistSpec(:MvLogNormal, (Diagonal(cov_vec),), norm_val_mat),
        DistSpec(:(cov_num -> MvLogNormal($ddim, cov_num)), (cov_num,), norm_val_mat),

        DistSpec(:Dirichlet, (alpha,), dir_val_mat),
    ]
    xmulti_cont_dists = [
        # Vector x
        filter(!isnothing, filldist_spec.(uni_cont_dists; n = 2, d = 1));
        filter(!isnothing, arraydist_spec.(uni_cont_dists; n = 2, d = 1));

        # Matrix case
        # Doesn't work for some reason now but not too important
        #filter(!isnothing, filldist_spec.(uni_cont_dists; n = 2, d = 2));
        #filter(!isnothing, arraydist_spec.(uni_cont_dists; n = 2, d = 2));

        multi_cont_dists;
    ]
    broken_multi_cont_dists = [
        # Dispatch error
        DistSpec(:MvNormalCanon, (dmean, cov_mat), norm_val_vec),
        DistSpec(:MvNormalCanon, (dmean, cov_vec), norm_val_vec),
        DistSpec(:MvNormalCanon, (dmean, cov_num), norm_val_vec),
        DistSpec(:MvNormalCanon, (cov_mat,), norm_val_vec),
        DistSpec(:MvNormalCanon, (cov_vec,), norm_val_vec),
        DistSpec(:(cov_num -> MvNormalCanon($ddim, cov_num)), (cov_num,), norm_val_vec),
        DistSpec(:MvNormalCanon, (dmean, cov_mat), norm_val_mat),
        DistSpec(:MvNormalCanon, (dmean, cov_vec), norm_val_mat),
        DistSpec(:MvNormalCanon, (dmean, cov_num), norm_val_mat),
        DistSpec(:MvNormalCanon, (cov_mat,), norm_val_mat),
        DistSpec(:MvNormalCanon, (cov_vec,), norm_val_mat),
        DistSpec(:(cov_num -> MvNormalCanon($ddim, cov_num)), (cov_num,), norm_val_mat),
        # Test failure
        DistSpec(:MvNormal, (dmean, cov_mat), norm_val_mat),
        DistSpec(:MvNormal, (cov_mat,), norm_val_mat),
        DistSpec(:MvLogNormal, (dmean, cov_mat), norm_val_mat),
        DistSpec(:MvLogNormal, (cov_mat,), norm_val_mat),
    ]

    matrix_cont_dists = [
        # Matrix x
        DistSpec(:((n1, n2)->MatrixBeta($ddim, n1, n2)), (3.0, 3.0), beta_mat),
        DistSpec(:Wishart, (3.0, cov_mat), cov_mat),
        DistSpec(:InverseWishart, (3.0, cov_mat), cov_mat),
        DistSpec(:TuringWishart, (3.0, cov_mat), cov_mat),
        DistSpec(:TuringInverseWishart, (3.0, cov_mat), cov_mat),
        # Vector of matrices x
        DistSpec(:((n1, n2)->MatrixBeta($ddim, n1, n2)), (3.0, 3.0), fill(beta_mat, 2)),
        DistSpec(:Wishart, (3.0, cov_mat), fill(cov_mat, 2)),
        DistSpec(:InverseWishart, (3.0, cov_mat), fill(cov_mat, 2)),
        DistSpec(:TuringWishart, (3.0, cov_mat), fill(cov_mat, 2)),
        DistSpec(:TuringInverseWishart, (3.0, cov_mat), fill(cov_mat, 2)),
    ]
    xmatrix_cont_dists = [
        # Matrix x
        filter(!isnothing, filldist_spec.(uni_cont_dists; n = (2, 2)));
        filter(!isnothing, filldist_spec.(multi_cont_dists; disttype = :multi, n = 2));
        filter(!isnothing, arraydist_spec.(uni_cont_dists; n = (2, 2)));
        filter(!isnothing, arraydist_spec.(multi_cont_dists; disttype = :multi, n = 2));

        # Vector of matrices x
        filter(!isnothing, filldist_spec.(uni_cont_dists; n = (2, 2), d = 2));
        filter(!isnothing, filldist_spec.(multi_cont_dists; disttype = :multi, n = 2, d = 2));
        filter(!isnothing, arraydist_spec.(uni_cont_dists; n = (2, 2), d = 2));
        filter(!isnothing, arraydist_spec.(multi_cont_dists; disttype = :multi, n = 2, d = 2));

        matrix_cont_dists;
    ]
    broken_matrix_cont_dists = [
        # Other
        DistSpec(:MatrixNormal, (cov_mat, cov_mat, cov_mat), cov_mat),
        DistSpec(:(()->MatrixNormal($ddim, $ddim)), (), cov_mat),
        DistSpec(:MatrixTDist, (1.0, cov_mat, cov_mat, cov_mat), cov_mat),
        DistSpec(:MatrixFDist, (3.0, 3.0, cov_mat), cov_mat),
    ]

    test_head(s) = println("\n"*s*"\n")
    separator() = println("\n"*"="^50)

    separator()
    @testset "Univariate discrete distributions" begin
        test_head("Testing: Univariate discrete distributions")
        for d in uni_disc_dists
            test_info(d.name)
            for testf in get_all_functions(d, false)
                test_ad(testf.f, testf.x)
            end
        end
    end
    separator()

    # Note: broadcasting logpdf with univariate distributions having tracked parameters breaks
    # Tracker. Ref: https://github.com/FluxML/Tracker.jl/issues/65
    # filldist works around this so it is the recommended way for AD-friendly "broadcasting"
    # of logpdf with univariate distributions

    @testset "Univariate continuous distributions" begin
        test_head("Testing: Univariate continuous distributions")
        for d in uni_cont_dists
            test_info(d.name)
            for testf in get_all_functions(d, true)
                test_ad(testf.f, testf.x)
            end
        end
    end
    separator()

    @testset "Multivariate discrete distributions" begin
        test_head("Testing: Multivariate discrete distributions")
        for d in xmulti_disc_dists
            test_info(d.name)
            for testf in get_all_functions(d, false)
                test_ad(testf.f, testf.x)
            end
        end
    end

    separator()
    @testset "Multivariate continuous distributions" begin
        test_head("Testing: Multivariate continuous distributions")
        for d in xmulti_cont_dists
            test_info(d.name)
            for testf in get_all_functions(d, true)
                test_ad(testf.f, testf.x)
            end
        end
    end
    separator()

    @testset "Matrix-variate continuous distributions" begin
        test_head("Testing: Matrix-variate continuous distributions")
        for d in xmatrix_cont_dists
            test_info(d.name)
            for testf in get_all_functions(d, true)
                test_ad(testf.f, testf.x)
            end
        end
    end
    separator()
end


text:
export regularize
export L1Regularize
export L2Regularize
export sparsityof
export sparsitypair

"""
    regularize(Ws::Vector; lambda::AbstractFloat, method::String="L2")

apply L1 or L2 regularization to unbiased weights, recurrent weights not included
"""
function regularize(Ws::Vector; lambda::AbstractFloat, method::String="L2")
    method = uppercase(method)
    method=="L2" && return L2Regularize(Ws, lambda)
    method=="L1" && return L1Regularize(Ws, lambda)
end

"""
    L1Regularize(Ws::Vector, λ::AbstractFloat)

apply L1 regularization to unbiased weights, recurrent weights not included
"""
function L1Regularize(Ws::Vector, λ::AbstractFloat)
    Threads.@threads for w in Ws
        @. w += -λ * sign(w)
    end
end

"""
    L2Regularize(Ws::Vector, λ::AbstractFloat)

apply L2 regularizationto unbiased weights, recurrent weights not included
"""
function L2Regularize(Ws::Vector, λ::AbstractFloat)
    Threads.@threads for w in Ws
        @. w += -λ * w
    end
end

"""
    sparsityof(x, a::AbstractFloat) -> sparserRatio

calculate the sparsity of Array `x` at given threshold `a > 0`
"""
function sparsityof(x, a::AbstractFloat)
    @assert a>0.0 "the threshold should > 0, but got $a"
    return sum(abs.(x) .< a) / length(x)
end

"""
    sparsitypair(x, a::AbstractFloat) -> (sparsity::Int, length(x)::Int)

calculate the #sparsity of Array `x` at given threshold `a > 0` and return (sparsity, length(x))
"""
function sparsitypair(x, a::AbstractFloat)
    @assert a>0.0 "the threshold should > 0, but got $a"
    return sum(abs.(x) .< a), length(x)
end

"""
    sparsityof(xs::Vector, a::AbstractFloat) -> sparserRatio

calculate the sparsity of a Vector of Arrays `xs` at given threshold `a > 0`
"""
function sparsityof(xs::Vector, a::AbstractFloat)
    @assert a>0.0 "the threshold should > 0, but got $a"
    S = 0
    T = 0
    for x in xs
        s, t = sparsitypair(x, a)
        S += s
        T += t
    end
    return S/T
end

ext: jl

max_stars_repo_name: sonosole/Delta

max_issues_repo_name: sonosole/Delta

max_forks_repo_name: sonosole/Delta


text:
"""
Type for Basic ABC Output
# Arguments
* `particles`
* `distances`
* `acceptbw`
* `nsampled`
"""
type ABCResult{T <: Number, G <: Real}
    particles::Union{Array{T, 1}, Array{T, 2}}
    distances::Union{Array{G, 1}, Array{G, 2}}
    acceptbw::Union{G, Array{G, 1}}
    nsampled::Int64
end

"""
Standard ABC algorithm for univariate parameter of interest
"""
function abc_standard1D(T::Type, summaryStatistics::Any,
                                   N::Int64,
                                   acceptbw::G,
                                   sample_prior::Function,
                                   forward_model::Function,
                                   compute_distance::Function,
                                   accept_reject::Function) where G <: Real
    #Allocate variables
    acceptedDraws = Array{T}(N)
    acceptedDistances = Array{G}(N)
    accepted = 0
    nsampled = 0

    #Repeat sampling until N particles accepted
    while accepted < N
        nsampled += 1

        #Draw sample and simulate data
        proposal = sample_prior()
        simulatedData = forward_model(proposal)
        proposalDistance = compute_distance(simulatedData, summaryStatistics)

        #Retain sample is distance is below acceptbw
        #if proposalDistance < acceptbw
        if accept_reject(proposalDistance, acceptbw)
            accepted += 1
            acceptedDraws[accepted] = proposal
            acceptedDistances[accepted] = proposalDistance
        end
    end

    #Return result
    return ABCResult(acceptedDraws, acceptedDistances, acceptbw, nsampled)
end

function abc_standard1D(T::Type, summaryStatistics::Any,
                                   N::Int64,
                                   acceptbw::Array{G, 1},
                                   sample_prior::Function,
                                   forward_model::Function,
                                   compute_distance::Function,
                                   accept_reject::Function) where G <: Real

    #Allocate variables
    acceptedDraws = Array{T}(N)
    acceptedDistances = Array{G}(N, length(acceptbw))
    accepted = 0
    nsampled = 0

    #Repeat sampling until N particles accepted
    while accepted < N
        nsampled += 1

        #Draw sample and simulate data
        proposal = sample_prior()
        simulatedData = forward_model(proposal)
        proposalDistance = compute_distance(simulatedData, summaryStatistics)

        #Retain sample is distance is below acceptbw
        #if all(proposalDistance .< acceptbw)
        if accept_reject(proposalDistance, acceptbw)
            accepted += 1
            acceptedDraws[accepted] = proposal
            acceptedDistances[accepted, :] = proposalDistance
        end
    end

    #Return result
    return ABCResult(acceptedDraws, acceptedDistances, acceptbw, nsampled)
end

"""
Standard ABC algorithm for multivariate parameter of interest
"""
function abc_standardMultiD(T::Type, summaryStatistics::Any,
                                       N::Int64,
                                       d::Int64,
                                       acceptbw::G,
                                       sample_prior::Function,
                                       forward_model::Function,
                                       compute_distance::Function,
                                       accept_reject::Function) where G <: Real

    #Allocate variables
    acceptedDraws = Array{T}(N, d)
    acceptedDistances = Array{G}(N)
    accepted = 0
    nsampled = 0

    #Repeat sampling until N particles accepted
    while accepted < N
        nsampled += 1

        #Draw sample and simulate data
        proposal = sample_prior()
        simulatedData = forward_model(proposal)
        proposalDistance = compute_distance(simulatedData, summaryStatistics)

        #Retain sample is distance is below acceptbw
        #if proposalDistance < acceptbw
        if accept_reject(proposalDistance, acceptbw)
            accepted += 1
            acceptedDraws[accepted, :] = proposal
            acceptedDistances[accepted] = proposalDistance
        end
    end

    #Return result
    return ABCResult(acceptedDraws, acceptedDistances, acceptbw, nsampled)
end

function abc_standardMultiD(T::Type, summaryStatistics::Any,
                                       N::Int64,
                                       d::Int64,
                                       acceptbw::Array{G, 1},
                                       sample_prior::Function,
                                       forward_model::Function,
                                       compute_distance::Function,
                                       accept_reject::Function) where G <: Real
    acceptedDraws = Array{T}(N, d)
    acceptedDistances = Array{G}(N, length(acceptbw))
    accepted = 0
    nsampled = 0

    #Repeat sampling until N particles accepted
    while accepted < N
        nsampled += 1

        #Draw sample and simulate data
        proposal = sample_prior()
        simulatedData = forward_model(proposal)
        proposalDistance = compute_distance(simulatedData, summaryStatistics)

        #Retain sample is distance is below acceptbw
        #if all(proposalDistance .< acceptbw)
        if accept_reject(proposalDistance, acceptbw)
            accepted += 1
            acceptedDraws[accepted, :] = proposal
            acceptedDistances[accepted, :] = proposalDistance
        end
    end

    #Return result
    return ABCResult(acceptedDraws, acceptedDistances, acceptbw, nsampled)
end

"""
Standard ABC Algorithm, wrapper that determines dimention of algorithm to be called and runs appropriate version

# Arguments
* `summaryStatistics::Any` summary statistics of observed data used in
distance comparison.
* `N::Int64` Number of samples to accept before termination
* `acceptbw::Float64 or ::Array{Float64, 1}` acceptance acceptbw
* `sample_prior::Function` Function to draw a sample from the prior distribution from either a univariate or multivariate distribution called by `sample_prior()`.
* `forward_model::Function` Function to simulate data based on draw from parameter prior.  The call `forward_model(sample_prior())` should run.
* `compute_distance::Function` Function to compute distance between output of the `forward_model` and the provided `summaryStatistics`. The call `compute_distance(forward_model(sample_prior()), summaryStatistics)` should run.

# Value
Returns an object of type `ABCResult` containing the accepted particles, corresponding distances, acceptance acceptbw, and the number of samples drawn from the prior before the desired number of particles were accepted

# Examples
```julia
using ApproximateBayesianComputation

data = rand(Distributions.Normal(2.0, 0.5), 100)

#Construct Appropriate Functions
sample_prior = make_sample_prior(Distributions.Normal(0.0, 0.5))
function forward_model(μ::Float64)
    return rand(Distributions.Normal(μ, 0.5), 100)
end
function compute_distance(x::Array{Float64, 1}, y::Float64)
    return abs(mean(x) - y)
end

result = abc_standard(mean(data), 200, 0.05, sample_prior, forward_model, compute_distance, boxcar_accept)
```
"""
function abc_standard(summaryStatistics::Any,
                                 N::Int64,
                                 acceptbw::Union{G, Array{G, 1}},
                                 sample_prior::Function,
                                 forward_model::Function,
                                 compute_distance::Function;
                                 accept_reject::Function = boxcar_accept) where G <: Real
    proposal = sample_prior()
    if typeof(proposal) <: Array
        return abc_standardMultiD(eltype(proposal), summaryStatistics, N, length(proposal), acceptbw, sample_prior, forward_model, compute_distance, accept_reject)
    else
        return abc_standard1D(typeof(proposal), summaryStatistics, N, acceptbw, sample_prior, forward_model, compute_distance, accept_reject)
    end
end


text:
module MockTests

using Fixtures
using FactCheck

facts("Mock tests") do
  context("Can make, call and inspect and reset a mock") do
    my_mock = mock()

    @fact calls(my_mock) --> []

    result = my_mock()

    @fact result --> nothing

    @fact calls(my_mock) --> [call()]

    reset(my_mock)

    @fact calls(my_mock) --> []
  end

  context("Multiple calls") do
    my_mock = mock()

    @fact calls(my_mock) --> []

    result = my_mock()

    @fact result --> nothing

    @fact calls(my_mock) --> [call()]

    result = my_mock()

    @fact result --> nothing

    @fact calls(my_mock) --> [call(), call()]

    result = my_mock()

    @fact result --> nothing

    @fact calls(my_mock) --> [call(), call(), call()]

    reset(my_mock)

    @fact calls(my_mock) --> []
  end

  context("Can specify a return value") do
    my_mock = mock(return_value=12)
    @fact my_mock() --> 12
  end

  context("Can specify a side effect") do
    my_mock = mock(side_effect=() -> 123456)
    @fact my_mock() --> 123456
  end

  context("If you specify a side effect and a return value the side effect wins") do
    my_mock = mock(return_value=12, side_effect=() -> 123456)
    @fact my_mock() --> 123456
  end

  context("Passing an argument to a mock") do
    my_mock = mock()

    @fact calls(my_mock) --> []

    result = my_mock(100)

    @fact result --> nothing

    @fact calls(my_mock) --> [call(100)]

    result = my_mock(200)

    @fact result --> nothing

    @fact calls(my_mock) --> [call(100), call(200)]

    reset(my_mock)

    @fact calls(my_mock) --> []
  end

  context("Passing a named argument to a mock") do
    my_mock = mock()

    @fact calls(my_mock) --> []

    result = my_mock(thingy=100)

    @fact result --> nothing

    @fact calls(my_mock) --> [call(thingy=100)]

    result = my_mock(whatnot=200)

    @fact result --> nothing

    @fact calls(my_mock) --> [call(thingy=100), call(whatnot=200)]

    reset(my_mock)

    @fact calls(my_mock) --> []
  end

  context("Passing named and unnamed arguments to a mock") do
    my_mock = mock()

    @fact calls(my_mock) --> []

    result = my_mock(1,2,3,thingy=100)

    @fact result --> nothing

    @fact calls(my_mock) --> [call(1,2,3,thingy=100)]

    result = my_mock(4,6, whatnot=200, doodah=300)

    @fact result --> nothing

    @fact calls(my_mock) --> [call(1,2,3,thingy=100), call(4,6, whatnot=200, doodah=300)]

    reset(my_mock)

    @fact calls(my_mock) --> []
  end

  context("Order of named arguments matters") do
    @fact call(foo=1, bar=2) --> not(call(bar=2, foo=1))

    # But the position of the unamed arguments reletive to the named ones does not
    @fact call(0, foo=1, bar=2) --> call(foo=1, 0, bar=2)
    @fact call(0, foo=1, bar=2) --> call(foo=1, bar=2, 0)
  end

  context("Calling one mock doesn't effect the other(s)") do
    mock1 = mock()
    mock2 = mock()
    mock3 = mock()

    @fact calls(mock1) --> []
    @fact calls(mock2) --> []
    @fact calls(mock3) --> []

    @fact mock1(100) --> nothing

    @fact calls(mock1) --> [call(100)]
    @fact calls(mock2) --> []
    @fact calls(mock3) --> []

    @fact mock1(200) --> nothing

    @fact calls(mock1) --> [call(100), call(200)]
    @fact calls(mock2) --> []
    @fact calls(mock3) --> []
  end

  context("Resetting one mock doesn't effect the other(s)") do
    mock1 = mock()
    mock2 = mock()
    mock3 = mock()

    @fact calls(mock1) --> []
    @fact calls(mock2) --> []
    @fact calls(mock3) --> []

    @fact mock1(100) --> nothing
    @fact mock2(200) --> nothing
    @fact mock3(300) --> nothing

    @fact mock1(400) --> nothing
    @fact mock2(500) --> nothing
    @fact mock3(600) --> nothing

    @fact calls(mock1) --> [call(100), call(400)]
    @fact calls(mock2) --> [call(200), call(500)]
    @fact calls(mock3) --> [call(300), call(600)]

    reset(mock1)

    @fact calls(mock1) --> []
    @fact calls(mock2) --> [call(200), call(500)]
    @fact calls(mock3) --> [call(300), call(600)]
  end

  #context("Using ANYTHING with a mock") do
  #  my_mock = mock()
  #  @fact my_mock(100) --> nothing
  #  @fact calls(my_mock) --> [call(ANYTHING)]
  #end

  context("Using a mock to patch Base.open") do
    function firstline(filename)
      f = open(filename)
      try
        return chomp(readlines(f)[1])
      finally
        close(f)
      end
    end

    mock_open = mock(return_value=IOBuffer("Hello Julia\nGoodbye Julia"))

    patch(Base, :open, mock_open) do
      @fact firstline("foobar.txt") --> "Hello Julia"
    end

    @fact calls(mock_open) --> [call("foobar.txt")]
  end
end

end


text:
""" b5 (spinel) lattice """
b5(T::Type=Float64; unit=u"nm") = Crystal(
    T[0 0.5 0.5; 0.5 0 0.5; 0.5 0.5 0] * unit,
    tposition=T[0.5 0.5 0.5;
                0.5 0.25 0.25;
                0.25 0.5 0.25;
                0.25 0.25 0.5;
                0.875 0.875 0.875;
                0.125 0.125 0.125;
                0.25 0.25 0.25;
                0.25 0.5 0.5;
                0.5 0.25 0.5;
                0.5 0.5 0.25;
                0.75 0.75 0.75;
                0.75 0.5 0.5;
                0.5 0.75 0.5;
                0.5 0.5 0.75] * unit,
    species=['A', 'A', 'A', 'A', 'B', 'B', 'X', 'X', 'X', 'X', 'X', 'X', 'X', 'X']
)


text:
using Dates
using Logging
using Printf
using LinearAlgebra

using MPI
using StaticArrays
using NCDatasets

using ClimateMachine
ClimateMachine.init(parse_clargs = true)

using ClimateMachine.BalanceLaws: vars_state_conservative, vars_state_auxiliary
using ClimateMachine.Mesh.Topologies
using ClimateMachine.Mesh.Grids
using ClimateMachine.Mesh.Filters
using ClimateMachine.DGMethods
using ClimateMachine.DGMethods.NumericalFluxes
using ClimateMachine.MPIStateArrays
using ClimateMachine.ODESolvers
using ClimateMachine.VariableTemplates: flattenednames
using ClimateMachine.SplitExplicit01
using ClimateMachine.GenericCallbacks
using ClimateMachine.VTK

using CLIMAParameters
using CLIMAParameters.Planet: grav

struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

import ClimateMachine.SplitExplicit01:
    ocean_init_aux!,
    ocean_init_state!,
    ocean_boundary_state!,
    CoastlineFreeSlip,
    CoastlineNoSlip,
    OceanFloorFreeSlip,
    OceanFloorNoSlip,
    OceanSurfaceNoStressNoForcing,
    OceanSurfaceStressNoForcing,
    OceanSurfaceNoStressForcing,
    OceanSurfaceStressForcing

import ClimateMachine.DGMethods:
    update_auxiliary_state!,
    update_auxiliary_state_gradient!,
    vars_state_conservative,
    vars_state_auxiliary,
    VerticalDirection

const ArrayType = ClimateMachine.array_type()

struct SimpleBox{T} <: AbstractOceanProblem
    Lˣ::T
    Lʸ::T
    H ::T
    τₒ::T
    λʳ::T
    θᴱ::T
end

@inline function ocean_boundary_state!(m::OceanModel, p::SimpleBox, bctype, x...)
    if bctype == 1
        ocean_boundary_state!(m, CoastlineNoSlip(), x...)
    elseif bctype == 2
        ocean_boundary_state!(m, OceanFloorNoSlip(), x...)
    elseif bctype == 3
        ocean_boundary_state!(m, OceanSurfaceStressForcing(), x...)
    end
end

@inline function ocean_boundary_state!(m::Continuity3dModel, p::SimpleBox, bctype, x...)
   #if bctype == 1
        ocean_boundary_state!(m, CoastlineNoSlip(), x...)
   #end
end

@inline function ocean_boundary_state!(m::BarotropicModel, p::SimpleBox, bctype, x...)
    return ocean_boundary_state!(m, CoastlineNoSlip(), x...)
end

function ocean_init_state!(p::SimpleBox, Q, A, coords, t)
    @inbounds y = coords[2]
    @inbounds z = coords[3]
    @inbounds H = p.H

    Q.u = @SVector [-0, -0]
    Q.η = -0
    Q.θ = (5 + 4 * cos(y * π / p.Lʸ)) * (1 + z / H)

    return nothing
end

function ocean_init_aux!(m::OceanModel, p::SimpleBox, A, geom)
    FT = eltype(A)
    @inbounds A.y = geom.coord[2]

    # not sure if this is needed but getting weird intialization stuff
    A.w = -0
    A.pkin = -0
    A.wz0 = -0
    A.u_d = @SVector [-0, -0]
    A.ΔGu = @SVector [-0, -0]

    return nothing
end

# A is Filled afer the state
function ocean_init_aux!(m::BarotropicModel, P::SimpleBox, A, geom)
    @inbounds A.y = geom.coord[2]

    A.Gᵁ = @SVector [-0, -0]
    A.U_c = @SVector [-0, -0]
    A.η_c = -0
    A.U_s = @SVector [-0, -0]
    A.η_s = -0
    A.Δu = @SVector [-0, -0]
    A.η_diag = -0
    A.Δη = -0

    return nothing
end

function make_callbacks(step, nout, mpicomm, odesolver, dg_slow, model_slow, Q_slow, dg_fast, model_fast, Q_fast)
    Np = N
    Nᵖ⁺¹ = Np + 1

    ΣNˣ = (Np+1)*Nˣ
    ΣNʸ = (Np+1)*Nʸ
    ΣNᶻ = (Np+1)*Nᶻ

    gnd = reshape(dg_slow.grid.vgeo, (Np+1, Np+1, Np+1, 16, Nᶻ, Nʸ, Nˣ))
    x = gnd[:, :, :, 13, :, :, :] |> Array
    y = gnd[:, :, :, 14, :, :, :] |> Array
    z = gnd[:, :, :, 15, :, :, :] |> Array

    ΔX = (maximum(x) - minimum(x)) / Nˣ
    ΔY = (maximum(y) - minimum(y)) / Nʸ
    ΔZ = (maximum(z) - minimum(z)) / Nᶻ

    ds = NCDataset("simple_box_ivd.nc", "c")

    defDim(ds, "time", Inf)
    defVar(ds, "time", Float64, ("time",))

    defDim(ds, "x", ΣNˣ)
    defDim(ds, "y", ΣNʸ)
    defDim(ds, "z", ΣNᶻ)

    defVar(ds, "x_nodal", Float64, ("x", "y", "z"))
    defVar(ds, "y_nodal", Float64, ("x", "y", "z"))
    defVar(ds, "z_nodal", Float64, ("x", "y", "z"))

    defVar(ds, "u", Float64, ("x", "y", "z", "time"))
    defVar(ds, "v", Float64, ("x", "y", "z", "time"))
    defVar(ds, "η", Float64, ("x", "y", "z", "time"))
    defVar(ds, "θ", Float64, ("x", "y", "z", "time"))

    netcdf_output = GenericCallbacks.EveryXSimulationSteps(1) do (init = false)
        time_index = length(ds["time"]) + 1

        if time_index == 1
            xs = zeros(ΣNˣ, ΣNʸ, ΣNᶻ)
            ys = zeros(ΣNˣ, ΣNʸ, ΣNᶻ)
            zs = zeros(ΣNˣ, ΣNʸ, ΣNᶻ)
        end

        Q3nd = reshape(Q_slow.realdata, (Np+1, Np+1, Np+1, 4, Nᶻ, Nʸ, Nˣ))
        u = Q3nd[:, :, :, 1, :, :, :] |> Array
        v = Q3nd[:, :, :, 2, :, :, :] |> Array
        η = Q3nd[:, :, :, 3, :, :, :] |> Array
        θ = Q3nd[:, :, :, 4, :, :, :] |> Array

        us = zeros(ΣNˣ, ΣNʸ, ΣNᶻ)
        vs = zeros(ΣNˣ, ΣNʸ, ΣNᶻ)
        ηs = zeros(ΣNˣ, ΣNʸ, ΣNᶻ)
        θs = zeros(ΣNˣ, ΣNʸ, ΣNᶻ)

        for I in 1:Nˣ, J in 1:Nʸ, K in 1:Nᶻ
            I′ = x[:, :, :, K, J, I] ./ ΔX |> maximum |> round |> Int
            J′ = y[:, :, :, K, J, I] ./ ΔY |> maximum |> round |> Int
            K′ = z[:, :, :, K, J, I] ./ ΔZ |> maximum |> round |> Int
            K′ = Nᶻ + K′

            @debug "($I, $J, $K) -> ($I′, $J′, $(Nᶻ+K′))"

            i_elem = (I′-1) * Nᵖ⁺¹ + 1 : I′ * Nᵖ⁺¹
            j_elem = (J′-1) * Nᵖ⁺¹ + 1 : J′ * Nᵖ⁺¹
            k_elem = (K′-1) * Nᵖ⁺¹ + 1 : K′ * Nᵖ⁺¹

            if time_index == 1
                xs[i_elem, j_elem, k_elem] .= x[:, :, :, K, J, I]
                ys[i_elem, j_elem, k_elem] .= y[:, :, :, K, J, I]
                zs[i_elem, j_elem, k_elem] .= z[:, :, :, K, J, I]
            end

            us[i_elem, j_elem, k_elem] .= u[:, :, :, K, J, I]
            vs[i_elem, j_elem, k_elem] .= v[:, :, :, K, J, I]
            ηs[i_elem, j_elem, k_elem] .= η[:, :, :, K, J, I]
            θs[i_elem, j_elem, k_elem] .= θ[:, :, :, K, J, I]
        end

        if time_index == 1
            ds["x_nodal"][:, :, :] = xs
            ds["y_nodal"][:, :, :] = ys
            ds["z_nodal"][:, :, :] = zs
        end

        ds["u"][:, :, :, time_index] = us
        ds["v"][:, :, :, time_index] = vs
        ds["η"][:, :, :, time_index] = ηs
        ds["θ"][:, :, :, time_index] = θs

        sync(ds)
    end

    starttime = Ref(now())
    cbinfo = GenericCallbacks.EveryXWallTimeSeconds(60, mpicomm) do (s = false)
        if s
            starttime[] = now()
        else
            energy = norm(Q_slow)
            @info @sprintf(
                """Update
                simtime = %8.2f / %8.2f
                runtime = %s
                norm(Q) = %.16e""",
                ODESolvers.gettime(odesolver), timeend,
                Dates.format(
                    convert(Dates.DateTime, Dates.now() - starttime[]),
                    Dates.dateformat"HH:MM:SS",
                ),
                energy)
        end
    end

    return (netcdf_output, cbinfo)
end

#################
# RUN THE TESTS #
#################
FT = Float64

const timeend = 5 * 24 * 3600 # s
const tout = 24 * 3600 # s

const N = 4
const Nˣ = 20
const Nʸ = 20
const Nᶻ = 20
const Lˣ = 4e6  # m
const Lʸ = 4e6  # m
const H = 1000  # m

xrange = range(FT(0); length = Nˣ + 1, stop = Lˣ)
yrange = range(FT(0); length = Nʸ + 1, stop = Lʸ)
zrange = range(FT(-H); length = Nᶻ + 1, stop = 0)

#const cʰ = sqrt(gravity * H)
const cʰ = 1  # typical of ocean internal-wave speed
const cᶻ = 0

#- inverse ratio of additional fast time steps (for weighted average)
#  --> do 1/add more time-steps and average from: 1 - 1/add up to: 1 + 1/add
# e.g., = 1 --> 100% more ; = 2 --> 50% more ; = 3 --> 33% more ...
add_fast_substeps = 2

#- number of Implicit vertical-diffusion sub-time-steps within one model full time-step
# default = 0 : disable implicit vertical diffusion
numImplSteps = 5

#const τₒ = 2e-1  # (Pa = N/m^2)
# since we are using old BC (with factor of 2), take only half:
const τₒ = 1e-1
const λʳ = 10 // 86400 # m/s
const θᴱ = 10    # deg.C

mpicomm = MPI.COMM_WORLD

ll = uppercase(get(ENV, "JULIA_LOG_LEVEL", "INFO"))
loglevel = ll == "DEBUG" ? Logging.Debug :
           ll == "WARN"  ? Logging.Warn  :
           ll == "ERROR" ? Logging.Error : Logging.Info
logger_stream = MPI.Comm_rank(mpicomm) == 0 ? stderr : devnull
global_logger(ConsoleLogger(logger_stream, loglevel))

brickrange_2D = (xrange, yrange)
topl_2D = BrickTopology(mpicomm, brickrange_2D, periodicity = (false, false))
grid_2D = DiscontinuousSpectralElementGrid(topl_2D, FloatType = FT, DeviceArray = ArrayType,
                                           polynomialorder = N)

brickrange_3D = (xrange, yrange, zrange)
topl_3D = StackedBrickTopology(mpicomm, brickrange_3D; periodicity = (false, false, false),
                               boundary = ((1, 1), (1, 1), (2, 3)))
grid_3D = DiscontinuousSpectralElementGrid(topl_3D, FloatType = FT, DeviceArray = ArrayType, polynomialorder = N)

prob = SimpleBox{FT}(Lˣ, Lʸ, H, τₒ, λʳ, θᴱ)
gravity = grav(param_set)

#- set model time-step:
dt_fast = 240
dt_slow = 5400

nout = ceil(Int64, tout / dt_slow)
dt_slow = tout / nout
numImplSteps > 0 ? ivdc_dt = dt_slow / FT(numImplSteps) : ivdc_dt = dt_slow

model = OceanModel{FT}(prob, grav = gravity, cʰ = cʰ, add_fast_substeps = add_fast_substeps,
                       numImplSteps = numImplSteps, ivdc_dt = ivdc_dt, κᶜ = FT(0.1))

barotropicmodel = BarotropicModel(model)

minΔx = min_node_distance(grid_3D, HorizontalDirection())
minΔz = min_node_distance(grid_3D, VerticalDirection())

#- 2 horiz directions
gravity_max_dT = 1 / ( 2 * sqrt(gravity * H) / minΔx )

#- 2 horiz directions + harmonic visc or diffusion: 2^2 factor in CFL:
viscous_max_dT = 1 / ( 2 * model.νʰ / minΔx^2 + model.νᶻ / minΔz^2 )/ 4
diffusive_max_dT = 1 / ( 2 * model.κʰ / minΔx^2 + model.κᶻ / minΔz^2 )/ 4

@info @sprintf(
    """Update
       Gravity Max-dT = %.1f
       Timestep       = %.1f""",
    gravity_max_dT, dt_fast)

@info @sprintf(
    """Update
   Viscous   Max-dT = %.1f
   Diffusive Max-dT = %.1f
   Timestep      = %.1f""",
    viscous_max_dT, diffusive_max_dT, dt_slow)

dg = OceanDGModel(
    model, grid_3D,
#   CentralNumericalFluxFirstOrder(),
    RusanovNumericalFlux(),
    CentralNumericalFluxSecondOrder(),
    CentralNumericalFluxGradient(),
)

barotropic_dg = DGModel(
    barotropicmodel, grid_2D,
#   CentralNumericalFluxFirstOrder(),
    RusanovNumericalFlux(),
    CentralNumericalFluxSecondOrder(),
    CentralNumericalFluxGradient(),
)

Q_3D = init_ode_state(dg, FT(0); init_on_cpu = true)
Q_2D = init_ode_state(barotropic_dg, FT(0); init_on_cpu = true)

lsrk_ocean = LSRK54CarpenterKennedy(dg, Q_3D, dt = dt_slow, t0 = 0)
lsrk_barotropic = LSRK54CarpenterKennedy(barotropic_dg, Q_2D, dt = dt_fast, t0 = 0)

odesolver = SplitExplicitLSRK2nSolver(lsrk_ocean, lsrk_barotropic)

# Set up State Check call back for config state arrays, called every ntFreq time steps
ntFreq=1
cbcs_dg=ClimateMachine.StateCheck.sccreate(
        [(Q_3D,"oce Q_3D",),
         (dg.state_auxiliary,"oce aux",),
    #    (dg.diffstate,"oce diff",),
    #    (lsrk_ocean.dQ,"oce_dQ",),
    #    (dg.modeldata.tendency_dg.state_auxiliary,"tend Int aux",),
    #    (dg.modeldata.conti3d_Q,"conti3d_Q",),
         (Q_2D,"baro Q_2D",),
         (barotropic_dg.state_auxiliary ,"baro aux",)
    ],
        ntFreq; prec=12)

step = [0, 0]
cbvector = make_callbacks(step, nout, mpicomm, odesolver, dg, model, Q_3D,
      barotropic_dg, barotropicmodel, Q_2D)

eng0 = norm(Q_3D)

@info @sprintf """Starting
norm(Q₀) = %.16e
ArrayType = %s""" eng0 ArrayType

# slow fast state tuple
Qvec = (slow = Q_3D, fast = Q_2D)
# solve!(Qvec, odesolver; timeend = timeend, callbacks = cbvector)
cbv=(cbvector...,cbcs_dg)
solve!(Qvec, odesolver; timeend = timeend, callbacks = cbv)

## Enable the code block below to print table for use in reference value code
## reference value code sits in a file named $(@__FILE__)_refvals.jl. It is hand
## edited using code generated by block below when reference values are updated.
regenRefVals = false
if regenRefVals
    ## Print state statistics in format for use as reference values
    println(
        "# SC ========== Test number ", 1,
        " reference values and precision match template. =======",
    )
    println("# SC ========== $(@__FILE__) test reference values ======================================")
    ClimateMachine.StateCheck.scprintref(cbcs_dg)
    println("# SC ====================================================================================")
end

## Check results against reference if present
checkRefVals = true
if checkRefVals
    include("simple_box_ivd_refvals.jl")
    refDat = (refVals[1], refPrecs[1])
    checkPass = ClimateMachine.StateCheck.scdocheck(cbcs_dg, refDat)
    checkPass ? checkRep = "Pass" : checkRep = "Fail"
    @info @sprintf("""Compare vs RefVals: %s""", checkRep )
end


text:
using Documenter
using FlightMechanicsUtils

makedocs(
    sitename = "Flight Mechanics Utils Documentation",
    modules = [FlightMechanicsUtils],
    pages = ["Home" => "index.md", "API" => "api.md"],
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
)

deploydocs(
    repo = "github.com:AlexS12/FlightMechanicsUtils.jl.git",
    push_preview=true,
    devbranch = "main"
)


text:
using Crypto
using Base.Test

##############################################################################
##
## Tests for Julia implementations
##
##############################################################################

init()

# TODO: add convert integer to byte array tests, including BigInt

test = bytes2hex(digest("SHA256", ""))
@test test == "e3b0c44298fc1c149afbf4c8996fb92427ae41e4649b934ca495991b7852b855"

test = bytes2hex(digest("RIPEMD160", ""))
@test test == "9c1185a5c5e9fc54612808977ee8f548b2258d31"

test = bytes2hex(digest("SHA256", "abcd", is_hex = true))
@test test == "123d4c7ef2d1600a1b3a0f6addc60a10f05a3495c9409f2ecbf4cc095d000a6b"

@test_throws ErrorException digest("asdf", "abcd")

test = ec_pub_key(hex2bytes("18E14A7B6A307F426A94F8114701E7C8E774E7F9A47E2C2035DB29A206321725"))
@test bytes2hex(test) == "0450863ad64a87ae8a2fe83c1af1a8403cb53f53e486d8511dad8a04887e5b23522cd470243453a299fa9e77237716103abc11a1df38855ed6f2ee187e9c582ba6"

@test length(random(244)) == 31

priv_key = random(256)
pub_key = ec_pub_key(priv_key)

a = "abc".data
ec_sign(a, priv_key)
@test ec_verify(a, ec_sign(a, priv_key), pub_key)

a = "abc"
ec_sign(a, priv_key)
@test ec_verify(a, ec_sign(a, priv_key), pub_key)

# Utilities
# hex2oct: empty string
@test hex2oct("") == []
# hex2oct: odd- and even-length strings
@test hex2oct("adfcef981") == [0x0a, 0xdf, 0xce, 0xf9, 0x81]
@test hex2oct("aadfcef981") == [0xaa, 0xdf, 0xce, 0xf9, 0x81]

# SHA2 tests
# Test empty string
# @test SHA2.sha256("") == "e3b0c44298fc1c149afbf4c8996fb92427ae41e4649b934ca495991b7852b855"

# @test SHA2.sha256("a", is_hex=false) == "ca978112ca1bbdcafac231b39a23dc4da786eff8147c4e72b9807785afee48bb"
# @test SHA2.sha256("Scientific progress goes 'boink'", is_hex=false) == "2f2ba2a09a66771bf1fdf541af6e9db4b443145f9935ddd5d4c323c21a8bdcee"
# @test SHA2.sha256("I'd hold you up to say to your mother, 'this kid's gonna be the best kid in the world. This kid's gonna be somebody better than anybody I ever knew.' And you grew up good and wonderful. It was great just watching you, every day was like a privilege. Then the time come for you to be your own man and take on the world, and you did. But somewhere along the line, you changed. You stopped being you. You let people stick a finger in your face and tell you you're no good. And when things got hard, you started looking for something to blame, like a big shadow. Let me tell you something you already know. The world ain't all sunshine and rainbows. It's a very mean and nasty place and I don't care how tough you are it will beat you to your knees and keep you there permanently if you let it. You, me, or nobody is gonna hit as hard as life. But it ain't about how hard ya hit. It's about how hard you can get hit and keep moving forward. How much you can take and keep moving forward. That's how winning is done! Now if you know what you're worth then go out and get what you're worth. But ya gotta be willing to take the hits, and not pointing fingers saying you ain't where you wanna be because of him, or her, or anybody! Cowards do that and that ain't you! You're better than that! I'm always gonna love you no matter what. No matter what happens. You're my son and you're my blood. You're the best thing in my life. But until you start believing in yourself, ya ain't gonna have a life. Don't forget to visit your mother.", is_hex=false) == "a5d8cfb99203ae8cd0c222e8aaef815a7a53493f650c5dec0d73de7f912e91f2"

# Test > 448 bits (> 56 characters)
# @test SHA2.sha256("asdfghjkqwasdfghjkqwasdfghjkqwasdfghjkqwasdfghjkqwasdfghjkqw", is_hex=false) == "07a95e647687cf0e8cd3d0ca78c9cc9b120ab41497f5f3be912c6c3f1ecd3a31"

# Test hex strings
# @test SHA2.sha256("800c28fca386c7a227600b2fe50b7cae11ec86d3bf1fbe471be89827e19d72aa1d") == "8147786c4d15106333bf278d71dadaf1079ef2d2440a4dde37d747ded5403592"
# @test SHA2.sha256("8147786c4d15106333bf278d71dadaf1079ef2d2440a4dde37d747ded5403592") == "507a5b8dfed0fc6fe8801743720cedec06aa5c6fca72b07c49964492fb98a714"

# RIPEMD tests
# @test RIPEMD.ripemd160("asdf", is_hex=false) == "0ef2aed6346def670a8019e4ea42cf4c76018139"
# @test RIPEMD.ripemd160("", is_hex=false) == "9c1185a5c5e9fc54612808977ee8f548b2258d31"
# @test RIPEMD.ripemd160("a", is_hex=false) == "0bdc9d2d256b3ee9daae347be6f4dc835a467ffe"
# @test RIPEMD.ripemd160("abc", is_hex=false) == "8eb208f7e05d987a9b044a8e98c6b087f15a0bfc"
# @test RIPEMD.ripemd160("message digest", is_hex=false) == "5d0689ef49d2fae572b881b123a85ffa21595f36"
# @test RIPEMD.ripemd160("ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789", is_hex=false) == "b0e20b6e3116640286ed3a87a5713079b21f5189"
# @test RIPEMD.ripemd160("I'd hold you up to say to your mother, 'this kid's gonna be the best kid in the world. This kid's gonna be somebody better than anybody I ever knew.' And you grew up good and wonderful. It was great just watching you, every day was like a privilege. Then the time come for you to be your own man and take on the world, and you did. But somewhere along the line, you changed. You stopped being you. You let people stick a finger in your face and tell you you're no good. And when things got hard, you started looking for something to blame, like a big shadow. Let me tell you something you already know. The world ain't all sunshine and rainbows. It's a very mean and nasty place and I don't care how tough you are it will beat you to your knees and keep you there permanently if you let it. You, me, or nobody is gonna hit as hard as life. But it ain't about how hard ya hit. It's about how hard you can get hit and keep moving forward. How much you can take and keep moving forward. That's how winning is done! Now if you know what you're worth then go out and get what you're worth. But ya gotta be willing to take the hits, and not pointing fingers saying you ain't where you wanna be because of him, or her, or anybody! Cowards do that and that ain't you! You're better than that! I'm always gonna love you no matter what. No matter what happens. You're my son and you're my blood. You're the best thing in my life. But until you start believing in yourself, ya ain't gonna have a life. Don't forget to visit your mother.", is_hex=false) == "fff55c23c197b4fded67e09424e5aef9dafad1c6"

# ECDSA tests
# a = ECDSA.FiniteFields.makeModular(324, 3851)
# b = ECDSA.FiniteFields.makeModular(1287, 3851)
# insecureCurve = ECDSA.EllipticCurves.Curve(a, b)
# x = ECDSA.FiniteFields.makeModular(920, 3851)
# y = ECDSA.FiniteFields.makeModular(303, 3851)
# basePoint = ECDSA.EllipticCurves.ConcretePoint(x, y, insecureCurve)

# aliceSecretKey = rand(0:255)
# bobSecretKey = rand(0:255)

# alicePublicKey = aliceSecretKey * basePoint
# bobPublicKey = bobSecretKey * basePoint

# sharedSecret1 = bobSecretKey * alicePublicKey
# sharedSecret2 = aliceSecretKey * bobPublicKey

# @test sharedSecret1 == sharedSecret2


text:
module IndexUtils
export insertsorted!, deletesorted!, set_attr_index!, unset_attr_index!

""" Insert into sorted vector, preserving the sorting.
"""
function insertsorted!(a::AbstractVector, x)
  insert!(a, searchsortedfirst(a, x), x)
end

""" Delete one occurrence of value from sorted vector, if present.

Returns whether an occurence was found and deleted.
"""
function deletesorted!(a::AbstractVector, x)
  i = searchsortedfirst(a, x)
  found = i <= length(a) && a[i] == x
  if found; deleteat!(a, i) end
  found
end

""" Set key and value for acset attribute index.
"""
function set_attr_index!(d::AbstractDict{K,<:AbstractVector{Int}},
                         k::K, v::Int) where K
  insertsorted!(get!(d, k) do; Int[] end, v)
end

""" Unset key and value from acset attribute index, if it is set.
"""
function unset_attr_index!(d::AbstractDict{K,<:AbstractVector{Int}},
                           k::K, v::Int) where K
  if haskey(d, k)
    vs = d[k]
    if deletesorted!(vs, v) && isempty(vs)
      delete!(d, k)
    end
  end
end

end


text:
using ReinforcementLearning
using Documenter
using Markdown
using DemoCards

open(joinpath(@__DIR__, "..", "README.md"), "r") do f_src
    open(joinpath(@__DIR__, "src", "index.md"), "w") do f_dest
        s_dest = read(f_src, String)
        s_dest = replace(s_dest, "<!-- ```@raw html -->" => "```@raw html")
        s_dest = replace(s_dest, "<!-- ``` -->" => "```")
        write(f_dest, s_dest)
    end
end

experiments, postprocess_cb, experiments_assets = makedemos("experiments")

assets = [
    "assets/favicon.ico",
    "assets/custom.css",
    experiments_assets
]

makedocs(
    modules = [
        ReinforcementLearning,
        ReinforcementLearningBase,
        ReinforcementLearningCore,
        ReinforcementLearningEnvironments,
        ReinforcementLearningZoo,
    ],
    format = Documenter.HTML(
        prettyurls = true,
        analytics = "UA-149861753-1",
        assets = assets,
    ),
    sitename = "ReinforcementLearning.jl",
    linkcheck = !("skiplinks" in ARGS),
    pages = [
        "Home" => "index.md",
        "Tutorial" => "tutorial.md",
        "Guides" => [
            "How to write a customized environment?" => "How_to_write_a_customized_environment.md",
            "How to implement a new algorithm?" => "How_to_implement_a_new_algorithm.md",
            "How to use hooks?" => "How_to_use_hooks.md",
            "Which algorithm should I use?" => "Which_algorithm_should_I_use.md",
        ],
        "FAQ" => "FAQ.md",
        experiments,
        "Tips for Developers" => "tips.md",
        "Manual" => [
            "RLBase" => "rlbase.md",
            "RLCore" => "rlcore.md",
            "RLEnvs" => "rlenvs.md",
            "RLZoo" => "rlzoo.md",
        ],
    ],
)

postprocess_cb()


text:
using MathOptInterface
const MOI = MathOptInterface
const MOIT = MOI.Test
const MOIU = MOI.Utilities
const MOIB = MOI.Bridges

MOIU.@model(PIPSModelData,
            (),
            (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan),
            (MOI.Zeros, MOI.Nonnegatives, MOI.Nonpositives),
            (),
            (MOI.SingleVariable,),
            (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction),
            (MOI.VectorOfVariables,),
            (MOI.VectorAffineFunction,))

const optimizer = PIPS.Optimizer()
const config = MOIT.TestConfig(atol=1e-4, rtol=1e-4,
                               optimal_status=MOI.LOCALLY_SOLVED)

@testset "SolverName" begin
    @test MOI.get(optimizer, MOI.SolverName()) == "PipsNlp"
end

@testset "MOI Linear tests" begin
    exclude = ["linear8a", # Behavior in infeasible case doesn't match test.
               "linear12", # Same as above.
               "linear8b", # Behavior in unbounded case doesn't match test.
               "linear8c", # Same as above.
               "linear7",  # VectorAffineFunction not supported.
               "linear15", # VectorAffineFunction not supported.
               # new excludes
               "linear14",
                  "linear6",
                  "linear4",
                  "linear3",
                  "linear9",
                  # "linear1",
                  "linear2",
                  "linear10",
                  "linear5",
                  "linear13",
                  "linear11"
               ]
    model_for_pips = MOIU.UniversalFallback(PIPSModelData{Float64}())
    linear_optimizer = MOI.Bridges.SplitInterval{Float64}(
                         MOIU.CachingOptimizer(model_for_pips, optimizer))
    MOIT.contlineartest(linear_optimizer, config, exclude)
end

# MOI.empty!(optimizer)
#
# @testset "MOI QP/QCQP tests" begin
#     qp_optimizer = MOIU.CachingOptimizer(PIPSModelData{Float64}(), optimizer)
#     MOIT.qptest(qp_optimizer, config)
#     exclude = ["qcp1", # VectorAffineFunction not supported.
#               ]
#     MOIT.qcptest(qp_optimizer, config, exclude)
# end
#
# MOI.empty!(optimizer)
#
# @testset "MOI NLP tests" begin
#     MOIT.nlptest(optimizer, config)
# end


text:
module Py2JlRuntime
export @noscope
export PyIterator, PyVector
export @jpy_yield, @jpy_yieldfrom, jpy_literal, jpy_addlist, jpy_slice, jpy_getiter, jpy_movenext, jpy_getcurrent, jpy_call, jpy_bool, jpy_none, @jpy_all, jpy_dict, jpy_set, jpy_list, @jpy_any, jpy_add, jpy_sub, jpy_mul, jpy_floordiv, jpy_div, jpy_iadd, jpy_isub, jpy_imul, jpy_ifloordiv, jpy_idiv, jpy_pos, jpy_neg, jpy_invert, jpy_not, jpy_eq, jpy_ne, jpy_lt, jpy_le, jpy_gt, jpy_ge, jpy_isnot, jpy_is, jpy_in, jpy_notin, @jpy_conjunctive_cmp

macro noscope(ex)
    Meta.isexpr(ex, :while) || error("noscope: only use for while")
    cond = ex.args[1]
    body = ex.args[2]
    l1 = gensym("loopsetup")
    l2 = gensym("loopend")

    Expr(:block,
        Expr(:symboliclabel, l1),
        Expr(:if, cond, Expr(:block, body), Expr(:symbolicgoto, l2)),
        Expr(:symbolicgoto, l1),
        Expr(:symboliclabel, l2),
    ) |> esc
end
mutable struct PyIterator{V, E}
    self :: V
    state :: Union{E, Nothing}
end

struct PyVector{T}
    inner :: Vector{T}
end

@inline Base.iterate(x::PyVector, args...) = iterate(x.inner, args...)
@inline Base.getindex(x::PyVector, i::Integer) = getindex(x.inner, i+1)
@inline Base.setindex!(x::PyVector, v, i::Integer) = setindex!(x.inner, v, i+1)
@inline Base.push!(x::PyVector, args...) = push!(x.inner, args...)
@inline Base.pop!(x::PyVector, args...) = pop!(x.inner, args...)
@inline Base.eltype(x::PyVector) = eltype(x.inner)
@inline Base.length(x::PyVector) = length(x.inner)

macro jpy_yield(v)
    esc(:($put!(PyContinuation, $v)))
end

macro jpy_yieldfrom(v)
    iter = gensym(:iter)
    esc(quote
        $iter = $jpy_getcurrent(v)
        while $jpy_movenext($iter)
            $jpy_yield($jpy_getcurrent($iter))
        end
    end)
end

@inline function jpy_literal(x)
    x
end

Base.@pure function _infer_state_type(@nospecialize(t))
    Base.Core.Compiler._return_type(Base.iterate, (t, ))
end

@inline function jpy_addlist(x, y)
    push!(x, y)
end

@inline function jpy_slice(x, start, stop, step)
    if start === jpy_none
        start = 0
    end
    if stop === jpy_none
        stop = length(x)-1
    end
    if step === jpy_none
        step = 1
    end
    return x[start:step:stop]

end

@inline function jpy_getiter(x::T) where T
    PyIterator{T, _infer_state_type(T)}(x, nothing)
end

@inline function jpy_movenext(x::PyIterator)
    state = if x.state === nothing
        iterate(x.self)
    else
        iterate(x.self, x.state[2])
    end
    if isnothing(state)
        return false
    end
    x.state = state
    return true
end

@inline function jpy_getcurrent(x::PyIterator)
    return x.state[1]
end

@inline function jpy_call(f, args, kwargs)
    f(args...; kwargs...)
end

@inline function jpy_bool(x::Bool)
    x
end

@inline function jpy_bool(x::Union{AbstractVector, AbstractSet, AbstractDict})
    isempty(x)
end

@inline function jpy_bool(x::Tuple)
    isempty(x)
end

const jpy_none = nothing

macro jpy_all(xs...)
    foldr(xs) do (l, r)
        left = gensym(:left)
        Expr(:block,
            :($left = $l),
            :(if $jpy_bool($left)
                $r
            else
                left
            end)
        )
    end |> esc
end

@inline function jpy_dict(xs...)
    Dict(xs...)
end

@inline function jpy_set(xs...)
    Set(xs)
end

@inline function jpy_list(xs...)
    PyVector(collect(xs))
end

macro jpy_any(xs...)
    foldr(xs) do (l, r)
        left = gensym(:left)
        Expr(:block,
            :($left = $l),
            :(if $jpy_bool($left)
                $left
            else
                $r
            end))
    end |> esc
end

@inline jpy_add(a, b) = a + b
@inline jpy_sub(a, b) = a - b
@inline jpy_mul(a, b) = a * b
@inline jpy_floordiv(a, b) = div(a, b)
@inline jpy_div(a, b) = a / b

@inline jpy_iadd(a, b) = jpy_add(a, b)
@inline jpy_isub(a, b) = jpy_sub(a, b)
@inline jpy_imul(a, b) = jpy_mul(a, b)
@inline jpy_ifloordiv(a, b) = jpy_floordiv(a, b)
@inline jpy_idiv(a, b) = jpy_div(a, b)

@inline jpy_pos(a) = +a
@inline jpy_neg(a) = -a
@inline jpy_invert(a) = ~a
@inline jpy_not(a) = !jpy_bool(a)
@inline jpy_eq(a, b) = a == b
@inline jpy_ne(a, b) = a != b
@inline jpy_lt(a, b) = a < b
@inline jpy_le(a, b) = a <= b
@inline jpy_gt(a, b) = a > b
@inline jpy_ge(a, b) = a >= b
@inline jpy_isnot(a, b) = a !== b
@inline jpy_is(a, b) = a === b
@inline jpy_in(a, b) = a in b
@inline jpy_notin(a, b) = !(a in b)

macro jpy_conjunctive_cmp(op, args...)
    xs = []
    for i = 1:length(args)-1
        push!(xs, :($op($(args[i]), $(args[i+1]))))
    end
    return esc(Expr(:&&, xs...))
end

end # module


text:
using GLAbstraction, GLPlot, Reactive

window = createdisplay(eyeposition=Vec3(2), lookat=Vec3(0.5))

##########################################################
# Volume

#So far just 1 dimensional color values are supported
N 		= 128
function func(x,y,z)
    R = sqrt(x^2 + y^2+z^2)
    sin(R)/R
end
volume 	= Float32[sin(x/15f0)+sin(y/15f0)+sin(z/15f0) for x=1:N, y=1:N, z=1:N]
max 	= maximum(volume)
min 	= minimum(volume)
volume 	= (volume .- min) ./ (max .- min)

#Filter keydown events
keypressed = window.inputs[:buttonspressed]

#Make some attributes intseractive
algorithm 	= foldl( (v0, v1) -> in('I', v1) ? 2f0 : in('M', v1) ? 1f0 : v0, 2f0, keypressed) # i for isosurface, m for MIP#
isovalue 	= foldl( (v0, v1) -> in(GLFW.KEY_UP, v1) ? (v0 + 0.01f0) : (in(GLFW.KEY_DOWN, v1) ? (v0 - 0.01f0) : v0), 0.5f0, keypressed)
stepsize 	= foldl( (v0, v1) -> in(GLFW.KEY_LEFT, v1) ? (v0 + 0.0001f0) : (in(GLFW.KEY_RIGHT, v1)  ? (v0 - 0.0001f0) : v0), 0.005f0, keypressed)
lift(println, isovalue)
glplot(volume, algorithm=algorithm, isovalue=isovalue, stepsize=stepsize, color=Vec3(1,0,0))

#glplot(imread("someexample.nrrd"), algorithm=algorithm, isovalue=isovalue, stepsize=stepsize, color=Vec3(1,0,0))

renderloop(window)


text:
function SubPlane(arg0::Hyperplane, arg1::Region)
    return SubPlane((Hyperplane, Region), arg0, arg1)
end

function apply_transform(obj::AbstractSubHyperplane, arg0::Transform)
    return jcall(obj, "applyTransform", AbstractSubHyperplane, (Transform,), arg0)
end

function copy_self(obj::AbstractSubHyperplane)
    return jcall(obj, "copySelf", AbstractSubHyperplane, ())
end

function get_hyperplane(obj::AbstractSubHyperplane)
    return jcall(obj, "getHyperplane", Hyperplane, ())
end

function get_remaining_region(obj::AbstractSubHyperplane)
    return jcall(obj, "getRemainingRegion", Region, ())
end

function get_size(obj::AbstractSubHyperplane)
    return jcall(obj, "getSize", jdouble, ())
end

function is_empty(obj::AbstractSubHyperplane)
    return jcall(obj, "isEmpty", jboolean, ())
end

function reunite(obj::AbstractSubHyperplane, arg0::SubHyperplane)
    return jcall(obj, "reunite", AbstractSubHyperplane, (SubHyperplane,), arg0)
end

function split(obj::SubPlane, arg0::Hyperplane)
    return jcall(obj, "split", SubHyperplane_SplitSubHyperplane, (Hyperplane,), arg0)
end


text:
module RulesAdditionIn1Test

using Test
using ReactiveMP
using Random
import ReactiveMP: @test_rules

@testset "rules:typeof(+):in1" begin

     @testset "Belief Propagation: (m_out::PointMass, m_in2::PointMass)" begin

        @test_rules [ with_float_conversions = true ] typeof(+)(:in1, Marginalisation) [
            (input = (m_out = PointMass(1.0),                   m_in2 = PointMass(-1.0)),                    output = PointMass(2.0)),
            (input = (m_out = PointMass([ 1.0 ]),               m_in2 = PointMass([ -2.0 ])),                output = PointMass([ 3.0 ])),
            (input = (m_out = PointMass([ 1.0 2.0; 3.0 4.0 ]),  m_in2 = PointMass([-2.0 -1.0; -4.0 -1.0])),  output = PointMass([ 3.0 3.0; 7.0 5.0 ])),
        ]

    end

     @testset "Belief Propagation: (m_out::UnivariateNormalDistributionsFamily, m_in2::PointMass)" begin

        @test_rules [ with_float_conversions = true ] typeof(+)(:in1, Marginalisation) [
            (input = (m_out = NormalMeanVariance(1.0, 2.0),   m_in2 = PointMass(2.0)),  output = NormalMeanVariance(-1.0, 2.0)),
            (input = (m_out = NormalMeanVariance(-1.0, 3.0),  m_in2 = PointMass(-3.0)), output = NormalMeanVariance(2.0, 3.0)),
            (input = (m_out = NormalMeanPrecision(4.0, 7.0),  m_in2 = PointMass(1.0)),  output = NormalMeanPrecision(3.0, 7.0)),
            (input = (m_out = NormalMeanPrecision(-1.0, 2.0), m_in2 = PointMass(-3.0)), output = NormalMeanPrecision(2.0, 2.0)),
            (input = (m_out = NormalWeightedMeanPrecision(4.0, 2.0), m_in2 = PointMass(1.0)),  output = NormalMeanVariance(1.0, 1/2)),
            (input = (m_out = NormalWeightedMeanPrecision(8.0, 4.0), m_in2 = PointMass(-3.0)), output = NormalMeanVariance(5.0, 1/4)),
        ]

    end

        @testset "Belief Propagation: (m_out::PointMass, m_in2::UnivariateNormalDistributionsFamily)" begin

        @test_rules [ with_float_conversions = true ] typeof(+)(:in1, Marginalisation) [
            (input = (m_out = PointMass(1.0),  m_in2 = NormalMeanVariance(1.0, 2.0)),   output = NormalMeanVariance(0.0, 2.0)),
            (input = (m_out = PointMass(-3.0), m_in2 = NormalMeanVariance(-1.0, 3.0)),  output = NormalMeanVariance(-2.0, 3.0)),
            (input = (m_out = PointMass(1.0),  m_in2 = NormalMeanPrecision(4.0, 7.0)),  output = NormalMeanPrecision(-3.0, 7.0)),
            (input = (m_out = PointMass(-3.0), m_in2 = NormalMeanPrecision(-1.0, 2.0)), output = NormalMeanPrecision(-2.0, 2.0)),
            (input = (m_out = PointMass(1.0),  m_in2 = NormalWeightedMeanPrecision(4.0, 2.0)),  output = NormalMeanVariance(-1.0, 1/2)),
            (input = (m_out = PointMass(-3.0), m_in2 = NormalWeightedMeanPrecision(8.0, 4.0)),  output = NormalMeanVariance(-5.0, 1/4)),
        ]

    end

        @testset "Belief Propagation: (m_out::MultivariateNormalDistributionsFamily, m_in2::PointMass)" begin

        @test_rules [ with_float_conversions = true ] typeof(+)(:in1, Marginalisation) [
            (input = (m_out = MvNormalMeanCovariance([1.0, 3.0], [3.0 2.0; 2.0 4.0]),  m_in2 = PointMass([1.0, 1.0])),  output = MvNormalMeanCovariance([0.0, 2.0], [3.0 2.0; 2.0 4.0])),
            (input = (m_out = MvNormalMeanCovariance([-4.0, 3.0], [3.0 2.0; 2.0 4.0]), m_in2 = PointMass([1.0, 1.0])),  output = MvNormalMeanCovariance([-5.0, 2.0], [3.0 2.0; 2.0 4.0])),
            (input = (m_out = MvNormalMeanPrecision([-4.0, 3.0], [3.0 2.0; 2.0 4.0]),  m_in2 = PointMass([3.0, 2.0])),  output = MvNormalMeanPrecision([-7.0, 1.0], [3.0 2.0; 2.0 4.0])),
            (input = (m_out = MvNormalMeanPrecision([-4.0, 3.0], [3.0 2.0; 2.0 4.0]),  m_in2 = PointMass([-2.0, 3.0])), output = MvNormalMeanPrecision([-2.0, 0.0], [3.0 2.0; 2.0 4.0])),
            (input = (m_out = MvNormalWeightedMeanPrecision([1.0, 3.0], [ 1.0 0.0; 0.0 1.0 ]), m_in2 = PointMass([ 2.0, 7.0 ])),  output = MvNormalMeanCovariance([ -1.0, -4.0 ], [ 1.0 0.0; 0.0 1.0 ])),
            (input = (m_out = MvNormalWeightedMeanPrecision([2.0, 4.0], [ 2.0 0.0; 0.0 2.0 ]), m_in2 = PointMass([ 1.0, 3.0 ])),  output = MvNormalMeanCovariance([ 0.0, -1.0 ], [ 1/2 0.0; 0.0 1/2 ])),
        ]

    end

    @testset "Belief Propagation: (m_out::PointMass, m_in2::MultivariateNormalDistributionsFamily)" begin

        @test_rules [ with_float_conversions = true ] typeof(+)(:in1, Marginalisation) [
            (input = (m_out = PointMass([1.0, 1.0]),    m_in2 = MvNormalMeanCovariance([1.0, 3.0], [3.0 2.0; 2.0 4.0])),    output = MvNormalMeanCovariance([0.0, -2.0], [3.0 2.0; 2.0 4.0])),
            (input = (m_out = PointMass([1.0, 1.0]),    m_in2 = MvNormalMeanCovariance([-4.0, 3.0], [3.0 2.0; 2.0 4.0])),   output = MvNormalMeanCovariance([5.0, -2.0], [3.0 2.0; 2.0 4.0])),
            (input = (m_out = PointMass([1.0, 1.0]),    m_in2 = MvNormalMeanPrecision([-4.0, 3.0], [3.0 2.0; 2.0 4.0])),    output = MvNormalMeanPrecision([5.0, -2.0], [3.0 2.0; 2.0 4.0])),
            (input = (m_out = PointMass([-2.0, 1.0]),   m_in2 = MvNormalMeanPrecision([-4.0, 3.0], [3.0 2.0; 2.0 4.0])),    output = MvNormalMeanPrecision([2.0, -2.0], [3.0 2.0; 2.0 4.0])),
            (input = (m_out = PointMass([ 2.0, 7.0 ]),  m_in2 = MvNormalWeightedMeanPrecision([1.0, 3.0], [ 1.0 0.0; 0.0 1.0 ])),  output = MvNormalMeanCovariance([ 1.0, 4.0 ], [ 1.0 0.0; 0.0 1.0 ])),
            (input = (m_out = PointMass([ 1.0, 3.0 ]),  m_in2 = MvNormalWeightedMeanPrecision([2.0, 4.0], [ 2.0 0.0; 0.0 2.0 ])),  output = MvNormalMeanCovariance([ 0.0, 1.0 ], [ 1/2 0.0; 0.0 1/2 ])),
        ]

    end

    @testset "Belief Propagation: (m_out::UnivariateNormalDistributionsFamily, m_in2::UnivariateNormalDistributionsFamily)" begin

        @test_rules [ with_float_conversions = true ] typeof(+)(:in1, Marginalisation) [
            (input = (m_out = NormalMeanVariance(1.0, 2.0),             m_in2 = NormalMeanVariance(3.0, 4.0)),          output = NormalMeanVariance(-2.0, 6.0)),
            (input = (m_out = NormalMeanVariance(-1.0, 2.0),            m_in2 = NormalMeanVariance(-2.0, 3.0)),         output = NormalMeanVariance(1.0, 5.0)),
            (input = (m_out = NormalMeanPrecision(2.0, 2.0),            m_in2 = NormalMeanPrecision(-1.0, 3.0)),        output = NormalMeanVariance(3.0, (2.0 + 3.0) / (2.0*3.0))),
            (input = (m_out = NormalMeanPrecision(-1.0, 2.0),           m_in2 = NormalMeanPrecision(-1.0, 3.0)),        output = NormalMeanVariance(0.0, (2.0 + 3.0) / (2.0*3.0))),
            (input = (m_out = NormalMeanPrecision(2.0, 2.0),            m_in2 = NormalMeanVariance(-1.0, 3.0)),         output = NormalMeanVariance(3.0, 3.5)),
            (input = (m_out = NormalWeightedMeanPrecision(8.0, 4.0),    m_in2 = NormalMeanVariance(-3.0, 1.0)),         output = NormalMeanVariance(5.0, 5/4)),
            (input = (m_out = NormalMeanVariance(-4.0, 2.0),            m_in2 = NormalWeightedMeanPrecision(6.0, 3.0)), output = NormalMeanVariance(-6.0, 7/3)),
        ]

    end

    @testset "Belief Propagation: (m_out::MultivariateNormalDistributionsFamily, m_in2::MultivariateNormalDistributionsFamily)" begin

        @test_rules [ with_float_conversions = true ] typeof(+)(:in1, Marginalisation) [
            (input = (m_out = MvNormalMeanCovariance([1.0, 3.0], [3.0 2.0; 2.0 4.0]),   m_in2 = MvNormalMeanCovariance([1.0, 3.0], [3.0 2.0; 2.0 4.0])),    output = MvNormalMeanCovariance([0.0, 0.0], [6.0 4.0; 4.0 8.0])),
            (input = (m_out = MvNormalMeanCovariance([-1.0, 1.0], [3.0 2.0; 2.0 4.0]),  m_in2 = MvNormalMeanCovariance([1.0, 3.0], [3.0 2.0; 2.0 4.0])),    output = MvNormalMeanCovariance([-2.0, -2.0], [6.0 4.0; 4.0 8.0])),
            (input = (m_out = MvNormalMeanPrecision([1.0, 0.0], [2.0 0.0; 0.0 1.0]),    m_in2 = MvNormalMeanPrecision([1.0, -7.0], [2.0 0.0; 0.0 3.0])),    output = MvNormalMeanCovariance([0.0, 7.0], [1.0 0.0; 0.0 4/3])),
            (input = (m_out = MvNormalMeanCovariance([1.0, -1.0], [3.0 1.0; 1.0 4.0]),  m_in2 = MvNormalMeanPrecision([1.0, 4.0], [2.0 1.0; 1.0 3.0])),     output = MvNormalMeanCovariance([0.0, -5.0], [36/10 4/5; 4/5 44/10])),
            (input = (m_out = MvNormalMeanPrecision([1.0, 0.0], [2.0 0.0; 0.0 1.0]),    m_in2 = MvNormalWeightedMeanPrecision([1.0, -7.0], [2.0 0.0; 0.0 3.0])),    output = MvNormalMeanCovariance([1/2, 7/3], [1.0 0.0; 0.0 4/3])),
            (input = (m_out = MvNormalMeanCovariance([1.0, 1.0], [3.0 1.0; 1.0 4.0]),   m_in2 = MvNormalWeightedMeanPrecision([1.0, 4.0], [2.0 1.0; 1.0 3.0])),     output = MvNormalMeanCovariance([6/5, -2/5], [36/10 4/5; 4/5 44/10])),
            (input = (m_out = MvNormalWeightedMeanPrecision([1.0, -7.0], [2.0 0.0; 0.0 3.0]),   m_in2 = MvNormalMeanPrecision([1.0, 0.0], [2.0 0.0; 0.0 1.0])),     output = MvNormalMeanCovariance([-1/2, -7/3], [1.0 0.0; 0.0 4/3])),
            (input = (m_out = MvNormalWeightedMeanPrecision([1.0, 4.0], [2.0 1.0; 1.0 3.0]),    m_in2 = MvNormalMeanCovariance([1.0, 1.0], [3.0 1.0; 1.0 4.0])),    output = MvNormalMeanCovariance([-6/5, 2/5], [36/10 4/5; 4/5 44/10])),
        ]

    end
end
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule AFLplusplus_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("AFLplusplus")
JLLWrappers.@generate_main_file("AFLplusplus", UUID("2a46d845-6c66-546d-9393-effbe7e72607"))
end  # module AFLplusplus_jll


text:
export pradius, pradiusb

function pradiusk(As, p, k, pnorm)
    mean(map(A->opnorm(A, pnorm)^p, As))^(1/(p*k))
end

function pradius(s::DiscreteSwitchedLinearSystem, p, lift::Function; pnorm=Inf, ɛ=1e-2, forceub=false)
    ρpmpp = ρ(sum((rm -> lift(rm.A, p)).(s.resetmaps)))
    if forceub
        ρpmpp^(1/p)
    else
        (ρpmpp / ρA(s))^(1/p)
    end
end

function kozyakinlift(s, t, As)
    n = nstates(s)
    A = As[symbol(s, t)]
    kron(sparse([source(s, t)], [target(s, t)], [one(eltype(A))], n, n), A)
end
function pradius(s::ConstrainedDiscreteSwitchedLinearSystem, p, lift::Function; pnorm=Inf, ɛ=1e-2, forceub=false)
    Alifted = (rm -> lift(rm.A, p)).(s.resetmaps)
    ρpmpp = ρ(sum(kozyakinlift(s, t, Alifted) for t in transitions(s)))
    if forceub
        ρpmpp^(1/p)
    else
        (ρpmpp / ρA(s))^(1/p)
    end
end

abstract type PRadiusAlgorithm end
abstract type ExactPRadiusAlgorithm <: PRadiusAlgorithm end
function checkeven(p, algo::ExactPRadiusAlgorithm)
    if !iseven(p)
        throw(ArgumentError("Odd p is not supported yet for pradius computation with $algo algorithm"))
    end
end

export VeroneseLift, KroneckerLift, BruteForce
struct VeroneseLift <: ExactPRadiusAlgorithm end
struct KroneckerLift <: ExactPRadiusAlgorithm end
struct BruteForce <: PRadiusAlgorithm end

liftfunction(::VeroneseLift) = veroneselift
liftfunction(::KroneckerLift) = kroneckerlift
function kroneckerlift(A::AbstractMatrix, p::Integer)
    @assert iseven(p)
    kronpow(veroneselift(A, 2), div(p, 2))
end

function pradius(s::AbstractDiscreteSwitchedSystem, p, algo::ExactPRadiusAlgorithm=VeroneseLift(); kws...)
    checkeven(p, algo)
    pradius(s, p, liftfunction(algo); kws...)
end
# The Inf- and 1-norm of a matrices are easier to compute than the 2-norm
function pradius(s::DiscreteSwitchedLinearSystem, p, ::BruteForce; pnorm=Inf, ɛ=1e-2, forceub=false)
    As = map(rm -> rm.A, s.resetmaps)
    Ascur = As
    ρp = Float64[]
    m = ntransitions(s)
    k = 0
    while k < 2 || !isapprox(ρp[end-1], ρp[end], rtol=ɛ)
        k += 1
        l = length(Ascur)
        Asnew = similar(Ascur, l * m)
        for (i,A) in enumerate(As)
            for (j,B) in enumerate(Ascur)
                Asnew[(i-1)*l+j] = A*B
            end
        end
        Ascur = Asnew
        push!(ρp, pradiusk(Ascur, p, k, pnorm))
    end
    if forceub
        ρp[end] * length(As)^(1/p)
    else
        ρp[end]
    end
end

# p-radius algo needs to be exact to compute bounds on JSR
function pradiusb(s::AbstractDiscreteSwitchedSystem, p, algo::ExactPRadiusAlgorithm=VeroneseLift())
    ρpmp = pradius(s, p, algo, forceub=true)
    ρp = ρpmp / ρA(s)^(1/p)
    updateb!(s, ρp, ρpmp)
end


text:
export SparseBinMatrixCSR, ParallelBinCSR
using Compat

###################### SparseBinMatrixCSR #######################

type SparseBinMatrixCSR
  m::Int
  n::Int
  col_ind::Vector{Int32}
  row_ptr::Vector{Int32}
end

@compat type ParallelBinCSR
  m::Int
  n::Int
  pids::Vector{Int64}
  csrs::Vector{Future}
  mranges::Vector{StepRange{Int64,Int64}}
  blocksize::Int64
end

function SparseBinMatrixCSR(rows::Vector{Int32}, cols::Vector{Int32})
  rsorted = sortperm(rows)
  m = convert(Int, maximum(rows))
  n = convert(Int, maximum(cols))
  row_ptr = zeros(Int32, m+1)
  row_ptr[1:end] = length(rows)+1
  rows2   = rows[rsorted]
  prev = zero(Int32)
  for i = one(Int32):convert(Int32,length(rows))
    while rows2[i] > prev
      prev += one(Int32)
      row_ptr[prev] = i
    end
  end
  return SparseBinMatrixCSR(m, n, cols[rsorted], row_ptr)
end

import Base.size
size(A::SparseBinMatrixCSR) = (A.m, A.n)
size(X::SparseBinMatrixCSR, d::Int) = d==1 ? X.m : X.n

import Base.show
function show(io::IO, csr::SparseBinMatrixCSR)
  println(io, "$(csr.m) x $(csr.m) binary sparse matrix (CSR) with $(length(csr.col_ind)) entries.")
end

import Base.A_mul_B!
function A_mul_B!{Tx}(y::AbstractArray{Tx,1}, A::SparseBinMatrixCSR, x::AbstractArray{Tx,1})
    A.n == length(x) || throw(DimensionMismatch("A.n=$(A.n) must equal length(x)=$(length(x))"))
    A.m == length(y) || throw(DimensionMismatch("A.m=$(A.m) must equal length(y)=$(length(y))"))
    zro = zero(Tx)
    row_ptr = A.row_ptr
    col_ind = A.col_ind
    @inbounds for row = 1:A.m
        tmp = zro
        for i = row_ptr[row]:row_ptr[row+1]-1
          tmp += x[col_ind[i]]
        end
        y[row] = tmp
    end
    return
end

######## ParallelBinCSR #########
function ParallelBinCSR(rows::Vector{Int32}, cols::Vector{Int32}, pids::Vector{Int}, blocksize=64)
  csr = SparseBinMatrixCSR(rows, cols)
  npids  = length(pids)
  ranges = StepRange{Int,Int}[ 1+(i-1)*blocksize:npids*blocksize:size(csr,1) for i in 1:npids ]
  @compat pcsr = ParallelBinCSR(csr.m, csr.n, pids, Future[], ranges, blocksize)
  for i in 1:npids
    ref = @spawnat pids[i] fetch(csr)
    push!(pcsr.csrs, ref)
  end
  return pcsr
end

function A_mul_B!{Tx}(y::SharedArray{Tx,1}, A::ParallelBinCSR, x::SharedArray{Tx,1})
  A.n == length(x) || throw(DimensionMismatch("A.n=$(A.n) must equal length(x)=$(length(x))"))
  A.m == length(y) || throw(DimensionMismatch("A.m=$(A.m) must equal length(y)=$(length(y))"))

  @sync begin
    for p in 1:length(A.pids)
      pid = A.pids[p]
      if pid != myid() || np == 1
        @async begin
          remotecall_wait(A_mul_B_part_ref, pid, y, A.csrs[p], x, A.mranges[p], A.blocksize)
        end
      end
    end
  end
  return nothing
end

@compat function A_mul_B_part_ref{Tx}(y::SharedArray{Tx,1}, Aref::Future, x::SharedArray{Tx,1}, range::StepRange{Int,Int}, blocksize::Int)
  A = fetch(Aref)::SparseBinMatrixCSR
  A_mul_B_range!(y, A, x, range, blocksize)
  return nothing
end

function A_mul_B_range!{Tx}(y::AbstractArray{Tx,1}, A::SparseBinMatrixCSR, x::AbstractArray{Tx,1}, range, blocksize::Int)
    A.n == length(x) || throw(DimensionMismatch("A.n=$(A.n) must equal length(x)=$(length(x))"))
    A.m == length(y) || throw(DimensionMismatch("A.m=$(A.m) must equal length(y)=$(length(y))"))
    zro = zero(Tx)
    row_ptr = A.row_ptr
    col_ind = A.col_ind
    @inbounds for b0 = range, row=b0:min(b0+blocksize-1, A.m)
        tmp = zro
        for i = row_ptr[row]:row_ptr[row+1]-1
          tmp += x[col_ind[i]]
        end
        y[row] = tmp
    end
    return
end


text:
# [SOREN] TRUST REGION MODIFICATION OF NEWTON'S METHOD
# [N&W] Numerical optimization
# [Yuan] A review of trust region algorithms for optimization
abstract type TRSPSolver end
abstract type NearlyExactTRSP <: TRSPSolver end

include("solvers/NWI.jl")
include("solvers/Dogleg.jl")
include("solvers/NTR.jl")
include("solvers/TCG.jl")
#include("subproblemsolvers/TRS.jl") just make an example instead of relying onTRS.jl

function tr_return(; λ, ∇f, H, s, interior, solved, hard_case, Δ, m = nothing)
    m = m isa Nothing ? dot(∇f, s) + dot(s, H * s) / 2 : m
    (
        p = s,
        mz = m,
        interior = interior,
        λ = λ,
        hard_case = hard_case,
        solved = solved,
        Δ = Δ,
    )
end

function update_H!(H, h, λ = nothing)
    T = eltype(h)
    n = length(h)
    if !(λ == T(0))
        for i = 1:n
            @inbounds H[i, i] = λ isa Nothing ? h[i] : h[i] + λ
        end
    end
    H
end

max_stars_repo_name: JuliaNLSolvers/NLSolvers

max_issues_repo_name: JuliaNLSolvers/NLSolvers

max_forks_repo_name: JuliaNLSolvers/NLSolvers


text:
using Includev
using Test
using GR

testfiletoinc = joinpath(@__DIR__,"plotrandomwalk.jl")
includev(testfiletoinc)


text:
function split_merge_sample!(
        model::SeqModel,
        spikes::Vector{Spike},
        num_samples::Int64,
        assignments::Vector{Int64},
        split_merge_window::Float64
    )

    # Return early if no samples.
    num_samples == 0 && return assignments

    # Arrays holding indices for split clusters.
    A = Int64[]
    B = Int64[]
    C = Int64[]

    # Spike indices not assigned to the background.
    non_bkgd_idx = findall(c->(c != -1), assignments)

    # Return early if everything is assigned to the background.
    length(non_bkgd_idx) < 3 && return assignments

    # Concentration parameter.
    α = model.priors.seq_event_amplitude.α

    for s = 1:num_samples

        # Sample random spike.
        i = rand(non_bkgd_idx)
        ti = spikes[i].timestamp

        # Sample another random spike, within 'split_merge_window' time
        # units from spike i. Do this by placing all eligible spikes in
        # set 'A' and then sampling a spike index j from A.
        for k in non_bkgd_idx
            tk = spikes[k].timestamp
            if (i != k) && (abs(ti - tk) < split_merge_window)
                push!(A, k)
            end
        end

        # Failed to identify any valid proposal.
        isempty(A) && continue

        # Sample second spike and empty set A.
        j = rand(A)
        empty!(A)

        # Cluster assignment for spike i and spike j.
        ci = assignments[i]
        cj = assignments[j]

        # Propose split.
        if ci == cj

            # println("proposing split...")

            # Put spikes[i] in set A, spikes[j] in set B.
            push!(A, i)
            push!(B, j)

            # Randomly assign other spikes in cluster to A or B.
            for k in non_bkgd_idx
                if (k != i) && (k != j) && (assignments[k] == ci)
                    if rand() > 0.5
                        push!(A, k)
                    else
                        push!(B, k)
                    end
                end
            end

            # Create combined list of all spikes.
            append!(C, A)
            append!(C, B)

            # Compute log acceptance probability.
            accept_log_prob = (
                event_log_like(model, view(spikes, A))
                + event_log_like(model, view(spikes, B))
                - event_log_like(model, view(spikes, C))
                + lgamma(length(A) + α)
                + lgamma(length(B) + α)
                - lgamma(length(C) + α)
                - lgamma(α)
                + model.new_cluster_log_prob
                - (length(C) - 2) * log(0.5)
            )
            # println("  accept_log_prob : ", accept_log_prob)

            # Accept move.
            if rand() < exp(accept_log_prob)

                # println("!!!!!! SPLIT ACCEPTED !!!!!!!!!!")
                # println("num events before = ", length(model.sequence_events))

                # Create new event for spike j and move it there.
                remove_datapoint!(model, spikes[j], assignments[j])
                k = add_event!(model, spikes[j])
                assignments[j] = k

                # Move other spikes in set B to the new event at index k.
                for b in B
                    if b != j
                        remove_datapoint!(model, spikes[b], assignments[b])
                        assignments[b] = add_datapoint!(model, spikes[b], k)
                    end
                end

                # @assert model.sequence_events[ci].spike_count == length(A)
                # @assert model.sequence_events[k].spike_count == length(B)
                # @assert sum(assignments .== ci) == length(A)
                # @assert sum(assignments .== k) == length(B)

                # println("num events after = ", length(model.sequence_events))

                # Spike i and other spikes in set A stay where they are now.
            end

            # Empty arrays for future split moves.
            empty!(A)
            empty!(B)
            empty!(C)

        # Propose merge.
        else

            # println("proposing merge...")

            # Put spikes[i] in A, spikes[j] in B.
            push!(A, i)
            push!(B, j)

            # Put spikes where (assignment == ci) in A.
            # Put spikes where (assignment == cj) in B.
            for k in non_bkgd_idx
                if (k != i) && (k != j)
                    if (assignments[k] == ci)
                        push!(A, k)
                    elseif (assignments[k] == cj)
                        push!(B, k)
                    end
                end
            end

            # Create combined list of all spikes.
            append!(C, A)
            append!(C, B)

            # Compute log acceptance probability.
            accept_log_prob = (
                event_log_like(model, view(spikes, C))
                - event_log_like(model, view(spikes, A))
                - event_log_like(model, view(spikes, B))
                + lgamma(length(C) + α)
                + lgamma(α)
                - lgamma(length(A) + α)
                - lgamma(length(B) + α)
                - model.new_cluster_log_prob
                + (length(C) - 2) * log(0.5)
            )
            # println("  accept_log_prob : ", accept_log_prob)

            # Accept move.
            if rand() < exp(accept_log_prob)

                # println("!!!!!! MERGE ACCEPTED !!!!!!!!!!")
                # println("num events before = ", length(model.sequence_events))

                # Move spike j and others in set B to cluster ci.
                for b in B
                    remove_datapoint!(model, spikes[b], assignments[b])
                    assignments[b] = add_datapoint!(model, spikes[b], ci)
                end

                # @assert model.sequence_events[cj].spike_count == 0
                # @assert sum(assignments .== cj) == 0
                # @assert sum(assignments .== ci) == length(C)

                # println("num events after = ", length(model.sequence_events))

                # Spike i and other spikes in set A stay where they are now.
            end

            # Empty arrays for future merge moves.
            empty!(A)
            empty!(B)
            empty!(C)

        end

    end

    return assignments

end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ExportTemplateRequest <: SwaggerModel
    resources::Any # spec type: Union{ Nothing, Vector{String} } # spec name: resources
    options::Any # spec type: Union{ Nothing, String } # spec name: options

    function ExportTemplateRequest(;resources=nothing, options=nothing)
        o = new()
        validate_property(ExportTemplateRequest, Symbol("resources"), resources)
        setfield!(o, Symbol("resources"), resources)
        validate_property(ExportTemplateRequest, Symbol("options"), options)
        setfield!(o, Symbol("options"), options)
        o
    end
end # type ExportTemplateRequest

const _property_map_ExportTemplateRequest = Dict{Symbol,Symbol}(Symbol("resources")=>Symbol("resources"), Symbol("options")=>Symbol("options"))
const _property_types_ExportTemplateRequest = Dict{Symbol,String}(Symbol("resources")=>"Vector{String}", Symbol("options")=>"String")
Base.propertynames(::Type{ ExportTemplateRequest }) = collect(keys(_property_map_ExportTemplateRequest))
Swagger.property_type(::Type{ ExportTemplateRequest }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ExportTemplateRequest[name]))}
Swagger.field_name(::Type{ ExportTemplateRequest }, property_name::Symbol) =  _property_map_ExportTemplateRequest[property_name]

function check_required(o::ExportTemplateRequest)
    true
end

function validate_property(::Type{ ExportTemplateRequest }, name::Symbol, val)
end


text:
@testset "linker" begin

Context() do ctx
Builder(ctx) do builder
    mod1 = let
        mod = LLVM.Module("SomeModule"; ctx)
        ft = LLVM.FunctionType(LLVM.VoidType(ctx))
        fn = LLVM.Function(mod, "SomeFunction", ft)

        entry = BasicBlock(fn, "entry"; ctx)
        position!(builder, entry)

        ret!(builder)

        mod
    end

    mod2 = let
        mod = LLVM.Module("SomeOtherModule"; ctx)
        ft = LLVM.FunctionType(LLVM.VoidType(ctx))
        fn = LLVM.Function(mod, "SomeOtherFunction", ft)

        entry = BasicBlock(fn, "entry"; ctx)
        position!(builder, entry)

        ret!(builder)

        mod
    end

    link!(mod1, mod2)
    @test haskey(functions(mod1), "SomeFunction")
    @test haskey(functions(mod1), "SomeOtherFunction")
end
end

end


text:
function openhandle(f, h::Handle)
    f(h)
    @wraperror LibClingo.clingo_solve_handle_close(h.p)
end

function getmodel(h::Handle)
    @wraperror LibClingo.clingo_solve_handle_resume(h.p)
    r = Ref(Ptr{LibClingo.clingo_model}())
    @wraperror LibClingo.clingo_solve_handle_model(h.p,r)
    m = Model(r[])
    return (m.p == C_NULL ? nothing : m)
end

function getmodelnumber(m::Model)
    n = Ref{Csize_t}(0)
    @wraperror LibClingo.clingo_model_number(m.p, n)
    return n[]
end
function getresult(h::Handle)
    r = Ref{LibClingo.clingo_solve_result_bitset_t}(0)
    @wraperror LibClingo.clingo_solve_handle_get(h.p, r)

    return LibClingo.clingo_solve_result_e(r[])
end
function modeltostrings(m::Model)
    numatoms = Ref{Csize_t}(0)
    @wraperror LibClingo.clingo_model_symbols_size(m.p, LibClingo.clingo_show_type_shown, numatoms)
    numatoms = numatoms[]
    atoms = Vector{LibClingo.clingo_symbol_t}(undef, numatoms)
    @wraperror LibClingo.clingo_model_symbols(m.p, LibClingo.clingo_show_type_shown, atoms, numatoms)
    symbols = Vector{String}()
    for i ∈ 1:numatoms
        strsize = Ref{Csize_t}(0)
        @wraperror LibClingo.clingo_symbol_to_string_size(atoms[i], strsize)
        strsize = strsize[]
        str = Vector{Cchar}(undef, strsize)
        @wraperror LibClingo.clingo_symbol_to_string(atoms[i], str, strsize)
        push!(symbols, unsafe_string(pointer(str)))
    end
    return symbols
end
function modelinfo(m::Model)
    d = Dict{Symbol,Any}()
    d[:number] = Int(getmodelnumber(m))
    s = modeltostrings(m)
    d[:symbols] = s
    return d
end

function storemodels(handle::Handle, v::Vector{Dict{Symbol, Any}})
    while true
        m = getmodel(h)
        m !== nothing || break
        push!(v,modelinfo(m))
    end
    return v
end


text:
using ThreeBodyProblem
using LinearAlgebra
# This particular state is from a meteor we sensed
r_enu = [191.0973115173407; 78.7745585024403; 57.3291938159787] # [km]
v_enu = [-52.199447002937305; -21.517740337533645; -15.659828241680069] #[km/s]
rv_enu = [r_enu; v_enu]
# coordinates of observer (ALTAIR radar)
ϕ = 9.3987 # [deg] latitude
λ = 167.482 # [deg] longitude
h = 0.083532 # [km] altitude
ϕ, λ, h

rv_ecef = enu2ecef(rv_enu, ϕ, λ, h, geodetic=true, ang_unit=:deg)

# Let's check that we get the same r_enu vector when working backwards
error_enu_ecef = norm(rv_enu - ecef2enu(rv_ecef, ϕ, λ, h))

# UT1 November 17, 2000 15:00:05.634
date = [2000,11,17,15,0,5.64] #[Y,M,D,h,m,s] observation date
mjd = date2mjd(date) # [days] modified julian date
θ = mjd2gmst(mjd) # [rad] Greenwich Mean Sidereal time

rv_eci = ecef2eci(rv_ecef, θ, ang_unit=:rad)
# Let's check that we get the same r_ecef vector when working backwards
error_ecef_eci = norm(rv_ecef - eci2ecef(rv_eci, θ))

sys = sun_earth()
a,e,i,Ω,ω,ν,Π,u,l,ℰ = cart2oe(rv_eci, 3.986004328969392e+05)

# Let's check that we get the same r_eci vector when working backwards
error_eci_oe = norm(rv_eci - oe2cart(a,e,i,Ω,ω,ν,G*sys.sec.m))


text:
function NelderMeadSimplex(arg0::Vector{Vector{jdouble}})
    return NelderMeadSimplex((Vector{Vector{jdouble}},), arg0)
end

function NelderMeadSimplex(arg0::Vector{Vector{jdouble}}, arg1::jdouble, arg2::jdouble, arg3::jdouble, arg4::jdouble)
    return NelderMeadSimplex((Vector{Vector{jdouble}}, jdouble, jdouble, jdouble, jdouble), arg0, arg1, arg2, arg3, arg4)
end

function NelderMeadSimplex(arg0::Vector{jdouble})
    return NelderMeadSimplex((Vector{jdouble},), arg0)
end

function NelderMeadSimplex(arg0::Vector{jdouble}, arg1::jdouble, arg2::jdouble, arg3::jdouble, arg4::jdouble)
    return NelderMeadSimplex((Vector{jdouble}, jdouble, jdouble, jdouble, jdouble), arg0, arg1, arg2, arg3, arg4)
end

function NelderMeadSimplex(arg0::jint)
    return NelderMeadSimplex((jint,), arg0)
end

function NelderMeadSimplex(arg0::jint, arg1::jdouble)
    return NelderMeadSimplex((jint, jdouble), arg0, arg1)
end

function NelderMeadSimplex(arg0::jint, arg1::jdouble, arg2::jdouble, arg3::jdouble, arg4::jdouble)
    return NelderMeadSimplex((jint, jdouble, jdouble, jdouble, jdouble), arg0, arg1, arg2, arg3, arg4)
end

function NelderMeadSimplex(arg0::jint, arg1::jdouble, arg2::jdouble, arg3::jdouble, arg4::jdouble, arg5::jdouble)
    return NelderMeadSimplex((jint, jdouble, jdouble, jdouble, jdouble, jdouble), arg0, arg1, arg2, arg3, arg4, arg5)
end

function iterate(obj::NelderMeadSimplex, arg0::MultivariateFunction, arg1::Comparator)
    return jcall(obj, "iterate", void, (MultivariateFunction, Comparator), arg0, arg1)
end