text:
@testset "distributions" begin
    Random.seed!(1234)

    # Create random vectors and matrices
    dim = 3
    a = rand(dim)
    b = rand(dim)
    c = rand(dim)
    A = rand(dim, dim)
    B = rand(dim, dim)
    C = rand(dim, dim)

    # Create random numbers
    alpha = rand()
    beta = rand()
    gamma = rand()

    # Create matrix `X` such that `X` and `I - X` are positive definite if `A ≠ 0`.
    function to_beta_mat(A)
        S = A * A' + I
        invL = inv(cholesky(S).L)
        return invL * invL'
    end

    # Create positive values.
    to_positive(x) = exp.(x)
    to_positive(x::AbstractArray{<:AbstractArray}) = to_positive.(x)

    # The following definition should not be needed
    # It seems there is a bug in the default `rand_tangent` that causes a
    # StackOverflowError though
    function ChainRulesTestUtils.rand_tangent(::Random.AbstractRNG, ::typeof(to_positive))
        return NoTangent()
    end

    # Tests that have a `broken` field can be executed but, according to FiniteDifferences,
    # fail to produce the correct result. These tests can be checked with `@test_broken`.
    univariate_distributions = DistSpec[
        ## Univariate discrete distributions

        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(TuringPoissonBinomial, ([0.5, 0.5],), 0),
        DistSpec(TuringPoissonBinomial, ([0.5, 0.5],), [0, 0]),

        ## Univariate continuous distributions

        DistSpec(Arcsine, (), 0.5),
        DistSpec(Arcsine, (1.0,), 0.5),
        DistSpec(Arcsine, (0.0, 2.0), 0.5),

        DistSpec(Beta, (), 0.4),
        DistSpec(Beta, (1.5,), 0.4),
        DistSpec(Beta, (1.5, 2.0), 0.4),

        DistSpec(BetaPrime, (), 0.4),
        DistSpec(BetaPrime, (1.5,), 0.4),
        DistSpec(BetaPrime, (1.5, 2.0), 0.4),

        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(Chernoff, (), 0.5, broken=(:Zygote,)),

        DistSpec(Chi, (1.0,), 0.5),

        DistSpec(Chisq, (1.0,), 0.5),

        DistSpec(Cosine, (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(Exponential, (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(Gamma, (), 0.4),
        DistSpec(Gamma, (1.5,), 0.4),
        DistSpec(Gamma, (1.5, 2.0), 0.4),

        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(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(InverseGaussian, (), 0.5),
        DistSpec(InverseGaussian, (1.0,), 0.5),
        DistSpec(InverseGaussian, (1.0, 2.0), 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(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(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(LogNormal, (), 0.5),
        DistSpec(LogNormal, (1.0,), 0.5),
        DistSpec(LogNormal, (1.0, 2.0), 0.5),

        # Dispatch error caused by ccall
        DistSpec(NoncentralBeta, (1.0, 2.0, 1.0), 0.5, broken=(:Tracker, :ForwardDiff, :Zygote, :ReverseDiff)),
        DistSpec(NoncentralChisq, (1.0, 2.0), 0.5, broken=(:Tracker, :ForwardDiff, :Zygote, :ReverseDiff)),
        DistSpec(NoncentralF, (1.0, 2.0, 1.0), 0.5, broken=(:Tracker, :ForwardDiff, :Zygote, :ReverseDiff)),
        DistSpec(NoncentralT, (1.0, 2.0), 0.5, broken=(:Tracker, :ForwardDiff, :Zygote, :ReverseDiff)),

        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(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(PGeneralizedGaussian, (), 0.5),
        DistSpec(PGeneralizedGaussian, (1.0, 1.0, 1.0), 0.5),

        DistSpec(Rayleigh, (), 0.5),
        DistSpec(Rayleigh, (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(TDist, (1.0,), 0.5),

        DistSpec(TriangularDist, (1.0, 3.0), 1.5),
        DistSpec(TriangularDist, (1.0, 3.0, 2.0), 1.5),

        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(Uniform, (), 0.5),
        DistSpec(Uniform, (alpha, alpha + beta), alpha + beta * gamma),

        DistSpec(TuringUniform, (), 0.5),
        DistSpec(TuringUniform, (alpha, alpha + beta), alpha + beta * gamma),

        DistSpec(VonMises, (), 1.0),

        DistSpec(Weibull, (), 1.0),
        DistSpec(Weibull, (1.0,), 1.0),
        DistSpec(Weibull, (1.0, 1.0), 1.0),
    ]

    # Tests cannot be executed, so cannot be checked with `@test_broken`.
    broken_univariate_distributions = DistSpec[
        # Broken in Distributions even without autodiff
        DistSpec(() -> KSDist(1), (), 0.5), # `pdf` method not defined
        DistSpec(() -> KSOneSided(1), (), 0.5), # `pdf` method not defined
        DistSpec(StudentizedRange, (1.0, 2.0), 0.5), # `srdistlogpdf` method not defined

        # Stackoverflow caused by SpecialFunctions.besselix
        DistSpec(VonMises, (1.0,), 1.0),
        DistSpec(VonMises, (1, 1), 1),

        # Some tests are broken on some Julia versions, therefore it can't be checked reliably
        DistSpec(PoissonBinomial, ([0.5, 0.5],), [0, 0]; broken=(:Zygote,)),
    ]

    # Tests that have a `broken` field can be executed but, according to FiniteDifferences,
    # fail to produce the correct result. These tests can be checked with `@test_broken`.
    multivariate_distributions = DistSpec[
        ## Multivariate discrete distributions

        # Vector x
        DistSpec(p -> Multinomial(2, p ./ sum(p)), (fill(0.5, 2),), [2, 0]),
        DistSpec(p -> Multinomial(2, p ./ sum(p)), (fill(0.5, 2),), [2 1; 0 1]),

        # Vector x
        DistSpec((m, A) -> MvNormal(m, to_posdef(A)), (a, A), b),
        DistSpec((m, s) -> MvNormal(m, to_posdef_diagonal(s)), (a, b), c),
        DistSpec((m, s) -> MvNormal(m, s^2 * I), (a, alpha), b),
        DistSpec(A -> MvNormal(to_posdef(A)), (A,), a),
        DistSpec(s -> MvNormal(to_posdef_diagonal(s)), (a,), b),
        DistSpec(s -> MvNormal(zeros(dim), s^2 * I), (alpha,), a),
        DistSpec((m, A) -> TuringMvNormal(m, to_posdef(A)), (a, A), b),
        DistSpec((m, s) -> TuringMvNormal(m, to_posdef_diagonal(s)), (a, b), c),
        DistSpec((m, s) -> TuringMvNormal(m, s^2 * I), (a, alpha), b),
        DistSpec(A -> TuringMvNormal(to_posdef(A)), (A,), a),
        DistSpec(s -> TuringMvNormal(to_posdef_diagonal(s)), (a,), b),
        DistSpec(s -> TuringMvNormal(zeros(dim), s^2 * I), (alpha,), a),
        DistSpec((m, A) -> MvLogNormal(m, to_posdef(A)), (a, A), b, to_positive),
        DistSpec((m, s) -> MvLogNormal(m, to_posdef_diagonal(s)), (a, b), c, to_positive),
        DistSpec((m, s) -> MvLogNormal(m, s^2 * I), (a, alpha), b, to_positive),
        DistSpec(A -> MvLogNormal(to_posdef(A)), (A,), a, to_positive),
        DistSpec(s -> MvLogNormal(to_posdef_diagonal(s)), (a,), b, to_positive),
        DistSpec(s -> MvLogNormal(zeros(dim), s^2 * I), (alpha,), a, to_positive),

        DistSpec(alpha -> Dirichlet(to_positive(alpha)), (a,), b, to_simplex),

        # Matrix case
        DistSpec((m, A) -> MvNormal(m, to_posdef(A)), (a, A), B),
        DistSpec((m, s) -> MvNormal(m, to_posdef_diagonal(s)), (a, b), A),
        DistSpec((m, s) -> MvNormal(m, s^2 * I), (a, alpha), A),
        DistSpec(A -> MvNormal(to_posdef(A)), (A,), B),
        DistSpec(s -> MvNormal(to_posdef_diagonal(s)), (a,), A),
        DistSpec(s -> MvNormal(zeros(dim), s^2 * I), (alpha,), A),
        DistSpec((m, A) -> TuringMvNormal(m, to_posdef(A)), (a, A), B),
        DistSpec((m, s) -> TuringMvNormal(m, to_posdef_diagonal(s)), (a, b), A),
        DistSpec((m, s) -> TuringMvNormal(m, s^2 * I), (a, alpha), A),
        DistSpec(A -> TuringMvNormal(to_posdef(A)), (A,), B),
        DistSpec(s -> TuringMvNormal(to_posdef_diagonal(s)), (a,), A),
        DistSpec(s -> TuringMvNormal(zeros(dim), s^2 * I), (alpha,), A),
        DistSpec((m, A) -> MvLogNormal(m, to_posdef(A)), (a, A), B, to_positive),
        DistSpec((m, s) -> MvLogNormal(m, to_posdef_diagonal(s)), (a, b), A, to_positive),
        DistSpec((m, s) -> MvLogNormal(m, s^2 * I), (a, alpha), A, to_positive),
        DistSpec(A -> MvLogNormal(to_posdef(A)), (A,), B, to_positive),
        DistSpec(s -> MvLogNormal(to_posdef_diagonal(s)), (a,), A, to_positive),
        DistSpec(s -> MvLogNormal(zeros(dim), s^2 * I), (alpha,), A, to_positive),

        DistSpec(alpha -> Dirichlet(to_positive(alpha)), (a,), A, to_simplex),
    ]

    # Tests cannot be executed, so cannot be checked with `@test_broken`.
    broken_multivariate_distributions = DistSpec[
        # Dispatch error
        DistSpec((m, A) -> MvNormalCanon(m, to_posdef(A)), (a, A), b),
        DistSpec((m, p) -> MvNormalCanon(m, to_posdef_diagonal(p)), (a, b), c),
        DistSpec((m, p) -> MvNormalCanon(m, p^2 * I), (a, alpha), b),
        DistSpec(A -> MvNormalCanon(to_posdef(A)), (A,), a),
        DistSpec(p -> MvNormalCanon(to_posdef_diagonal(p)), (a,), b),
        DistSpec(p -> MvNormalCanon(zeros(dim), p^2 * I), (alpha,), a),
        DistSpec((m, A) -> MvNormalCanon(m, to_posdef(A)), (a, A), B),
        DistSpec((m, p) -> MvNormalCanon(m, to_posdef_diagonal(p)), (a, b), A),
        DistSpec((m, p) -> MvNormalCanon(m, p^2 * I), (a, alpha), A),
        DistSpec(A -> MvNormalCanon(to_posdef(A)), (A,), B),
        DistSpec(p -> MvNormalCanon(to_posdef_diagonal(p)), (a,), A),
        DistSpec(p -> MvNormalCanon(zeros(dim), p^2 * I), (alpha,), A),
    ]

    # Tests that have a `broken` field can be executed but, according to FiniteDifferences,
    # fail to produce the correct result. These tests can be checked with `@test_broken`.
    matrixvariate_distributions = DistSpec[
        # Matrix x
        # We should use
        # DistSpec((n1, n2) -> MatrixBeta(dim, n1, n2), (3.0, 3.0), A, to_beta_mat),
        # but the default implementation of `rand_tangent` causes a StackOverflowError
        # Thus we use the following workaround
        DistSpec((n1, n2) -> MatrixBeta(3, n1, n2), (3.0, 3.0), A, to_beta_mat),
        DistSpec(() -> MatrixNormal(dim, dim), (), A, to_posdef, broken=(:Zygote,)),
        DistSpec((df, A) -> Wishart(df, to_posdef(A)), (3.0, A), B, to_posdef),
        DistSpec((df, A) -> InverseWishart(df, to_posdef(A)), (3.0, A), B, to_posdef),
        DistSpec((df, A) -> TuringWishart(df, to_posdef(A)), (3.0, A), B, to_posdef),
        DistSpec((df, A) -> TuringInverseWishart(df, to_posdef(A)), (3.0, A), B, to_posdef),

        # Vector of matrices x
        # Also here we should use
        # DistSpec(
        #    (n1, n2) -> MatrixBeta(dim, n1, n2),
        #    (3.0, 3.0),
        #    [A, B],
        #    x -> map(to_beta_mat, x),
        #),
        # but the default implementation of `rand_tangent` causes a StackOverflowError
        # Thus we use the following workaround
        DistSpec(
            (n1, n2) -> MatrixBeta(3, n1, n2),
            (3.0, 3.0),
            [A, B],
            x -> map(to_beta_mat, x),
        ),
        DistSpec(
            (df, A) -> Wishart(df, to_posdef(A)),
            (3.0, A),
            [B, C],
            x -> map(to_posdef, x),
        ),
        DistSpec(
            (df, A) -> InverseWishart(df, to_posdef(A)),
            (3.0, A),
            [B, C],
            x -> map(to_posdef, x),
        ),
        DistSpec(
            (df, A) -> TuringWishart(df, to_posdef(A)),
            (3.0, A),
            [B, C],
            x -> map(to_posdef, x),
        ),
        DistSpec(
            (df, A) -> TuringInverseWishart(df, to_posdef(A)),
            (3.0, A),
            [B, C],
            x -> map(to_posdef, x),
        ),
    ]

    # Tests cannot be executed, so cannot be checked with `@test_broken`.
    broken_matrixvariate_distributions = DistSpec[
        # Other
        # TODO different tests are broken on different combinations of backends
        DistSpec(
            (A, B, C) -> MatrixNormal(A, to_posdef(B), to_posdef(C)),
            (A, B, B),
            C,
            to_posdef,
        ),
        # TODO different tests are broken on different combinations of backends
        DistSpec(
            (df, A, B, C) -> MatrixTDist(df, A, to_posdef(B), to_posdef(C)),
            (1.0, A, B, B),
            C,
            to_posdef,
        ),
        # TODO different tests are broken on different combinations of backends
        DistSpec(
            (n1, n2, A) -> MatrixFDist(n1, n2, to_posdef(A)),
            (3.0, 3.0, A),
            B,
            to_posdef,
        ),
    ]

    @testset "Univariate distributions" begin
        println("\nTesting: Univariate distributions\n")

        for d in univariate_distributions
            @info "Testing: $(nameof(dist_type(d)))"
            test_ad(d)
        end
    end

    @testset "Multivariate distributions" begin
        println("\nTesting: Multivariate distributions\n")

        for d in multivariate_distributions
            @info "Testing: $(nameof(dist_type(d)))"
            test_ad(d)
        end

        # Test `filldist` and `arraydist` distributions of univariate distributions
        n = 2 # always use two distributions
        for d in univariate_distributions
            d.x isa Number || continue

            # Broken distributions
            D = dist_type(d)
            D <: Union{VonMises,TriangularDist} && continue

            # Skellam only fails in these tests with ReverseDiff
            # Ref: https://github.com/TuringLang/DistributionsAD.jl/issues/126
            # PoissonBinomial fails with Zygote
            # Matrix case does not work with Skellam:
            # https://github.com/TuringLang/DistributionsAD.jl/pull/172#issuecomment-853721493
            filldist_broken = if D <: Skellam
                ((d.broken..., :Zygote, :ReverseDiff), (d.broken..., :Zygote, :ReverseDiff))
            elseif D <: PoissonBinomial
                ((d.broken..., :Zygote), (d.broken..., :Zygote))
            elseif D <: Chernoff
                # Zygote is not broken with `filldist`
                ((), ())
            else
                (d.broken, d.broken)
            end
            arraydist_broken = if D <: PoissonBinomial
                ((d.broken..., :Zygote), (d.broken..., :Zygote))
            else
                (d.broken, d.broken)
            end

            # Create `filldist` distribution
            f = d.f
            f_filldist = (θ...,) -> filldist(f(θ...), n)
            d_filldist = f_filldist(d.θ...)

            # Create `arraydist` distribution
            f_arraydist = (θ...,) -> arraydist([f(θ...) for _ in 1:n])
            d_arraydist = f_arraydist(d.θ...)

            for (i, sz) in enumerate(((n,), (n, 2)))
                # Matrix case doesn't work for continuous distributions for some reason
                # now but not too important (?!)
                if length(sz) == 2 && D <: ContinuousDistribution
                    continue
                end

                # Compute compatible sample
                x = fill(d.x, sz)

                # Test AD
                @info "Testing: filldist($(nameof(D)), $sz)"
                test_ad(
                    DistSpec(
                        f_filldist,
                        d.θ,
                        x,
                        d.xtrans;
                        broken=filldist_broken[i],
                    )
                )

                @info "Testing: arraydist($(nameof(D)), $sz)"
                test_ad(
                    DistSpec(
                        f_arraydist,
                        d.θ,
                        x,
                        d.xtrans;
                        broken=arraydist_broken[i],
                    )
                )
            end
        end
    end

    @testset "Matrixvariate distributions" begin
        println("\nTesting: Matrixvariate distributions\n")

        for d in matrixvariate_distributions
            @info "Testing: $(nameof(dist_type(d)))"
            test_ad(d)
        end

        # Test `filldist` and `arraydist` distributions of univariate distributions
        n = (2, 2) # always use 2 x 2 distributions
        for d in univariate_distributions
            d.x isa Number || continue
            D = dist_type(d)
            D <: DiscreteDistribution && continue

            # Broken distributions
            D <: Union{VonMises,TriangularDist} && continue

            # Create `filldist` distribution
            f = d.f
            f_filldist = (θ...,) -> filldist(f(θ...), n...)

            # Create `arraydist` distribution
            # Zygote's fill definition does not like non-numbers, so we use a workaround
            f_arraydist = (θ...,) -> arraydist(reshape([f(θ...) for _ in 1:prod(n)], n))

            # Matrix `x`
            x_mat = fill(d.x, n)

            # Zygote is not broken with `filldist` + Chernoff
            filldist_broken = D <: Chernoff ? () : d.broken

            # Test AD
            @info "Testing: filldist($(nameof(D)), $n)"
            test_ad(
                DistSpec(
                    f_filldist,
                    d.θ,
                    x_mat,
                    d.xtrans;
                    broken=filldist_broken,
                )
            )
            @info "Testing: arraydist($(nameof(D)), $n)"
            test_ad(
                DistSpec(
                    f_arraydist,
                    d.θ,
                    x_mat,
                    d.xtrans;
                    broken=d.broken,
                )
            )

            # Vector of matrices `x`
            x_vec_of_mat = [fill(d.x, n) for _ in 1:2]

            # Test AD
            @info "Testing: filldist($(nameof(D)), $n, 2)"
            test_ad(
                DistSpec(
                    f_filldist,
                    d.θ,
                    x_vec_of_mat,
                    d.xtrans;
                    broken=filldist_broken,
                )
            )
            @info "Testing: arraydist($(nameof(D)), $n, 2)"
            test_ad(
                DistSpec(
                    f_arraydist,
                    d.θ,
                    x_vec_of_mat,
                    d.xtrans;
                    broken=d.broken,
                )
            )
        end

        # test `filldist` and `arraydist` distributions of multivariate distributions
        n = 2 # always use two distributions
        for d in multivariate_distributions
            d.x isa AbstractVector || continue
            D = dist_type(d)
            D <: DiscreteDistribution && continue

            # Tests are failing for matrix covariance vectorized MvNormal
            if D <: Union{
                MvNormal,MvLogNormal,
                DistributionsAD.TuringDenseMvNormal,
                DistributionsAD.TuringDiagMvNormal,
                DistributionsAD.TuringScalMvNormal,
                TuringMvLogNormal
            }
                any(x isa Matrix for x in d.θ) && continue
            end

            # Create `filldist` distribution
            f = d.f
            f_filldist = (θ...,) -> filldist(f(θ...), n)

            # Create `arraydist` distribution
            f_arraydist = (θ...,) -> arraydist([f(θ...) for _ in 1:n])

            # Matrix `x`
            x_mat = repeat(d.x, 1, n)

            # Test AD
            @info "Testing: filldist($(nameof(D)), $n)"
            test_ad(
                DistSpec(
                    f_filldist,
                    d.θ,
                    x_mat,
                    d.xtrans;
                    broken=d.broken,
                )
            )
            @info "Testing: arraydist($(nameof(D)), $n)"
            test_ad(
                DistSpec(
                    f_arraydist,
                    d.θ,
                    x_mat,
                    d.xtrans;
                    broken=d.broken,
                )
            )

            # Vector of matrices `x`
            x_vec_of_mat = [repeat(d.x, 1, n) for _ in 1:2]

            # Test AD
            @info "Testing: filldist($(nameof(D)), $n, 2)"
            test_ad(
                DistSpec(
                    f_filldist,
                    d.θ,
                    x_vec_of_mat,
                    d.xtrans;
                    broken=d.broken,
                )
            )
            @info "Testing: arraydist($(nameof(D)), $n, 2)"
            test_ad(
                DistSpec(
                    f_arraydist,
                    d.θ,
                    x_vec_of_mat,
                    d.xtrans;
                    broken=d.broken,
                )
            )
        end
    end
end


text:
"""
    Affine_Exp!(x::Vector{SMCg{N,T}},p::Vector{SMCg{N,T}},p_ref::Vector{SMCg{N,T}},
               xa::Vector{SMCg{N,T}},xA::Vector{SMCg{N,T}},z::Vector{SMCg{N,T}},
               opt::Array{Any})

Computates the affine relaxations of the state variable. Inputs are:
* `x::Vector{SMCg{N,T}}`: State variable relaxation
* `p::Vector{SMCg{N,T}}`: Decision variable relaxation
* `p_ref::Vector{SMCg{N,T}}`: Reference variable relaxation
* `xa::Vector{SMCg{N,T}}`: Lower affine relaxation of the state variable
* `xA::Vector{SMCg{N,T}}`: Upper affine relaxation of the state variable
* `z::Vector{SMCg{N,T}}`: Affine function in `X`
* `opt::Array{Any,1}`: `[np,nx,lambda]` values for relaxation
Returns the tuple `(xa,xA,z)`:
* `xa::Vector{SMCg{N,T}}`: Lower affine relaxation of the state variable
* `xA::Vector{SMCg{N,T}}`: Upper affine relaxation of the state variable
* `z::Vector{SMCg{N,T}}`: Affine function in X
--------------------------------------------------------------------------------
"""
function Affine_Exp!(x::Vector{SMCg{N,T}}, p::Vector{SMCg{N,T}}, p_ref::Vector{SMCg{N,T}},
                     xa::Vector{SMCg{N,T}}, xA::Vector{SMCg{N,T}}, z::Vector{SMCg{N,T}},
                     opt::Vector{Any}) where {N,T<:AbstractFloat}

  nx::Int64,np::Int64,lambda::Float64 = opt
  S1::SMCg{N,T},S2::SMCg{N,T},S3::SMCg{N,T} = zero(SMCg{N,T}),zero(SMCg{N,T}),zero(SMCg{N,T})
  for i = 1:nx
   S1 = zero(SMCg{N,T})
   S2 = zero(SMCg{N,T})
   S3 = zero(SMCg{N,T})
   for j = 1:N
      S1 = S1 + (p[j]-p_ref[j])*x[i].cv_grad[j]
      S2 = S2 + (p[j]-p_ref[j])*x[i].cc_grad[j]
      S3 = S3 + (lambda*x[i].cv_grad[j]+(1.0-lambda)*x[i].cc_grad[j])*(p[j]-p_ref[j])
   end
   temp1::SMCg{N,T} = x[i].cv + S1
   temp2::SMCg{N,T} = x[i].cc + S2
   temp3::SMCg{N,T} = lambda*x[i].cv+(1.0-lambda)*x[i].cc+S3
   xa[i] = SMCg{N,T}(temp1.cc,temp1.cv,x[i].cv_grad,x[i].cv_grad,Interval{T}(temp1.cv,temp1.cc),S1.cnst,x[i].IntvBox,x[i].xref)
   xA[i] = SMCg{N,T}(temp2.cc,temp2.cv,x[i].cc_grad,x[i].cc_grad,Interval{T}(temp2.cv,temp2.cc),S2.cnst,x[i].IntvBox,x[i].xref)
   z[i] = SMCg{N,T}(temp3.cc,temp3.cv,
                    lambda*x[i].cv_grad+(1.0-lambda)*x[i].cc_grad,
                    lambda*x[i].cv_grad+(1.0-lambda)*x[i].cc_grad,
                    Interval{T}(temp3.cv,temp3.cc),S3.cnst,x[i].IntvBox,x[i].xref)
  end
end

"""
    Correct_Exp!(z_mc::Vector{SMCg{N,T}},x_mc::Vector{SMCg{N,T}},xL::Vector{Float64},
                 xU::Vector{Float64},nx::Int64,np::Int64,epsv::Float64)

Corrects the relaxation of the state variable `x_mc` if the affine relaxation,
'z_mc', exceeds the interval bounds `xL` or `xU`.
* `z_mc::Vector{SMCg{N,T}}`: Affine relaxation
* `x_mc::Vector{SMCg{N,T}}`: Relaxation of state variable
* `xL::Vector{Float64}`: Lower bound on state vector
* `xU::Vector{Float64}`: Upper bound on state vector
* `nx::Int64`: Size of the state vector
* `np::Int64`: Size of the decision vector
* `epsv::Float64`: Tolerance for checking that subgradient exceeds bound
"""
function Correct_Exp!(z_mc::Vector{SMCg{N,T}},x_mc::Vector{SMCg{N,T}},
                     X::Vector{Interval{T}},nx::Int64,np::Int64,
                     epsv::Float64) where {N,T<:AbstractFloat}
  zero_grad::SVector{N,T} = @SVector zeros(T,N)
  for i = 1:nx
    if (z_mc[i].Intv.lo-epsv < X[i].lo)
      x_mc[i] = SMCg{N,T}(x_mc[i].cc,X[i].lo,x_mc[i].cc_grad,zero_grad,x_mc[i].Intv,true,x_mc[i].IntvBox,x_mc[i].xref)
    end
    if (z_mc[i].Intv.hi+epsv > X[i].hi)
      x_mc[i] = SMCg(X[i].hi,x_mc[i].cv,zero_grad,x_mc[i].cv_grad,x_mc[i].Intv,true,x_mc[i].IntvBox,x_mc[i].xref)
    end
  end
end

ext: jl


text:
function FST_global(data::PopData)
    # observed/expected het per locus per pop
    het_df = @inbounds DataFrames.combine(
        groupby(data.loci, [:locus, :population]),
        :genotype => nonmissing => :n,
        :genotype => hetero_o => :het_obs,
        :genotype => (i -> hetero_e(i)) => :het_exp,
        :genotype => allele_freq => :alleles
    )

    # collapse down to retrieve averages and counts
    n_df = @inbounds DataFrames.combine(
        groupby(het_df, :locus),
        :n => (n -> sum(.!iszero.(n))) => :count,
        :n => (n -> sum(.!iszero.(n)) ./ sum(reciprocal.(n))) => :mn,
        [:het_obs, :het_exp, :n] => ((o, e, n) -> mean(skipmissing(gene_diversity_nei87.(e,o,sum(.!iszero.(n)) ./ sum(reciprocal.(n)))))) => :HS,
        :het_obs => (o -> mean(skipmissing(o))) => :Het_obs,
        :alleles => (alleles -> sum(values(avg_allele_freq(alleles)).^2)) => :avg_freq
    )

    HT = mean(@inbounds 1.0 .- n_df.avg_freq .+ (n_df.HS ./ n_df.mn ./ n_df.count) - (n_df.Het_obs ./ 2.0 ./ n_df.mn ./ n_df.count))
    DST = mean(@inbounds HT .- n_df.HS)
    DST′ = mean(@inbounds n_df.count ./ (n_df.count .- 1) .* DST)
    HT′ = mean(@inbounds n_df.HS .+ DST′)

    return Dict{Symbol, Float64}(
        :FST => DST / HT,
        :FST′ => DST′ / HT′
    )
end

function FST_locus(data::PopData)
    # observed/expected het per locus per pop
    het_df = DataFrames.combine(
        groupby(data.loci, [:locus, :population]),
        :genotype => nonmissing => :n,
        :genotype => hetero_o => :het_obs,
        :genotype => (i -> hetero_e(i)) => :het_exp,
        :genotype => allele_freq => :alleles
    )
    # collapse down to retrieve averages and counts
    n_df = DataFrames.combine(
        [:het_obs, :het_exp, :n, :alleles] => (o,e,n,alleles) -> (
            count = sum(.!iszero.(n)),
            mn = sum(.!iszero.(n)) ./ sum(reciprocal.(n)),
            HS = mean(skipmissing(gene_diversity_nei87.(e,o,sum(.!iszero.(n)) ./ sum(reciprocal.(n))))),
            Het_obs = mean(skipmissing(o)),
            avg_freq = sum(values(avg_allele_freq(alleles)).^2)
            ),
        groupby(het_df, :locus)
    )

    HT = 1.0 .- n_df.avg_freq .+ (n_df.HS ./ n_df.mn ./ n_df.count) - (n_df.Het_obs ./ 2.0 ./ n_df.mn ./ n_df.count)
    DST = HT .- n_df.HS
    DST′ = n_df.count ./ (n_df.count .- 1) .* DST
    HT′ = n_df.HS .+ DST′

    return Dict{Symbol, Vector{Float64}}(
        :FST => DST ./ HT,
        :FST′ => DST′ ./ HT′
        )
end

function FST_global(data::AbstractDataFrame)
    # observed/expected het per locus per pop
    het_df = @inbounds DataFrames.combine(
        groupby(data, [:locus, :population]),
        :genotype => nonmissing => :n,
        :genotype => hetero_o => :het_obs,
        :genotype => (i -> hetero_e(i)) => :het_exp,
        :genotype => allele_freq => :alleles
    )

    # collapse down to retrieve averages and counts
    n_df = @inbounds DataFrames.combine(
        groupby(het_df, :locus),
        :n => (n -> sum(.!iszero.(n))) => :count,
        :n => (n -> sum(.!iszero.(n)) ./ sum(reciprocal.(n))) => :mn,
        [:het_obs, :het_exp, :n] => ((o, e, n) -> mean(skipmissing(gene_diversity_nei87.(e,o,sum(.!iszero.(n)) ./ sum(reciprocal.(n)))))) => :HS,
        :het_obs => (o -> mean(skipmissing(o))) => :Het_obs,
        :alleles => (alleles -> sum(values(avg_allele_freq(alleles)).^2)) => :avg_freq
    )

    HT = mean(@inbounds 1.0 .- n_df.avg_freq + (n_df.HS ./ n_df.mn ./ n_df.count) - (n_df.Het_obs ./ 2.0 ./ n_df.mn ./ n_df.count))
    DST = mean(@inbounds HT .- n_df.HS)
    DST′ = mean(@inbounds n_df.count ./ (n_df.count .- 1) .* DST)
    HT′ = mean(@inbounds n_df.HS .+ DST′)

    return Dict{Symbol, Float64}(
        :FST => DST / HT,
        :FST′ => DST′ / HT′
    )
end

function f_stat_p(data::PopData; nperm::Int = 1000)
    observed = FST_global(data)
    popnames = data.meta.population
    perm_dict = Dict{Symbol,Vector{Float64}}(
        :FST => Vector{Float64}(undef, nperm-1),
        :FST′ => Vector{Float64}(undef, nperm-1)
        )
    p = Progress(nperm-1, dt = 1, color = :blue)
    @inbounds @sync for n in 1:nperm-1
        Base.Threads.@spawn begin
            tmp = FST_global(permute_samples!(copy(data.loci), popnames))
            perm_dict[:FST][n] = tmp[:FST]
            perm_dict[:FST′][n] = tmp[:FST′]
            next!(p)
        end
    end
    @inbounds for (k,v) in perm_dict
        observed[k] = (1 + sum(observed[k] .<= v))/nperm
    end

    return observed
end


text:
"""
Provides the [`doctest`](@ref) function that makes sure that the `jldoctest` code blocks
in the documents and docstrings run and are up to date.
"""
module DocTests
using DocStringExtensions

import ..Documenter:
    DocSystem,
    DocMeta,
    Documenter,
    Documents,
    Expanders,
    Utilities,
    IdDict

import Markdown, REPL
import .Utilities: Markdown2
import IOCapture

# Julia code block testing.
# -------------------------

mutable struct MutableMD2CodeBlock
    language :: String
    code :: String
end
MutableMD2CodeBlock(block :: Markdown2.CodeBlock) = MutableMD2CodeBlock(block.language, block.code)

struct DocTestContext
    file :: String
    doc :: Documents.Document
    meta :: Dict{Symbol, Any}
    DocTestContext(file::String, doc::Documents.Document) = new(file, doc, Dict())
end

"""
$(SIGNATURES)

Traverses the pages and modules in the documenter blueprint, searching and
executing doctests.

Will abort the document generation when an error is thrown. Use `doctest = false`
keyword in [`Documenter.makedocs`](@ref) to disable doctesting.
"""
function doctest(blueprint::Documents.DocumentBlueprint, doc::Documents.Document)
    @debug "Running doctests."
    # find all the doctest blocks in the pages
    for (src, page) in blueprint.pages
        doctest(page, doc)
    end

    # find all the doctest block in all the docstrings (within specified modules)
    for mod in blueprint.modules
        for (binding, multidoc) in DocSystem.getmeta(mod)
            for signature in multidoc.order
                doctest(multidoc.docs[signature], mod, doc)
            end
        end
    end
end

function doctest(page::Documents.Page, doc::Documents.Document)
    ctx = DocTestContext(page.source, doc) # FIXME
    ctx.meta[:CurrentFile] = page.source
    doctest(ctx, page.md2ast)
end

function doctest(docstr::Docs.DocStr, mod::Module, doc::Documents.Document)
    md = DocSystem.parsedoc(docstr)
    # Note: parsedocs / formatdoc in Base is weird. It double-wraps the docstring Markdown
    # in a Markdown.MD object..
    @assert isa(md, Markdown.MD) # relying on Julia internals here
    while length(md.content) == 1 && isa(first(md.content), Markdown.MD)
        md = first(md.content)
    end
    md2ast = try
        Markdown2.convert(Markdown2.MD, md)
    catch err
        @error """
            Markdown2 conversion error for a docstring in $(mod).
            This is a bug — please report this on the Documenter issue tracker
            """ docstr.data
        rethrow(err)
    end
    ctx = DocTestContext(docstr.data[:path], doc)
    merge!(ctx.meta, DocMeta.getdocmeta(mod))
    ctx.meta[:CurrentFile] = get(docstr.data, :path, nothing)
    doctest(ctx, md2ast)
end

function parse_metablock(ctx::DocTestContext, block::Markdown2.CodeBlock)
    @assert startswith(block.language, "@meta")
    meta = Dict{Symbol, Any}()
    for (ex, str) in Utilities.parseblock(block.code, ctx.doc, ctx.file)
        if Utilities.isassign(ex)
            try
                meta[ex.args[1]] = Core.eval(Main, ex.args[2])
            catch err
                push!(ctx.doc.internal.errors, :meta_block)
                @warn "Failed to evaluate `$(strip(str))` in `@meta` block." err
            end
        end
    end
    return meta
end

function doctest(ctx::DocTestContext, md2ast::Markdown2.MD)
    Markdown2.walk(md2ast) do node
        isa(node, Markdown2.CodeBlock) || return true
        if startswith(node.language, "jldoctest")
            doctest(ctx, node)
        elseif startswith(node.language, "@meta")
            merge!(ctx.meta, parse_metablock(ctx, node))
        else
            return true
        end
        return false
    end
end

function doctest(ctx::DocTestContext, block_immutable::Markdown2.CodeBlock)
    lang = block_immutable.language
    if startswith(lang, "jldoctest")
        # Define new module or reuse an old one from this page if we have a named doctest.
        name = match(r"jldoctest[ ]?(.*)$", split(lang, ';', limit = 2)[1])[1]
        sandbox = Utilities.get_sandbox_module!(ctx.meta, "doctest", name)

        # Normalise line endings.
        block = MutableMD2CodeBlock(block_immutable)
        block.code = replace(block.code, "\r\n" => "\n")

        # parse keyword arguments to doctest
        d = Dict()
        idx = findfirst(c -> c == ';', lang)
        if idx !== nothing
            kwargs = try
                Meta.parse("($(lang[nextind(lang, idx):end]),)")
            catch e
                e isa Meta.ParseError || rethrow(e)
                push!(ctx.doc.internal.errors, :doctest)
                file = ctx.meta[:CurrentFile]
                lines = Utilities.find_block_in_file(block.code, file)
                @warn("""
                    Unable to parse doctest keyword arguments in $(Utilities.locrepr(file, lines))
                    Use ```jldoctest name; key1 = value1, key2 = value2

                    ```$(lang)
                    $(block.code)
                    ```
                    """, parse_error = e)
                return false
            end
            for kwarg in kwargs.args
                if !(isa(kwarg, Expr) && kwarg.head === :(=) && isa(kwarg.args[1], Symbol))
                    push!(ctx.doc.internal.errors, :doctest)
                    file = ctx.meta[:CurrentFile]
                    lines = Utilities.find_block_in_file(block.code, file)
                    @warn("""
                        invalid syntax for doctest keyword arguments in $(Utilities.locrepr(file, lines))
                        Use ```jldoctest name; key1 = value1, key2 = value2

                        ```$(lang)
                        $(block.code)
                        ```
                        """)
                    return false
                end
                d[kwarg.args[1]] = Core.eval(sandbox, kwarg.args[2])
            end
        end
        ctx.meta[:LocalDocTestArguments] = d

        for expr in [get(ctx.meta, :DocTestSetup, []); get(ctx.meta[:LocalDocTestArguments], :setup, [])]
            Meta.isexpr(expr, :block) && (expr.head = :toplevel)
            try
                Core.eval(sandbox, expr)
            catch e
                push!(ctx.doc.internal.errors, :doctest)
                @error("could not evaluate expression from doctest setup.",
                    expression = expr, exception = e)
                return false
            end
        end
        if occursin(r"^julia> "m, block.code)
            eval_repl(block, sandbox, ctx.meta, ctx.doc, ctx.file)
        elseif occursin(r"^# output$"m, block.code)
            eval_script(block, sandbox, ctx.meta, ctx.doc, ctx.file)
        else
            push!(ctx.doc.internal.errors, :doctest)
            file = ctx.meta[:CurrentFile]
            lines = Utilities.find_block_in_file(block.code, file)
            @warn("""
                invalid doctest block in $(Utilities.locrepr(file, lines))
                Requires `julia> ` or `# output`

                ```$(lang)
                $(block.code)
                ```
                """)
        end
        delete!(ctx.meta, :LocalDocTestArguments)
    end
    false
end
doctest(ctx::DocTestContext, block) = true

# Doctest evaluation.

mutable struct Result
    block  :: MutableMD2CodeBlock # The entire code block that is being tested.
    input  :: String # Part of `block.code` representing the current input.
    output :: String # Part of `block.code` representing the current expected output.
    file   :: String # File in which the doctest is written. Either `.md` or `.jl`.
    value  :: Any        # The value returned when evaluating `input`.
    hide   :: Bool       # Semi-colon suppressing the output?
    stdout :: IOBuffer   # Redirected stdout/stderr gets sent here.
    bt     :: Vector     # Backtrace when an error is thrown.

    function Result(block, input, output, file)
        new(block, input, rstrip(output, '\n'), file, nothing, false, IOBuffer())
    end
end

function eval_repl(block, sandbox, meta::Dict, doc::Documents.Document, page)
    for (input, output) in repl_splitter(block.code)
        result = Result(block, input, output, meta[:CurrentFile])
        for (ex, str) in Utilities.parseblock(input, doc, page; keywords = false, raise=false)
            # Input containing a semi-colon gets suppressed in the final output.
            result.hide = REPL.ends_with_semicolon(str)
            if VERSION >= v"1.5.0-DEV.178"
                # Use the REPL softscope for REPL jldoctests,
                # see https://github.com/JuliaLang/julia/pull/33864
                ex = REPL.softscope!(ex)
            end
            c = IOCapture.capture(rethrow = InterruptException) do
                Core.eval(sandbox, ex)
            end
            Core.eval(sandbox, Expr(:global, Expr(:(=), :ans, QuoteNode(c.value))))
            result.value = c.value
            print(result.stdout, c.output)
            if c.error
                result.bt = c.backtrace
            end
            # don't evaluate further if there is a parse error
            isa(ex, Expr) && ex.head === :error && break
        end
        checkresult(sandbox, result, meta, doc)
    end
end

function eval_script(block, sandbox, meta::Dict, doc::Documents.Document, page)
    # TODO: decide whether to keep `# output` syntax for this. It's a bit ugly.
    #       Maybe use double blank lines, i.e.
    #
    #
    #       to mark `input`/`output` separation.
    input, output = split(block.code, r"^# output$"m, limit = 2)
    input  = rstrip(input, '\n')
    output = lstrip(output, '\n')
    result = Result(block, input, output, meta[:CurrentFile])
    for (ex, str) in Utilities.parseblock(input, doc, page; keywords = false, raise=false)
        c = IOCapture.capture(rethrow = InterruptException) do
            Core.eval(sandbox, ex)
        end
        result.value = c.value
        print(result.stdout, c.output)
        if c.error
            result.bt = c.backtrace
            break
        end
    end
    checkresult(sandbox, result, meta, doc)
end

function filter_doctests(strings::NTuple{2, AbstractString},
                         doc::Documents.Document, meta::Dict)
    meta_block_filters = get(meta, :DocTestFilters, [])
    meta_block_filters == nothing && meta_block_filters == []
    doctest_local_filters = get(meta[:LocalDocTestArguments], :filter, [])
    for r in [doc.user.doctestfilters; meta_block_filters; doctest_local_filters]
        if all(occursin.((r,), strings))
            strings = replace.(strings, (r => "",))
        end
    end
    return strings
end

# Regex used here to replace gensym'd module names could probably use improvements.
function checkresult(sandbox::Module, result::Result, meta::Dict, doc::Documents.Document)
    sandbox_name = nameof(sandbox)
    mod_regex = Regex("(Main\\.)?(Symbol\\(\"$(sandbox_name)\"\\)|$(sandbox_name))[,.]")
    mod_regex_nodot = Regex(("(Main\\.)?$(sandbox_name)"))
    outio = IOContext(result.stdout, :module => sandbox)
    if isdefined(result, :bt) # An error was thrown and we have a backtrace.
        # To avoid dealing with path/line number issues in backtraces we use `[...]` to
        # mark ignored output from an error message. Only the text prior to it is used to
        # test for doctest success/failure.
        head = replace(split(result.output, "\n[...]"; limit = 2)[1], mod_regex  => "")
        head = replace(head, mod_regex_nodot => "Main")
        str  = replace(error_to_string(outio, result.value, result.bt), mod_regex => "")
        str  = replace(str, mod_regex_nodot => "Main")

        str, head = filter_doctests((str, head), doc, meta)
        # Since checking for the prefix of an error won't catch the empty case we need
        # to check that manually with `isempty`.
        if isempty(head) || !startswith(str, head)
            if doc.user.doctest === :fix
                fix_doctest(result, str, doc)
            else
                report(result, str, doc)
                @debug "Doctest metadata" meta
                push!(doc.internal.errors, :doctest)
            end
        end
    else
        value = result.hide ? nothing : result.value # `;` hides output.
        output = replace(rstrip(sanitise(IOBuffer(result.output))), mod_regex => "")
        str = replace(result_to_string(outio, value), mod_regex => "")
        # Replace a standalone module name with `Main`.
        str = rstrip(replace(str, mod_regex_nodot => "Main"))
        filteredstr, filteredoutput = filter_doctests((str, output), doc, meta)
        if filteredstr != filteredoutput
            if doc.user.doctest === :fix
                fix_doctest(result, str, doc)
            else
                report(result, str, doc)
                @debug "Doctest metadata" meta
                push!(doc.internal.errors, :doctest)
            end
        end
    end
    return nothing
end

# Display doctesting results.

function result_to_string(buf, value)
    value === nothing || Base.invokelatest(show, IOContext(buf, :limit => true), MIME"text/plain"(), value)
    return sanitise(buf)
end

function error_to_string(buf, er, bt)
    # Remove unimportant backtrace info.
    bt = remove_common_backtrace(bt, backtrace())
    # Remove everything below the last eval call (which should be the one in IOCapture.capture)
    index = findlast(ptr -> Base.ip_matches_func(ptr, :eval), bt)
    bt = (index === nothing) ? bt : bt[1:(index - 1)]
    # Print a REPL-like error message.
    print(buf, "ERROR: ")
    Base.invokelatest(showerror, buf, er, bt)
    return sanitise(buf)
end

function remove_common_backtrace(bt, reference_bt)
    cutoff = nothing
    # We'll start from the top of the backtrace (end of the array) and go down, checking
    # if the backtraces agree
    for ridx in 1:length(bt)
        # Cancel search if we run out the reference BT or find a non-matching one frames:
        if ridx > length(reference_bt) || bt[length(bt) - ridx + 1] != reference_bt[length(reference_bt) - ridx + 1]
            cutoff = length(bt) - ridx + 1
            break
        end
    end
    # It's possible that the loop does not find anything, i.e. that all BT elements are in
    # the reference_BT too. In that case we'll just return an empty BT.
    bt[1:(cutoff === nothing ? 0 : cutoff)]
end

# Strip trailing whitespace from each line and return resulting string
function sanitise(buffer)
    out = IOBuffer()
    for line in eachline(seekstart(Base.unwrapcontext(buffer)[1]))
        println(out, rstrip(line))
    end
    return rstrip(String(take!(out)), '\n')
end

import .Utilities.TextDiff

function report(result::Result, str, doc::Documents.Document)
    diff = TextDiff.Diff{TextDiff.Words}(result.output, rstrip(str))
    lines = Utilities.find_block_in_file(result.block.code, result.file)
    @error("""
        doctest failure in $(Utilities.locrepr(result.file, lines))

        ```$(result.block.language)
        $(result.block.code)
        ```

        Subexpression:

        $(result.input)

        Evaluated output:

        $(rstrip(str))

        Expected output:

        $(result.output)

        """, diff)
end

function fix_doctest(result::Result, str, doc::Documents.Document)
    code = result.block.code
    filename = Base.find_source_file(result.file)
    # read the file containing the code block
    content = read(filename, String)
    # output stream
    io = IOBuffer(sizehint = sizeof(content))
    # first look for the entire code block
    # make a regex of the code that matches leading whitespace
    rcode = "(\\h*)" * replace(Utilities.regex_escape(code), "\\n" => "\\n\\h*")
    r = Regex(rcode)
    codeidx = findfirst(r, content)
    if codeidx === nothing
        @warn "could not find code block in source file"
        return
    end
    # use the capture group to identify indentation
    indent = match(r, content).captures[1]
    # write everything up until the code block
    write(io, content[1:prevind(content, first(codeidx))])
    # next look for the particular input string in the given code block
    # make a regex of the input that matches leading whitespace (for multiline input)
    rinput = "\\h*" * replace(Utilities.regex_escape(result.input), "\\n" => "\\n\\h*")
    r = Regex(rinput)
    inputidx = findfirst(r, code)
    if inputidx === nothing
        @warn "could not find input line in code block"
        return
    end
    # construct the new code-snippet (without indent)
    # first part: everything up until the last index of the input string
    newcode = code[1:last(inputidx)]
    isempty(result.output) && (newcode *= '\n') # issue #772
    # second part: the rest, with the old output replaced with the new one
    if result.output == ""
        # This works around a regression in Julia 1.5.0 (https://github.com/JuliaLang/julia/issues/36953)
        # Technically, it is only necessary if VERSION >= v"1.5.0-DEV.826"
        newcode *= str
        newcode *= code[nextind(code, last(inputidx)):end]
    else
        newcode *= replace(code[nextind(code, last(inputidx)):end], result.output => str, count = 1)
    end
    # replace internal code block with the non-indented new code, needed if we come back
    # looking to replace output in the same code block later
    result.block.code = newcode
    # write the new code snippet to the stream, with indent
    newcode = replace(newcode, r"^(.+)$"m => Base.SubstitutionString(indent * "\\1"))
    write(io, newcode)
    # write rest of the file
    write(io, content[nextind(content, last(codeidx)):end])
    # write to file
    write(filename, seekstart(io))
    return
end

# REPL doctest splitter.

const PROMPT_REGEX = r"^julia> (.*)$"
const SOURCE_REGEX = r"^       (.*)$"

function repl_splitter(code)
    lines  = split(string(code, "\n"), '\n')
    input  = String[]
    output = String[]
    buffer = IOBuffer() # temporary buffer for doctest inputs and outputs
    found_first_prompt = false
    while !isempty(lines)
        line = popfirst!(lines)
        prompt = match(PROMPT_REGEX, line)
        # We allow comments before the first julia> prompt
        !found_first_prompt && startswith(line, '#') && continue
        if prompt === nothing
            source = match(SOURCE_REGEX, line)
            if source === nothing
                savebuffer!(input, buffer)
                println(buffer, line)
                takeuntil!(PROMPT_REGEX, buffer, lines)
            else
                println(buffer, source[1])
            end
        else
            found_first_prompt = true
            savebuffer!(output, buffer)
            println(buffer, prompt[1])
        end
    end
    savebuffer!(output, buffer)
    zip(input, output)
end

function savebuffer!(out, buf)
    n = bytesavailable(seekstart(buf))
    n > 0 ? push!(out, rstrip(String(take!(buf)))) : out
end

function takeuntil!(r, buf, lines)
    while !isempty(lines)
        line = lines[1]
        if !occursin(r, line)
            println(buf, popfirst!(lines))
        else
            break
        end
    end
end

end


text:
module PlantGeom

using MultiScaleTreeGraph
import MeshViz: viz, viz!, Viz

# For 3D (OPF):
import Meshes
import GeometryBasics
import Makie
import Colors: RGBA, Colorant, RGB
import ColorSchemes: get, rainbow, colorschemes, ColorScheme

# Read OPF:
import EzXML: readxml, root, StreamReader, attributes, expandtree # For reading OPF files
import EzXML: eachelement, nodename, nodecontent, elements
import EzXML: XMLDocument, ElementNode, setroot!, addelement!
import EzXML: prettyprint # to remove
import StaticArrays: SMatrix, SVector
import LinearAlgebra: I, UniformScaling # Used for geometry parsing in OPF
import CoordinateTransformations: Transformation, Translation, LinearMap, ∘, IdentityTransformation
import RecipesBase
import Base

include("structs.jl")
include("equality.jl")
include("helpers.jl")
include("opf/read_opf.jl")
include("opf/reference_meshes.jl")
include("tapering.jl")
include("opf/refmesh_to_mesh.jl")
include("opf/write_opf.jl")
include("colors.jl")
include("makie_recipes/RefMeshes_recipes.jl")
include("makie_recipes/opf_recipe.jl")
include("makie_recipes/mtg_recipe_helpers.jl")
include("makie_recipes/mtg_tree_recipe.jl")
include("meshes/summary_coordinates.jl")
include("meshes/transformations.jl")
include("ref_meshes/cylinder_refmesh.jl")
include("plots_recipes/plots_recipe.jl")
include("makie_recipes/colorbar.jl")

# 3D Plotting (OPF):
export get_ref_meshes
export viz, viz!
export nvertices, nelements
export read_opf, write_opf
export taper
export refmesh_to_mesh, get_ref_meshes_color
export xmax, ymax, zmax, xmin, ymin, zmin
export refmesh_to_mesh!
export transform_mesh!
export Material, Phong
export RefMesh, RefMeshes
export (==), names
export Diagram, diagram, diagram!
export cylinder
export colorbar

end


text:
using Test
using QuantumOpticsBase

@testset "subspace" begin

b = FockBasis(3)

u = Ket[fockstate(b, 1), fockstate(b, 2)]
v = Ket[fockstate(b, 2), fockstate(b, 1)]

bu = SubspaceBasis(u)
bv = SubspaceBasis(v)

T1 = projector(bu, b)
T2 = projector(bv, b)
T12 = projector(bu, bv)

state = fockstate(b, 2)
state_u = Ket(bu, [0, 1])
state_v = Ket(bv, [1., 0])

@test T1*state == state_u
@test T2*state == state_v

state_v = Ket(bv, [1, -1])
state_u = Ket(bu, [-1, 1])

@test T12*state_v == state_u

u2 = Ket[1.5*fockstate(b, 1), fockstate(b, 1) + fockstate(b, 2)]
bu2_orth = QuantumOpticsBase.orthonormalize(SubspaceBasis(u))
@test bu2_orth.basisstates == bu.basisstates

# Test sparse version
T1 = sparseprojector(bu, b)
T2 = sparseprojector(bv, b)
T12 = sparseprojector(bu, bv)

@test isa(T1, SparseOpType)

state_u = Ket(bu, [0, 1])
state_v = Ket(bv, [1., 0])

@test T1*state == state_u
@test T2*state == state_v

state_v = Ket(bv, [1, -1])
state_u = Ket(bu, [-1, 1])

@test T12*state_v == state_u

# Test errors
b2 = FockBasis(4)
@test_throws ArgumentError SubspaceBasis(b2, u)
@test_throws ArgumentError projector(bu, b2)
@test_throws ArgumentError projector(b2, bu)
b2_sub = SubspaceBasis(b2, [fockstate(b2, 1)])
@test_throws ArgumentError projector(bu, b2_sub)

end # testset


text:
using Printf, ProgressMeter
@everywhere using MDToolbox

using BSON: @save, @load
@load "afm_generate_test_data.bson" afms models models_template qs myseed sigma_noise nframe

configs = [AfmizeConfig(10.0 * (pi / 180),
    r,
    MDToolbox.Point2D(-250, -200),
    MDToolbox.Point2D(250, 200),
    MDToolbox.Point2D(6.25, 6.25),
    MDToolbox.defaultParameters())
    for r in [10.0, 15.0, 20.0, 25.0, 30.0]]

# single frame
#r = getafmposterior(afms[1], models_template, qs, configs)

# all frames with parallel map
@time r = @showprogress pmap(x -> getafmposterior(x, models_template, qs, configs), afms[1:20])

@save "afm_analyze_test_data_parallel.bson" r configs

max_stars_repo_name: matsunagalab/notebooks

max_issues_repo_name: matsunagalab/notebooks

max_forks_repo_name: matsunagalab/notebooks


text:
module Applicative

export Callable
export apply, <|
export fimap, fmap, iflip, trpos, make_indices

using SimpleTraits
using IterTools: imap

@traitdef Callable{F}

@traitimpl Callable{Function}
@traitimpl Callable{Type}

apply(f, x) = f(x)
<|(f, x) = f(x)

fimap(f) = a -> imap(f, a)
fmap(f) = a -> map(f, a)

iflip(x, y) = (y, x)
iflip((x, y)) = (y, x)
iflip(i::CartesianIndex{2}) = Tuple(i) |> iflip |> CartesianIndex

trpos(t) = p -> CartesianIndex(Tuple(p) .+ t)

function make_indices(rows, cols, itr)
    (
        fimap(trpos((rows ÷ 2 + 1, cols ÷ 2 + 1))) ∘ fimap(iflip) ∘
        fimap(x -> CartesianIndex(round.(Int, Tuple(x))))
    )(
        itr,
    )
end

end # module


text:
using Documenter, MinimalPerfectHash

makedocs(;
    modules=[MinimalPerfectHash],
    doctest=true,
    sitename="MinimalPerfectHash.jl",
    format=Documenter.HTML(),
    authors="Kyungmin Lee",
    checkdocs=:all,
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/kyungminlee/MinimalPerfectHash.jl",
    devbranch="dev",
)

max_stars_repo_name: kyungminlee/MinimalPerfectHashing

max_issues_repo_name: kyungminlee/MinimalPerfectHashing

max_forks_repo_name: kyungminlee/MinimalPerfectHashing


text:
"""
Any StaticInjection struct that wants to support dynamic injectors must implement this
method to set the value.

The method is only for internal uses.
"""
function set_dynamic_injector!(
    static_injector::T,
    dynamic_injector::U,
) where {T <: StaticInjection, U <: Union{Nothing, DynamicInjection}}
    current_dynamic_injector = get_dynamic_injector(static_injector)
    if !isnothing(current_dynamic_injector) && !isnothing(dynamic_injector)
        throw(ArgumentError("cannot assign a dynamic injector on a device that already has one"))
    end

    # All of these types implement this field.
    static_injector.dynamic_injector = dynamic_injector
end