text:
# TODO: uncomment vexity checks once SDP on vars/constraints changes vexity of problem
@add_problem sdp function sdp_sdp_variables(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    y = Variable((2, 2), :Semidefinite)
    p = minimize(y[1, 1]; numeric_type = T)

    # @fact vexity(p) --> ConvexVexity()
    handle_problem!(p)
    if test
        @test p.optval ≈ 0 atol = atol rtol = rtol
    end

    y = Variable((3, 3), :Semidefinite)
    p = minimize(y[1, 1], y[2, 2] == 1; numeric_type = T)

    # @fact vexity(p) --> ConvexVexity()
    handle_problem!(p)
    if test
        @test p.optval ≈ 0 atol = atol rtol = rtol
    end

    # Solution is obtained as y[2,2] -> infinity
    # This test fails on Mosek. See
    # https://github.com/JuliaOpt/Mosek.jl/issues/29
    # y = Variable((2, 2), :Semidefinite)
    # p = minimize(y[1, 1], y[1, 2] == 1; numeric_type = T)

    # # @fact vexity(p) --> ConvexVexity()
    # handle_problem!(p)
    # @fact p.optval --> roughly(0, atol)

    y = Semidefinite(3)
    p = minimize(sum(diag(y)), y[1, 1] == 1; numeric_type = T)

    # @fact vexity(p) --> ConvexVexity()
    handle_problem!(p)
    if test
        @test p.optval ≈ 1 atol = atol rtol = rtol
    end

    y = Variable((3, 3), :Semidefinite)
    p = minimize(tr(y), y[2, 1] <= 4, y[2, 2] >= 3; numeric_type = T)

    # @fact vexity(p) --> ConvexVexity()
    handle_problem!(p)
    if test
        @test p.optval ≈ 3 atol = atol rtol = rtol
    end

    x = Variable(Positive())
    y = Semidefinite(3)
    p = minimize(y[1, 2], y[2, 1] == 1; numeric_type = T)

    # @fact vexity(p) --> ConvexVexity()
    handle_problem!(p)
    if test
        @test p.optval ≈ 1 atol = atol rtol = rtol
    end
end

@add_problem sdp function sdp_sdp_constraints(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    x = Variable(Positive())
    y = Variable((3, 3))
    p = minimize(
        x + y[1, 1],
        isposdef(y),
        x >= 1,
        y[2, 1] == 1;
        numeric_type = T,
    )

    # @fact vexity(p) --> ConvexVexity()
    handle_problem!(p)
    if test
        @test p.optval ≈ 1 atol = atol rtol = rtol
    end
end

@add_problem sdp function sdp_nuclear_norm_atom(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    y = Semidefinite(3)
    p = minimize(
        nuclearnorm(y),
        y[2, 1] <= 4,
        y[2, 2] >= 3,
        y[3, 3] <= 2;
        numeric_type = T,
    )

    if test
        @test vexity(p) == ConvexVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 3 atol = atol rtol = rtol
        @test evaluate(nuclearnorm(y)) ≈ 3 atol = atol rtol = rtol
    end
end

@add_problem sdp function sdp_nuclear_norm_atom_complex(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    A = [1 2im 3 4; 4im 3im 2 1; 4 5 6 7]
    y = ComplexVariable(3, 4)
    p = minimize(nuclearnorm(y), y == A; numeric_type = T)

    if test
        @test vexity(p) == ConvexVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ sum(svdvals(A)) atol = atol rtol = rtol
        @test evaluate(nuclearnorm(y)) ≈ sum(svdvals(A)) atol = atol rtol = rtol
    end
end

@add_problem sdp function sdp_operator_norm_atom(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    y = Variable((3, 3))
    p = minimize(
        opnorm(y),
        y[2, 1] <= 4,
        y[2, 2] >= 3,
        sum(y) >= 12;
        numeric_type = T,
    )

    if test
        @test vexity(p) == ConvexVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 4 atol = atol rtol = rtol
        @test evaluate(opnorm(y)) ≈ 4 atol = atol rtol = rtol
    end
end

@add_problem sdp function sdp_operator_norm_atom_complex(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    A = [1 2im 3 4; 4im 3im 2 1; 4 5 6 7]
    y = ComplexVariable(3, 4)
    p = minimize(opnorm(y), y == A; numeric_type = T)

    if test
        @test vexity(p) == ConvexVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ maximum(svdvals(A)) atol = atol rtol = rtol
        @test evaluate(opnorm(y)) ≈ maximum(svdvals(A)) atol = atol rtol = rtol
    end
end

@add_problem sdp function sdp_sigma_max_atom(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    y = Variable((3, 3))
    p = minimize(
        sigmamax(y),
        y[2, 1] <= 4,
        y[2, 2] >= 3,
        sum(y) >= 12;
        numeric_type = T,
    )

    if test
        @test vexity(p) == ConvexVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 4 atol = atol rtol = rtol
        @test evaluate(sigmamax(y)) ≈ 4 atol = atol rtol = rtol
    end
end

@add_problem sdp function sdp_dual_lambda_max_atom(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    y = Semidefinite(3)
    p = minimize(eigmax(y), y[1, 1] >= 4; numeric_type = T)

    if test
        @test vexity(p) == ConvexVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 4 atol = atol rtol = rtol
        @test evaluate(eigmax(y)) ≈ 4 atol = atol rtol = rtol
    end

    # https://github.com/jump-dev/Convex.jl/issues/337
    x = ComplexVariable(2, 2)
    p = minimize(
        eigmax(x),
        [x[1, 2] == im, x[2, 2] == 1.0, x ⪰ -eye(2)];
        numeric_type = T,
    )
    handle_problem!(p)
    if test
        @test p.optval ≈ 1.5 atol = atol rtol = rtol
        @test p.constraints[1].dual ≈ im atol = atol rtol = rtol
        @test p.constraints[2].dual ≈ 0.75 atol = atol rtol = rtol
    end
end

@add_problem sdp function sdp_lambda_min_atom(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    y = Semidefinite(3)
    p = maximize(eigmin(y), tr(y) <= 6; numeric_type = T)

    if test
        @test vexity(p) == ConvexVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 2 atol = atol rtol = rtol
        @test evaluate(eigmin(y)) ≈ 2 atol = atol rtol = rtol
    end
end

@add_problem sdp function sdp_matrix_frac_atom(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    x = [1, 2, 3]
    P = Variable(3, 3)
    p = minimize(
        matrixfrac(x, P),
        P <= 2 * eye(3),
        P >= 0.5 * eye(3);
        numeric_type = T,
    )

    if test
        @test vexity(p) == ConvexVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 7 atol = atol rtol = rtol
        @test (evaluate(matrixfrac(x, P)))[1] ≈ 7 atol = atol rtol = rtol
    end
end

@add_problem sdp function sdp_matrix_frac_atom_both_arguments_variable(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    x = Variable(3)
    P = Variable(3, 3)
    p = minimize(matrixfrac(x, P), eigmax(P) <= 2, x[1] >= 1; numeric_type = T)

    if test
        @test vexity(p) == ConvexVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 0.5 atol = atol rtol = rtol
        @test (evaluate(matrixfrac(x, P)))[1] ≈ 0.5 atol = atol rtol = rtol
    end
end

@add_problem sdp function sdp_sum_largest_eigs(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    x = Semidefinite(3)
    p = minimize(sumlargesteigs(x, 2), x >= 1; numeric_type = T)

    handle_problem!(p)

    if test
        @test p.optval ≈ 3 atol = atol rtol = rtol
        @test evaluate(x) ≈ ones(3, 3) atol = atol rtol = rtol
    end

    x = Semidefinite(3)
    p = minimize(
        sumlargesteigs(x, 2) + sumlargesteigs(x, 0),
        [x[i, :] >= i for i in 1:3]...;
        numeric_type = T,
    )

    handle_problem!(p)

    if test
        @test p.optval ≈ 8.4853 atol = atol rtol = rtol
    end

    A = [
        1 -2im 3 4
        2im 7 3im 5
        3 -3im 2 9
        4 5 9 4
    ]

    x = ComplexVariable(4, 4)
    p = minimize(sumlargesteigs(x, 3), x == A; numeric_type = T)

    handle_problem!(p)

    if test
        @test p.optval ≈ sum(eigvals(A)[2:end]) atol = atol rtol = rtol
    end

    x1 = Semidefinite(3)
    p1 = minimize(eigmax(x1), x1[1, 1] >= 4; numeric_type = T)

    handle_problem!(p1)

    x2 = Semidefinite(3)
    p2 = minimize(sumlargesteigs(x2, 1), x2[1, 1] >= 4; numeric_type = T)

    handle_problem!(p2)

    if test
        @test p1.optval ≈ p2.optval atol = atol rtol = rtol
    end

    x1 = Semidefinite(3)
    p1 = minimize(eigmax(x1), [x1[i, :] >= i for i in 1:3]...; numeric_type = T)

    handle_problem!(p1)

    x2 = Semidefinite(3)
    p2 = minimize(
        sumlargesteigs(x2, 1),
        [x2[i, :] >= i for i in 1:3]...;
        numeric_type = T,
    )

    handle_problem!(p2)

    if test
        @test p1.optval ≈ p2.optval atol = atol rtol = rtol
    end
end

@add_problem sdp function sdp_kron_atom(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    id = eye(4)
    X = Semidefinite(4)
    W = kron(id, X)
    p = maximize(tr(W), tr(X) ≤ 1; numeric_type = T)

    if test
        @test vexity(p) == AffineVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 4 atol = atol rtol = rtol
    end
end

@add_problem sdp function sdp_Partial_trace(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    A = Semidefinite(2)
    B = [1 0; 0 0]
    ρ = kron(B, A)
    constraints = [
        partialtrace(ρ, 1, [2; 2]) ==
        [0.09942819 0.29923607; 0.29923607 0.90057181],
        ρ in :SDP,
    ]
    p = satisfy(constraints; numeric_type = T)

    handle_problem!(p)
    if test
        @test evaluate(ρ) ≈ [
            0.09942819 0.29923607 0 0
            0.299237 0.900572 0 0
            0 0 0 0
            0 0 0 0
        ] atol = atol rtol = rtol
        @test evaluate(partialtrace(ρ, 1, [2; 2])) ≈
              [0.09942819 0.29923607; 0.29923607 0.90057181] atol = atol rtol =
            rtol
    end

    function rand_normalized(n)
        A = 5 * randn(n, n) + im * 5 * randn(n, n)
        return A / tr(A)
    end

    As = [rand_normalized(3) for _ in 1:5]
    Bs = [rand_normalized(2) for _ in 1:5]
    p = rand(5)

    AB = sum(i -> p[i] * kron(As[i], Bs[i]), 1:5)
    if test
        @test partialtrace(AB, 2, [3, 2]) ≈ sum(p .* As) atol = atol rtol = rtol
        @test partialtrace(AB, 1, [3, 2]) ≈ sum(p .* Bs) atol = atol rtol = rtol
    end

    A, B, C = rand(5, 5), rand(4, 4), rand(3, 3)
    ABC = kron(kron(A, B), C)
    if test
        @test kron(A, B) * tr(C) ≈ partialtrace(ABC, 3, [5, 4, 3]) atol = atol rtol =
            rtol
    end

    # Test 281
    A = rand(6, 6)
    expr = partialtrace(Constant(A), 1, [2, 3])
    if test
        @test size(expr) == size(evaluate(expr))

        @test_throws ArgumentError partialtrace(rand(6, 6), 3, [2, 3])
        @test_throws ArgumentError partialtrace(rand(6, 6), 1, [2, 4])
        @test_throws ArgumentError partialtrace(rand(3, 4), 1, [2, 3])
    end
end

@add_problem sdp function sdp_Real_Variables_with_complex_equality_constraints(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    n = 10 # variable dimension (parameter)
    m = 5 # number of constraints (parameter)
    xo = rand(n)
    A = randn(m, n) + im * randn(m, n)
    b = A * xo
    x = Variable(n)
    p1 = minimize(sum(x), A * x == b, x >= 0; numeric_type = T)

    handle_problem!(p1)

    if test
        x1 = copy(evaluate(x))
    end

    p2 = minimize(
        sum(x),
        real(A) * x == real(b),
        imag(A) * x == imag(b),
        x >= 0;
        numeric_type = T,
    )

    handle_problem!(p2)
    if test
        x2 = evaluate(x)
        @test x1 == x2
    end
end

@add_problem sdp function sdp_Complex_Variable_with_complex_equality_constraints(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    n = 10 # variable dimension (parameter)
    m = 5 # number of constraints (parameter)
    xo = rand(n) + im * rand(n)
    A = randn(m, n) + im * randn(m, n)
    b = A * xo
    x = ComplexVariable(n)
    p1 = minimize(
        real(sum(x)),
        A * x == b,
        real(x) >= 0,
        imag(x) >= 0;
        numeric_type = T,
    )

    handle_problem!(p1)

    if test
        x1 = evaluate(x)
    end

    xr = Variable(n)
    xi = Variable(n)
    p2 = minimize(
        sum(xr),
        real(A) * xr - imag(A) * xi == real(b),
        imag(A) * xr + real(A) * xi == imag(b),
        xr >= 0,
        xi >= 0;
        numeric_type = T,
    )

    handle_problem!(p2)

    if test
        real_diff = real(x1) - evaluate(xr)
        @test real_diff ≈ zeros(10, 1) atol = atol rtol = rtol

        imag_diff = imag(x1) - evaluate(xi)
        @test imag_diff ≈ zeros(10, 1) atol = atol rtol = rtol
    end
end

@add_problem sdp function sdp_Issue_198(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    ρ = HermitianSemidefinite(2)
    constraints = [ρ == [1.0 0.0; 0.0 1.0]]
    p = satisfy(constraints; numeric_type = T)

    handle_problem!(p)
    if test
        @test p.status == MOI.OPTIMAL
        @test evaluate(ρ) ≈ [1.0 0.0; 0.0 1.0] atol = atol rtol = rtol
        @test p.optval ≈ 0 atol = atol rtol = rtol
    end
end

@add_problem sdp function sdp_socp_norm2_atom(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    a = 2 + 4im
    x = ComplexVariable()
    objective = norm2(a - x)
    c1 = real(x) >= 0
    p = minimize(objective, c1; numeric_type = T)

    handle_problem!(p)
    if test
        @test p.optval ≈ 0 atol = atol rtol = rtol
        @test evaluate(objective) ≈ 0 atol = atol rtol = rtol

        real_diff = real(evaluate(x)) - real(a)
        imag_diff = imag(evaluate(x)) - imag(a)
        @test real_diff ≈ 0 atol = atol rtol = rtol
        @test imag_diff ≈ 0 atol = atol rtol = rtol
    end
end

@add_problem sdp function sdp_socp_sumsquares_atom(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    a = [2 + 4im; 4 + 6im]
    x = ComplexVariable(2)
    objective = sumsquares(a - x)
    c1 = real(x) >= 0
    p = minimize(objective, c1; numeric_type = T)

    handle_problem!(p)
    if test
        @test p.optval ≈ 0 atol = atol rtol = rtol
        @test evaluate(objective) ≈ 0.0 atol = atol rtol = rtol

        real_diff = real.(evaluate(x)) - real.(a)
        imag_diff = imag.(evaluate(x)) - imag.(a)
        @test real_diff ≈ zeros(2, 1) atol = atol rtol = rtol
        @test imag_diff ≈ zeros(2, 1) atol = atol rtol = rtol
    end
end

@add_problem sdp function sdp_socp_abs_atom(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    a = 5 - 4im
    x = ComplexVariable()
    objective = abs(a - x)
    c1 = real(x) >= 0
    p = minimize(objective, c1; numeric_type = T)

    handle_problem!(p)
    if test
        @test p.optval ≈ 0 atol = atol rtol = rtol
        @test evaluate(objective) ≈ 0.0 atol = atol rtol = rtol

        real_diff = real(evaluate(x)) .- real(a)
        imag_diff = imag(evaluate(x)) .- imag(a)
        @test real_diff ≈ 0.0 atol = atol rtol = rtol
        @test imag_diff ≈ 0.0 atol = atol rtol = rtol
    end
end

@add_problem sdp function sdp_Complex_Semidefinite_constraint(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    n = 10
    A = rand(n, n) + im * rand(n, n)
    A = A + A' # now A is hermitian
    x = ComplexVariable(n, n)
    objective = sumsquares(A - x)
    c1 = x in :SDP
    p = minimize(objective, c1; numeric_type = T)

    handle_problem!(p)
    # test that X is approximately equal to posA:
    l, v = eigen(A)
    posA = v * Diagonal(max.(l, 0)) * v'

    if test
        real_diff = real.(evaluate(x)) - real.(posA)
        imag_diff = imag.(evaluate(x)) - imag.(posA)
        @test real_diff ≈ zeros(n, n) atol = atol rtol = rtol
        @test imag_diff ≈ zeros(n, n) atol = atol rtol = rtol
    end
end

@add_problem sdp function sdp_geom_mean_hypocone_argcheck(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    if test
        @test_throws DimensionMismatch GeomMeanHypoCone(
            zeros(2, 3),
            zeros(2, 3),
            1 // 2,
        )
        @test_throws DimensionMismatch GeomMeanHypoCone(
            zeros(2, 3),
            Variable(2, 3),
            1 // 2,
        )
        @test_throws DimensionMismatch GeomMeanHypoCone(
            Variable(2, 3),
            zeros(2, 3),
            1 // 2,
        )
        @test_throws DimensionMismatch GeomMeanHypoCone(
            Variable(2, 3),
            Variable(2, 3),
            1 // 2,
        )

        @test_throws DimensionMismatch GeomMeanHypoCone(
            zeros(2, 2),
            zeros(3, 3),
            1 // 2,
        )
        @test_throws DimensionMismatch GeomMeanHypoCone(
            zeros(2, 2),
            Variable(3, 3),
            1 // 2,
        )
        @test_throws DimensionMismatch GeomMeanHypoCone(
            Variable(2, 2),
            zeros(3, 3),
            1 // 2,
        )
        @test_throws DimensionMismatch GeomMeanHypoCone(
            Variable(2, 2),
            Variable(2, 3),
            1 // 2,
        )

        @test_throws DimensionMismatch Variable(2, 2) in GeomMeanHypoCone(
            Variable(3, 3),
            Variable(3, 3),
            1 // 2,
        )

        @test_throws DomainError GeomMeanHypoCone(
            Variable(3, 3),
            Variable(3, 3),
            -1 // 2,
        )
        @test_throws DomainError GeomMeanHypoCone(
            Variable(3, 3),
            Variable(3, 3),
            3 // 2,
        )
    end
end

@add_problem sdp function sdp_geom_mean_epicone_argcheck(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    if test
        @test_throws DimensionMismatch GeomMeanEpiCone(
            zeros(2, 3),
            zeros(2, 3),
            -1 // 2,
        )
        @test_throws DimensionMismatch GeomMeanEpiCone(
            zeros(2, 3),
            Variable(2, 3),
            -1 // 2,
        )
        @test_throws DimensionMismatch GeomMeanEpiCone(
            Variable(2, 3),
            zeros(2, 3),
            -1 // 2,
        )
        @test_throws DimensionMismatch GeomMeanEpiCone(
            Variable(2, 3),
            Variable(2, 3),
            -1 // 2,
        )

        @test_throws DimensionMismatch GeomMeanEpiCone(
            zeros(2, 2),
            zeros(3, 3),
            -1 // 2,
        )
        @test_throws DimensionMismatch GeomMeanEpiCone(
            zeros(2, 2),
            Variable(3, 3),
            -1 // 2,
        )
        @test_throws DimensionMismatch GeomMeanEpiCone(
            Variable(2, 2),
            zeros(3, 3),
            -1 // 2,
        )
        @test_throws DimensionMismatch GeomMeanEpiCone(
            Variable(2, 2),
            Variable(2, 3),
            -1 // 2,
        )

        @test_throws DimensionMismatch Variable(2, 2) in GeomMeanEpiCone(
            Variable(3, 3),
            Variable(3, 3),
            -1 // 2,
        )
        @test_throws DomainError GeomMeanEpiCone(
            Variable(3, 3),
            Variable(3, 3),
            -3 // 2,
        )
        @test_throws DomainError GeomMeanEpiCone(
            Variable(3, 3),
            Variable(3, 3),
            1 // 2,
        )
        @test_throws DomainError GeomMeanEpiCone(
            Variable(3, 3),
            Variable(3, 3),
            5 // 2,
        )
    end
end

@add_problem sdp function sdp_relative_entropy_argcheck(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    if test
        @test_throws DimensionMismatch RelativeEntropyEpiCone(
            zeros(2, 3),
            zeros(2, 3),
            3,
            3,
            eye(2),
        )
        @test_throws DimensionMismatch RelativeEntropyEpiCone(
            zeros(2, 3),
            Variable(2, 3),
            3,
            3,
            eye(2),
        )
        @test_throws DimensionMismatch RelativeEntropyEpiCone(
            Variable(2, 3),
            zeros(2, 3),
            3,
            3,
            eye(2),
        )
        @test_throws DimensionMismatch RelativeEntropyEpiCone(
            Variable(2, 3),
            Variable(2, 3),
            3,
            3,
            eye(2),
        )

        @test_throws DimensionMismatch RelativeEntropyEpiCone(
            zeros(2, 2),
            zeros(3, 3),
            3,
            3,
        )
        @test_throws DimensionMismatch RelativeEntropyEpiCone(
            zeros(2, 2),
            Variable(3, 3),
            3,
            3,
        )
        @test_throws DimensionMismatch RelativeEntropyEpiCone(
            Variable(2, 2),
            zeros(3, 3),
            3,
            3,
        )
        @test_throws DimensionMismatch RelativeEntropyEpiCone(
            Variable(2, 2),
            Variable(2, 3),
            3,
            3,
        )

        @test_throws DimensionMismatch RelativeEntropyEpiCone(
            zeros(3, 3),
            zeros(3, 3),
            3,
            3,
            zeros(2, 2),
        )
        @test_throws DimensionMismatch RelativeEntropyEpiCone(
            zeros(3, 3),
            Variable(3, 3),
            3,
            3,
            zeros(2, 2),
        )
        @test_throws DimensionMismatch RelativeEntropyEpiCone(
            Variable(3, 3),
            zeros(3, 3),
            3,
            3,
            zeros(2, 2),
        )
        @test_throws DimensionMismatch RelativeEntropyEpiCone(
            Variable(3, 3),
            Variable(2, 3),
            3,
            3,
            zeros(2, 2),
        )

        @test_throws DimensionMismatch RelativeEntropyEpiCone(
            zeros(3, 3),
            zeros(3, 3),
            3,
            3,
            zeros(2, 3),
        )
        @test_throws DimensionMismatch RelativeEntropyEpiCone(
            zeros(3, 3),
            Variable(3, 3),
            3,
            3,
            zeros(2, 3),
        )
        @test_throws DimensionMismatch RelativeEntropyEpiCone(
            Variable(3, 3),
            zeros(3, 3),
            3,
            3,
            zeros(2, 3),
        )
        @test_throws DimensionMismatch RelativeEntropyEpiCone(
            Variable(3, 3),
            Variable(2, 3),
            3,
            3,
            zeros(2, 3),
        )

        @test_throws DimensionMismatch RelativeEntropyEpiCone(
            zeros(3, 3),
            zeros(3, 3),
            3,
            3,
            zeros(2),
        )
        @test_throws DimensionMismatch RelativeEntropyEpiCone(
            zeros(3, 3),
            Variable(3, 3),
            3,
            3,
            zeros(2),
        )
        @test_throws DimensionMismatch RelativeEntropyEpiCone(
            Variable(3, 3),
            zeros(3, 3),
            3,
            3,
            zeros(2),
        )
        @test_throws DimensionMismatch RelativeEntropyEpiCone(
            Variable(3, 3),
            Variable(2, 3),
            3,
            3,
            zeros(2),
        )

        @test_throws DimensionMismatch Variable(2, 2) in RelativeEntropyEpiCone(
            Variable(3, 3),
            Variable(3, 3),
            3,
            3,
        )
    end
end

@add_problem sdp function sdp_relative_entropy(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    for cplx in [false, true]
        n = 4
        X = eye(n)
        if cplx
            Y = ComplexVariable(n, n)
        else
            Y = Variable(n, n)
        end

        c1 = eye(n) in RelativeEntropyEpiCone(X, Y)
        objective = real(tr(Y))
        p = minimize(objective, c1; numeric_type = T)

        handle_problem!(p)

        if test
            @test Y.value ≈ eye(n) * exp(-1) atol = atol rtol = rtol
        end
    end
end

@add_problem sdp function sdp_geom_mean_hypocone_real_0(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    n = 4
    A = Variable(n, n)
    B = randn(n, n)
    B = B * B' # now A is positive semidefinite
    B += 0.2 * I # prevent numerical instability

    c1 = eye(n) in GeomMeanHypoCone(A, B, 0)
    objective = tr(A)
    p = minimize(objective, c1; numeric_type = T)

    handle_problem!(p)

    if test
        @test A.value ≈ eye(n) atol = atol rtol = rtol
    end
end

@add_problem sdp function sdp_geom_mean_hypocone_real_1(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    n = 4
    A = Variable(n, n)
    B = randn(n, n)
    B = B * B' # now A is positive semidefinite
    B += 0.2 * I # prevent numerical instability

    c1 = eye(n) in GeomMeanHypoCone(B, A, 1)
    objective = tr(A)
    p = minimize(objective, c1; numeric_type = T)

    handle_problem!(p)

    if test
        @test A.value ≈ eye(n) atol = atol rtol = rtol
    end
end

@add_problem sdp function sdp_geom_mean_hypocone_real_1_2(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    n = 4
    A = randn(n, n)
    A = A * A' # now A is positive semidefinite
    A += 0.2 * I # prevent numerical instability
    B = Variable(n, n)

    c1 = eye(n) in GeomMeanHypoCone(A, B, 1 // 2)
    objective = tr(B)
    p = minimize(objective, c1; numeric_type = T)

    handle_problem!(p)

    if test
        @test B.value ≈ A^-1 atol = atol rtol = rtol
    end
end

@add_problem sdp function sdp_geom_mean_hypocone_real_3_8(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    n = 4
    A = randn(n, n)
    A = A * A' # now A is positive semidefinite
    A += 0.2 * I # prevent numerical instability
    B = Variable(n, n)

    c1 = eye(n) in GeomMeanHypoCone(A, B, 3 // 8)
    objective = tr(B)
    p = minimize(objective, c1; numeric_type = T)

    handle_problem!(p)

    # A #_t B = I  =>  B = A^(1-1/t)
    if test
        @test B.value ≈ A^(-5 // 3) atol = atol rtol = rtol
    end
end

@add_problem sdp function sdp_geom_mean_hypocone_real_3_5(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    n = 4
    A = randn(n, n)
    A = A * A' # now A is positive semidefinite
    A += 0.2 * I # prevent numerical instability
    B = Variable(n, n)

    c1 = eye(n) in GeomMeanHypoCone(A, B, 3 // 5)
    objective = tr(B)
    p = minimize(objective, c1; numeric_type = T)

    handle_problem!(p)

    # A #_t B = I  =>  B = A^(1-1/t)
    if test
        @test B.value ≈ A^(-2 // 3) atol = atol rtol = rtol
    end
end

@add_problem sdp function sdp_geom_mean_hypocone_cplx_1_2(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    n = 4
    A = randn(ComplexF64, n, n)
    A = A * A' # now A is positive semidefinite
    A += 0.2 * I # prevent numerical instability
    B = ComplexVariable(n, n)

    c1 = eye(n) in GeomMeanHypoCone(A, B, 1 // 2)
    objective = real(tr(B))
    p = minimize(objective, c1; numeric_type = T)

    handle_problem!(p)

    if test
        @test real.(B.value) ≈ real.(A^-1) atol = atol rtol = rtol
        @test imag.(B.value) ≈ imag.(A^-1) atol = atol rtol = rtol
    end
end

@add_problem sdp function sdp_geom_mean_hypocone_cplx_3_8(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    n = 4
    A = randn(ComplexF64, n, n)
    A = A * A' # now A is positive semidefinite
    A += 0.2 * I # prevent numerical instability
    B = ComplexVariable(n, n)

    c1 = eye(n) in GeomMeanHypoCone(A, B, 3 // 8)
    objective = real(tr(B))
    p = minimize(objective, c1; numeric_type = T)

    handle_problem!(p)

    # A #_t B = I  =>  B = A^(1-1/t)
    if test
        @test real.(B.value) ≈ real.(A^(-5 // 3)) atol = atol rtol = rtol
        @test imag.(B.value) ≈ imag.(A^(-5 // 3)) atol = atol rtol = rtol
    end
end

@add_problem sdp function sdp_geom_mean_hypocone_cplx_3_5(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    n = 4
    A = randn(ComplexF64, n, n)
    A = A * A' # now A is positive semidefinite
    A += 0.2 * I # prevent numerical instability
    B = ComplexVariable(n, n)

    c1 = eye(n) in GeomMeanHypoCone(A, B, 3 // 5)
    objective = real(tr(B))
    p = minimize(objective, c1; numeric_type = T)

    handle_problem!(p)

    # A #_t B = I  =>  B = A^(1-1/t)
    if test
        @test real.(B.value) ≈ real.(A^(-2 // 3)) atol = atol rtol = rtol
        @test imag.(B.value) ≈ imag.(A^(-2 // 3)) atol = atol rtol = rtol
    end
end

@add_problem sdp function sdp_geom_mean_hypocone_fullhyp(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    # These matrices satisfy A^(1/2) ⪰ B and so (A,eye(2),B) \in hyp_{1/2}
    # A and B however do not satisfy A ⪰ B^2 and so [A B; B eye(2)] not
    # psd, and using fullhyp=false will give an infeasible SDP.
    A = [6.25 0; 0 16]
    B = [2 1; 1 2]
    S = Semidefinite(2)

    p = minimize(
        0,
        [B in GeomMeanHypoCone(A, eye(2), 1 // 2, false)];
        numeric_type = T,
    )
    handle_problem!(p)
    if test
        @test p.status == MOI.INFEASIBLE
    end

    p = minimize(
        0,
        [B in GeomMeanHypoCone(A, eye(2), 1 // 2)];
        numeric_type = T,
    )
    handle_problem!(p)
    if test
        @test p.status == MOI.OPTIMAL
    end
end

@add_problem sdp function sdp_geom_mean_epicone_real_neg1_optA(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    n = 3
    B = randn(n, n)
    B = B * B' # now B is positive semidefinite
    B += 0.2 * I # prevent numerical instability
    A = Variable(n, n)

    c1 = eye(n) in GeomMeanEpiCone(A, B, -1)
    objective = tr(A)
    p = maximize(objective, c1; numeric_type = T)

    handle_problem!(p)

    # A #_t B = I  =>  B = A^(1-1/t)
    if test
        @test B ≈ A.value^2 atol = atol rtol = rtol
    end
end

@add_problem sdp function sdp_geom_mean_epicone_real_neg1_optB(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    n = 3
    A = randn(n, n)
    A = A * A' # now A is positive semidefinite
    A += 0.2 * I # prevent numerical instability
    A /= tr(A) # solver has problems if B is large
    B = Variable(n, n)

    c1 = eye(n) in GeomMeanEpiCone(A, B, -1)
    objective = tr(B)
    p = minimize(objective, c1; numeric_type = T)

    handle_problem!(p)

    # A #_t B = I  =>  B = A^(1-1/t)
    if test
        @test B.value ≈ A^2 atol = atol rtol = rtol
    end
end

@add_problem sdp function sdp_geom_mean_epicone_real_neg3_5_optA(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    n = 3
    B = randn(n, n)
    B = B * B' # now B is positive semidefinite
    B += 0.2 * I # prevent numerical instability
    A = Variable(n, n)

    c1 = eye(n) in GeomMeanEpiCone(A, B, -3 // 5)
    objective = tr(A)
    p = maximize(objective, c1; numeric_type = T)

    handle_problem!(p)

    # A #_t B = I  =>  B = A^(1-1/t)
    if test
        @test B ≈ A.value^(8 // 3) atol = atol rtol = rtol
    end
end

@add_problem sdp function sdp_geom_mean_epicone_real_neg3_5_optB(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    n = 3
    A = randn(n, n)
    A = A * A' # now A is positive semidefinite
    A += 0.2 * I # prevent numerical instability
    A /= tr(A) # solver has problems if B is large
    B = Variable(n, n)

    c1 = eye(n) in GeomMeanEpiCone(A, B, -3 // 5)
    objective = tr(B)
    p = minimize(objective, c1; numeric_type = T)

    handle_problem!(p)

    # A #_t B = I  =>  B = A^(1-1/t)
    if test
        @test B.value ≈ A^(8 // 3) atol = atol rtol = rtol
    end
end

@add_problem sdp function sdp_geom_mean_epicone_real_8_5_optA(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    n = 3
    B = randn(n, n)
    B = B * B' # now B is positive semidefinite
    B += 0.2 * I # prevent numerical instability
    B /= tr(B) # solver has problems if B is large
    A = Variable(n, n)

    c1 = eye(n) in GeomMeanEpiCone(A, B, 8 // 5)
    objective = tr(A)
    p = minimize(objective, c1; numeric_type = T)

    handle_problem!(p)

    # A #_t B = I  =>  B = A^(1-1/t)
    if test
        @test B ≈ A.value^(3 // 8) atol = atol rtol = rtol
    end
end

@add_problem sdp function sdp_geom_mean_epicone_real_8_5_optB(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    n = 3
    A = randn(n, n)
    A = A * A' # now A is positive semidefinite
    A += 0.2 * I # prevent numerical instability
    B = Variable(n, n)

    c1 = eye(n) in GeomMeanEpiCone(A, B, 8 // 5)
    objective = tr(B)
    p = maximize(objective, c1; numeric_type = T)

    handle_problem!(p)

    # A #_t B = I  =>  B = A^(1-1/t)
    if test
        @test B.value ≈ A^(3 // 8) atol = atol rtol = rtol
    end
end

@add_problem sdp function sdp_quantum_relative_entropy_argcheck(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    if test
        @test_throws DimensionMismatch quantum_relative_entropy(
            zeros(2, 3),
            zeros(2, 3),
        )
        @test_throws DimensionMismatch quantum_relative_entropy(
            zeros(2, 3),
            Variable(2, 3),
        )
        @test_throws DimensionMismatch quantum_relative_entropy(
            Variable(2, 3),
            zeros(2, 3),
        )
        @test_throws DimensionMismatch quantum_relative_entropy(
            Variable(2, 3),
            Variable(2, 3),
        )

        @test_throws DimensionMismatch quantum_relative_entropy(
            zeros(2, 2),
            zeros(3, 3),
        )
        @test_throws DimensionMismatch quantum_relative_entropy(
            zeros(2, 2),
            Variable(3, 3),
        )
        @test_throws DimensionMismatch quantum_relative_entropy(
            Variable(2, 2),
            zeros(3, 3),
        )
        @test_throws DimensionMismatch quantum_relative_entropy(
            Variable(2, 2),
            Variable(2, 3),
        )

        z = zeros(2, 2)
        v = Variable(2, 2)
        nh = [1 1; 0 1] # not hermitian
        np = [1 2; 2 1] # not positive semidefinite
        @test_throws DomainError quantum_relative_entropy(z, nh)
        @test_throws DomainError quantum_relative_entropy(z, np)
        @test_throws DomainError quantum_relative_entropy(nh, z)
        @test_throws DomainError quantum_relative_entropy(np, z)
        @test_throws DomainError quantum_relative_entropy(v, nh)
        @test_throws DomainError quantum_relative_entropy(v, np)
        @test_throws DomainError quantum_relative_entropy(nh, v)
        @test_throws DomainError quantum_relative_entropy(np, v)
    end
end

function sdp_quantum_relative_entropy_impl(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
    lowrank::Bool,
    mode::Integer,
) where {T,test}
    # Too slow for n>3.  Anything that can be done to speed it up?
    n = 3
    X = randn(ComplexF64, n, n)
    if lowrank
        X[:, 1] .= 0
    end
    X = X * X' # now X is positive semidefinite
    X /= tr(X) # make it a quantum state

    A = ComplexVariable(n, n)
    B = ComplexVariable(n, n)
    constraints = Array{Constraint,1}([A == X])
    if mode != 5
        push!(constraints, tr(B) == 1)
    end

    if mode == 1
        objective = quantum_relative_entropy(A, B)
    elseif mode == 2
        objective = quantum_relative_entropy(X, B)
    elseif mode == 3
        objective = quantum_relative_entropy(B, A)
    elseif mode == 4
        objective = quantum_relative_entropy(B, X)
    elseif mode == 5
        objective = quantum_relative_entropy(X, X)
    end

    p = minimize(objective, constraints; numeric_type = T)
    handle_problem!(p)

    if test
        @test real.(evaluate(A)) ≈ real.(X) atol = atol rtol = rtol
        @test imag.(evaluate(A)) ≈ imag.(X) atol = atol rtol = rtol
        if mode != 5
            @test real.(evaluate(B)) ≈ real.(X) atol = atol rtol = rtol
            @test imag.(evaluate(B)) ≈ imag.(X) atol = atol rtol = rtol
        end
        @test p.optval ≈ 0 atol = atol rtol = rtol
        if mode == 1
            @test p.optval ≈ evaluate(quantum_relative_entropy(A, B)) atol =
                atol rtol = rtol
        elseif mode == 2
            @test p.optval ≈ evaluate(quantum_relative_entropy(X, B)) atol =
                atol rtol = rtol
        elseif mode == 3
            @test p.optval ≈ evaluate(quantum_relative_entropy(B, A)) atol =
                atol rtol = rtol
        elseif mode == 4
            @test p.optval ≈ evaluate(quantum_relative_entropy(B, X)) atol =
                atol rtol = rtol
        elseif mode == 5
            @test p.optval ≈ evaluate(quantum_relative_entropy(X, X)) atol =
                atol rtol = rtol
        end
    end
end

@add_problem sdp function sdp_quantum_relative_entropy1_fullrank(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    return sdp_quantum_relative_entropy_impl(
        handle_problem!,
        Val(test),
        atol,
        rtol,
        T,
        false,
        1,
    )
end

@add_problem sdp function sdp_quantum_relative_entropy2_fullrank(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    return sdp_quantum_relative_entropy_impl(
        handle_problem!,
        Val(test),
        atol,
        rtol,
        T,
        false,
        2,
    )
end

@add_problem sdp function sdp_quantum_relative_entropy3_fullrank(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    return sdp_quantum_relative_entropy_impl(
        handle_problem!,
        Val(test),
        atol,
        rtol,
        T,
        false,
        3,
    )
end

@add_problem sdp function sdp_quantum_relative_entropy4_fullrank(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    return sdp_quantum_relative_entropy_impl(
        handle_problem!,
        Val(test),
        atol,
        rtol,
        T,
        false,
        4,
    )
end

@add_problem sdp function sdp_quantum_relative_entropy5_fullrank(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    return sdp_quantum_relative_entropy_impl(
        handle_problem!,
        Val(test),
        atol,
        rtol,
        T,
        false,
        5,
    )
end

@add_problem sdp function sdp_quantum_relative_entropy1_lowrank(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    return sdp_quantum_relative_entropy_impl(
        handle_problem!,
        Val(test),
        atol,
        rtol,
        T,
        true,
        1,
    )
end

@add_problem sdp function sdp_quantum_relative_entropy2_lowrank(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    return sdp_quantum_relative_entropy_impl(
        handle_problem!,
        Val(test),
        atol,
        rtol,
        T,
        true,
        2,
    )
end

@add_problem sdp function sdp_quantum_relative_entropy3_lowrank(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    return sdp_quantum_relative_entropy_impl(
        handle_problem!,
        Val(test),
        atol,
        rtol,
        T,
        true,
        3,
    )
end

@add_problem sdp function sdp_quantum_relative_entropy4_lowrank(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    return sdp_quantum_relative_entropy_impl(
        handle_problem!,
        Val(test),
        atol,
        rtol,
        T,
        true,
        4,
    )
end

@add_problem sdp function sdp_quantum_relative_entropy5_lowrank(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    return sdp_quantum_relative_entropy_impl(
        handle_problem!,
        Val(test),
        atol,
        rtol,
        T,
        true,
        5,
    )
end

@add_problem sdp function sdp_quantum_relative_entropy_const(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    if test
        @test quantum_relative_entropy(
            diagm(0 => [0.5, 0.5]),
            diagm(0 => [0.5, 0.5]),
        ) ≈ 0 atol = atol rtol = rtol
        @test quantum_relative_entropy(
            diagm(0 => [1.0, 0.0]),
            diagm(0 => [1.0, 0.0]),
        ) ≈ 0 atol = atol rtol = rtol
        @test quantum_relative_entropy(
            diagm(0 => [1.0, 0.0]),
            diagm(0 => [0.5, 0.5]),
        ) ≈ log(2) atol = atol rtol = rtol
        @test isinf(
            quantum_relative_entropy(
                diagm(0 => [0.5, 0.5]),
                diagm(0 => [1.0, 0.0]),
            ),
        )
    end
end

@add_problem sdp function sdp_quantum_channel_capacity(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    # Capacity of classical-quantum channel.  See chapter 3 of arxiv:1705.06671,
    # Fawzi, Fawzi "Efficient optimization of the quantum relative entropy".
    # For channels with two symbols, both pure states, we have an exact formula to compare against.

    n = 3

    for cplx in [false, true]
        if cplx
            ψ1 = normalize(randn(ComplexF64, n))
            ψ2 = normalize(randn(ComplexF64, n))
        else
            ψ1 = normalize(randn(Float64, n))
            ψ2 = normalize(randn(Float64, n))
        end
        ρ1 = ψ1 * ψ1'
        ρ2 = ψ2 * ψ2'

        p1 = Variable()
        p2 = Variable()

        objective =
            (
                quantum_entropy(p1 * ρ1 + p2 * ρ2) - p1 * quantum_entropy(ρ1) -
                p2 * quantum_entropy(ρ2)
            ) / log(2)
        constraints = [p1 >= 0, p2 >= 0, p1 + p2 == 1]
        p = maximize(objective, constraints; numeric_type = T)
        handle_problem!(p)

        ϵ = abs(ψ1' * ψ2)
        q = (1 + ϵ) / 2
        r = [q, 1 - q]
        v = -r' * log2.(r)

        if test
            @test p.optval ≈ v atol = atol rtol = rtol
            @test evaluate(objective) ≈ v atol = atol rtol = rtol

            @test_throws DimensionMismatch quantum_entropy(Variable(2, 3))
        end
    end
end

@add_problem sdp function sdp_trace_mpower_argcheck(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    if test
        @test_throws DimensionMismatch trace_mpower(
            zeros(2, 3),
            1 // 2,
            zeros(2, 3),
        )
        @test_throws DimensionMismatch trace_mpower(
            Variable(2, 3),
            1 // 2,
            zeros(2, 3),
        )

        @test_throws DimensionMismatch trace_mpower(
            zeros(2, 2),
            1 // 2,
            zeros(3, 3),
        )
        @test_throws DimensionMismatch trace_mpower(
            Variable(2, 2),
            1 // 2,
            zeros(3, 3),
        )

        @test_throws DomainError trace_mpower(
            Variable(3, 3),
            5 // 2,
            zeros(3, 3),
        )

        nh = [1 1; 0 1] # not hermitian
        np = [1 2; 2 1] # not positive semidefinite
        @test_throws DomainError trace_mpower(Variable(2, 2), 1 // 2, nh)
        @test_throws DomainError trace_mpower(Variable(2, 2), 1 // 2, np)
    end
end

@add_problem sdp function sdp_trace_mpower_real_2_3(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    # Too slow for n>3.  Anything that can be done to speed it up?
    n = 3
    A = randn(n, n)
    A = A * A' # now A is positive semidefinite
    C = randn(n, n)
    C = C * C' # now C is positive semidefinite
    t = 2 // 3

    B = Semidefinite(n, n)

    constraints = [B ⪯ A]
    objective = trace_mpower(B, t, C)
    p = maximize(objective, constraints; numeric_type = T)

    handle_problem!(p)

    if test
        @test real.(B.value) ≈ real.(A) atol = atol rtol = rtol
        @test imag.(B.value) ≈ imag.(A) atol = atol rtol = rtol
        @test p.optval ≈ tr(C * A^t) atol = atol rtol = rtol
        @test p.optval ≈ trace_mpower(A, t, C) atol = atol rtol = rtol
        @test p.optval ≈ evaluate(trace_mpower(B, t, C)) atol = atol rtol = rtol
    end
end

@add_problem sdp function sdp_trace_mpower_cplx_2_3(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    # Too slow for n>3.  Anything that can be done to speed it up?
    n = 3
    A = randn(ComplexF64, n, n)
    A = A * A' # now A is positive semidefinite
    C = randn(ComplexF64, n, n)
    C = C * C' # now C is positive semidefinite
    t = 2 // 3

    B = HermitianSemidefinite(n, n)

    constraints = [B ⪯ A]
    objective = trace_mpower(B, t, C)
    p = maximize(objective, constraints; numeric_type = T)

    handle_problem!(p)

    if test
        @test real.(B.value) ≈ real.(A) atol = atol rtol = rtol
        @test imag.(B.value) ≈ imag.(A) atol = atol rtol = rtol
        @test p.optval ≈ tr(C * A^t) atol = atol rtol = rtol
        @test p.optval ≈ trace_mpower(A, t, C) atol = atol rtol = rtol
        @test p.optval ≈ evaluate(trace_mpower(B, t, C)) atol = atol rtol = rtol
    end
end

@add_problem sdp function sdp_trace_mpower_real_5_4(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    # Too slow for n>3.  Anything that can be done to speed it up?
    n = 3
    A = randn(n, n)
    A = A * A' # now A is positive semidefinite
    C = randn(n, n)
    C = C * C' # now C is positive semidefinite
    t = 5 // 4

    B = Semidefinite(n, n)

    constraints = [
        # note: B ⪰ A doesn't work for this unit test because trace_mpower is not operator monotone
        # unless t ∈ [0,1].
        B == A,
    ]
    objective = trace_mpower(B, t, C)
    p = minimize(objective, constraints; numeric_type = T)

    handle_problem!(p)

    if test
        @test real.(B.value) ≈ real.(A) atol = atol rtol = rtol
        @test imag.(B.value) ≈ imag.(A) atol = atol rtol = rtol
        @test p.optval ≈ tr(C * A^t) atol = atol rtol = rtol
        @test p.optval ≈ trace_mpower(A, t, C) atol = atol rtol = rtol
        @test p.optval ≈ evaluate(trace_mpower(B, t, C)) atol = atol rtol = rtol
    end
end

@add_problem sdp function sdp_trace_mpower_cplx_5_4(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    # Too slow for n>3.  Anything that can be done to speed it up?
    n = 3
    A = randn(ComplexF64, n, n)
    A = A * A' # now A is positive semidefinite
    C = randn(ComplexF64, n, n)
    C = C * C' # now C is positive semidefinite
    t = 5 // 4

    B = HermitianSemidefinite(n, n)

    constraints = [
        # note: B ⪰ A doesn't work for this unit test because trace_mpower is not operator monotone
        # unless t ∈ [0,1].
        B == A,
    ]
    objective = trace_mpower(B, t, C)
    p = minimize(objective, constraints; numeric_type = T)

    handle_problem!(p)

    if test
        @test real.(B.value) ≈ real.(A) atol = atol rtol = rtol
        @test imag.(B.value) ≈ imag.(A) atol = atol rtol = rtol
        @test p.optval ≈ tr(C * A^t) atol = atol rtol = rtol
        @test p.optval ≈ trace_mpower(A, t, C) atol = atol rtol = rtol
        @test p.optval ≈ evaluate(trace_mpower(B, t, C)) atol = atol rtol = rtol
    end
end

@add_problem sdp function sdp_trace_mpower_real_neg1_4(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    # Too slow for n>3.  Anything that can be done to speed it up?
    n = 3
    A = randn(n, n)
    A = A * A' # now A is positive semidefinite
    A += 0.2 * I # prevent numerical instability
    C = randn(n, n)
    C = C * C' # now C is positive semidefinite
    t = -1 // 4

    B = Semidefinite(n, n)

    constraints = [
        # -trace_mpower is operator monotone for t ∈ [-1, 0].
        B ⪯ A,
    ]
    objective = trace_mpower(B, t, C)
    p = minimize(objective, constraints; numeric_type = T)

    handle_problem!(p)

    if test
        @test real.(B.value) ≈ real.(A) atol = atol rtol = rtol
        @test imag.(B.value) ≈ imag.(A) atol = atol rtol = rtol
        @test p.optval ≈ tr(C * A^t) atol = atol rtol = rtol
        @test p.optval ≈ trace_mpower(A, t, C) atol = atol rtol = rtol
        @test p.optval ≈ evaluate(trace_mpower(B, t, C)) atol = atol rtol = rtol
    end
end

@add_problem sdp function sdp_trace_mpower_cplx_neg1_4(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    # Too slow for n>3.  Anything that can be done to speed it up?
    n = 3
    A = randn(ComplexF64, n, n)
    A = A * A' # now A is positive semidefinite
    A += 0.2 * I # prevent numerical instability
    C = randn(ComplexF64, n, n)
    C = C * C' # now C is positive semidefinite
    t = -1 // 4

    B = HermitianSemidefinite(n, n)

    constraints = [
        # -trace_mpower is operator monotone for t ∈ [-1, 0].
        B ⪯ A,
    ]
    objective = trace_mpower(B, t, C)
    p = minimize(objective, constraints; numeric_type = T)

    handle_problem!(p)

    if test
        @test real.(B.value) ≈ real.(A) atol = atol rtol = rtol
        @test imag.(B.value) ≈ imag.(A) atol = atol rtol = rtol
        @test p.optval ≈ tr(C * A^t) atol = atol rtol = rtol
        @test p.optval ≈ trace_mpower(A, t, C) atol = atol rtol = rtol
        @test p.optval ≈ evaluate(trace_mpower(B, t, C)) atol = atol rtol = rtol
    end
end

@add_problem sdp function sdp_trace_logm_argcheck(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    if test
        @test_throws DimensionMismatch trace_logm(zeros(2, 3), zeros(2, 3))
        @test_throws DimensionMismatch trace_logm(Variable(2, 3), zeros(2, 3))

        @test_throws DimensionMismatch trace_logm(zeros(2, 2), zeros(3, 3))
        @test_throws DimensionMismatch trace_logm(Variable(2, 2), zeros(3, 3))

        nh = [1 1; 0 1] # not hermitian
        np = [1 2; 2 1] # not positive semidefinite
        @test_throws DomainError trace_logm(Variable(2, 2), nh)
        @test_throws DomainError trace_logm(Variable(2, 2), np)
    end
end

@add_problem sdp function sdp_trace_logm_real(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    n = 3
    C = randn(Float64, n, n)
    C = C * C' # now C is positive semidefinite

    X = Variable(n, n)

    constraints = [X ⪯ eye(n)]
    objective = trace_logm(X, C)
    p = maximize(objective, constraints; numeric_type = T)

    handle_problem!(p)

    if test
        @test real.(X.value) ≈ real.(eye(n)) atol = atol rtol = rtol
        @test imag.(X.value) ≈ imag.(eye(n)) atol = atol rtol = rtol
        @test p.optval ≈ tr(C * log(evaluate(X))) atol = atol rtol = rtol
        @test p.optval ≈ trace_logm(evaluate(X), C) atol = atol rtol = rtol
        @test p.optval ≈ evaluate(trace_logm(X, C)) atol = atol rtol = rtol
    end
end

@add_problem sdp function sdp_trace_logm_cplx(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    n = 3
    C = randn(ComplexF64, n, n)
    C = C * C' # now C is positive semidefinite

    X = ComplexVariable(n, n)

    constraints = [X ⪯ eye(n)]
    objective = trace_logm(X, C)
    p = maximize(objective, constraints; numeric_type = T)

    handle_problem!(p)

    if test
        @test real.(X.value) ≈ real.(eye(n)) atol = atol rtol = rtol
        @test imag.(X.value) ≈ imag.(eye(n)) atol = atol rtol = rtol
        @test p.optval ≈ tr(C * log(evaluate(X))) atol = atol rtol = rtol
        @test p.optval ≈ trace_logm(evaluate(X), C) atol = atol rtol = rtol
        @test p.optval ≈ evaluate(trace_logm(X, C)) atol = atol rtol = rtol
    end
end

@add_problem sdp function sdp_lieb_ando(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    if test
        @test_throws DomainError lieb_ando(
            zeros(3, 3),
            zeros(3, 3),
            eye(3),
            -3 // 2,
        )
        @test_throws DomainError lieb_ando(
            zeros(3, 3),
            zeros(3, 3),
            eye(3),
            5 // 2,
        )
        @test_throws DomainError lieb_ando(
            zeros(3, 3),
            Variable(3, 3),
            eye(3),
            -3 // 2,
        )
        @test_throws DomainError lieb_ando(
            zeros(3, 3),
            Variable(3, 3),
            eye(3),
            5 // 2,
        )
        @test_throws DomainError lieb_ando(
            Variable(3, 3),
            zeros(3, 3),
            eye(3),
            -3 // 2,
        )
        @test_throws DomainError lieb_ando(
            Variable(3, 3),
            zeros(3, 3),
            eye(3),
            5 // 2,
        )
        @test_throws DomainError lieb_ando(
            Variable(3, 3),
            Variable(3, 3),
            eye(3),
            -3 // 2,
        )
        @test_throws DomainError lieb_ando(
            Variable(3, 3),
            Variable(3, 3),
            eye(3),
            5 // 2,
        )
    end

    for n in [2, 3]
        for t in [1 // 2, 1 // 4, 3 // 4, 1 // 8, 3 // 2, 5 // 4]
            for cplx in [false, true]
                #@show n,t,cplx

                if cplx
                    A = randn(ComplexF64, n, n)
                    B = randn(ComplexF64, n, n)
                    X = ComplexVariable(n, n)
                    Y = ComplexVariable(n, n)
                else
                    A = randn(n, n)
                    B = randn(n, n)
                    X = Variable(n, n)
                    Y = Variable(n, n)
                end
                A = A * A'
                B = B * B'
                A = A / tr(A)
                B = B / tr(B)

                QtAB = lieb_ando(A, B, eye(n), t)

                objective = lieb_ando(X, B, eye(n), t)
                if t >= 0 && t <= 1
                    p = maximize(objective, [X == A]; numeric_type = T)
                else
                    p = minimize(objective, [X == A]; numeric_type = T)
                end
                handle_problem!(p)
                if test
                    @test p.optval ≈ QtAB atol = atol * 5 rtol = rtol
                end

                objective = lieb_ando(A, Y, eye(n), t)
                if t >= 0 && t <= 1
                    p = maximize(objective, [Y == B]; numeric_type = T)
                else
                    p = minimize(objective, [Y == B]; numeric_type = T)
                end
                handle_problem!(p)
                if test
                    @test p.optval ≈ QtAB atol = atol * 5 rtol = rtol
                end

                objective = lieb_ando(X, Y, eye(n), t)
                if t >= 0 && t <= 1
                    p = maximize(objective, [X == A, Y == B]; numeric_type = T)
                else
                    p = minimize(objective, [X == A, Y == B]; numeric_type = T)
                end
                handle_problem!(p)
                if test
                    @test p.optval ≈ QtAB atol = atol * 5 rtol = rtol
                end
            end
        end
    end
end

@add_problem sdp function sdp_min_maxeig_canon_lmi(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    # Minimize the maximum eigenvalue of an affine matrix function of a vector argument. Formulated
    # as a linear optimization subject to a canonical LMI (SDP) constraint:
    #
    # minimize      λ
    # subject to    A(x) ⪯ λ*I
    # where x is a vector in Rⁿ and the matrix function A(x) = A₀ + A₁x₁ + A₂x₂ + … + Aₙxₙ.
    #
    # Besides serving as a benchmark, it also tests the solution of the Issue
    # https://github.com/jump-dev/Convex.jl/issues/447
    # (the issues boils down to  evaluation of `A₁x₁`, because `x₁` was a matrix of size (1,1)).

    A₀ = [1.0 2.0; 2.0 3.0]
    A₁ = [1.0 0.0; 0.0 -1.0]
    A₂ = [4.0 5.0; 5.0 -6.0]

    x = Variable(2)
    λ = Variable()

    A = A₀ + A₁ * x[1] + A₂ * x[2]

    p = minimize(λ, A ⪯ λ * eye(2); numeric_type = T)

    handle_problem!(p)

    if test
        @test λ.value ≈ 2.400000051025101 atol = atol rtol = rtol
        @test x.value ≈ [3.0000000535867315, -0.4000000018585541] atol = atol rtol =
            rtol
        @test evaluate(A) ≈ [
            2.400000046152515 -9.292770553059881e-9
            -9.292770553059881e-9 2.399999957564593
        ] atol = atol rtol = rtol
    end
end


text:
"""
Stream ciphers produce a pseudorandom stream of bits called the
keystream. The keystream is XORed to a plaintext to encrypt it and then
XORed again to the ciphertext to decrypt it.

From a high-level perspective, there are two types of stream ciphers:
stateful and counter based. Stateful stream ciphers have a secret internal
state that evolves throughout keystream generation.
Counter-based stream ciphers produce chunks of keystream from a key, a
nonce, and a counter value.

A feedback shift register (FSR) is simply an array of bits equipped with an update
feedback function, which I’ll denote as f. FSRs have a period.
Linear feedback shift registers (LFSRs) are FSRs with a linear feedback
function - namely, a function that’s the XOR of some bits of the state.
To mitigate the insecurity of LFSRs, you can hide their linearity by
passing their output bits through a nonlinear function before returning
them to produce what is called a filtered LFSR.
Nonlinear FSRs (NFSRs) are like LFSRs but with a nonlinear feedback
function instead of a linear one. That is, instead of just bitwise XORs, the
feedback function can include bitwise AND and OR operations.

    Grain-128a combines NFSR and LFSR

Salsa20
RC4

To be primitive, the
polynomial must have the following qualities:
The polynomial must be irreducible, meaning that it can’t be
factorized; that is, written as a product of smaller polynomials. For
example, X + X 3 is not irreducible because it’s equal to (1 + X)(X +
X2):
(1 + X)(X + X2) = X + X2 + X2 + X3 = X + X3
The polynomial must satisfy certain other mathematical properties
that cannot be easily explained without nontrivial mathematical
notions but are easy to test.
"""
struct SteamCipher
    key::Vector{UInt8}
    nonce::Vector{UInt8}
 end

 struct ChaCha20
    key::Vector{UInt32}  # 256 bit
    nonce::Vector{UInt32}  # 96 bit
 end

 ChaCha20(key::Vector{UInt8}, nonce::Vector{UInt8}) = ChaCha20(reinterpret(UInt32, key), reinterpret(UInt32, nonce))

 function encrypt(cipher::ChaCha20, plaintext::Vector{UInt8}, rounds::Int64=10)
    function quarter_round!(s, a, b, c, d)
        s[a] += s[b]; s[d] ⊻= s[a]; s[d] = (s[d] << 16) + (s[d] >> (32-16));
        s[c] += s[d]; s[b] ⊻= s[c]; s[b] = (s[b] << 12) + (s[b] >> (32-12));
        s[a] += s[b]; s[d] ⊻= s[a]; s[d] = (s[d] << 8) + (s[d] >> (32-8));
        s[c] += s[d]; s[b] ⊻= s[c]; s[b] = (s[b] << 7) + (s[b] >> (32-7));
    end

    function double_round!(s)
        quarter_round!(s, 1, 5,  9, 13)	# 1st column
        quarter_round!(s, 2, 6, 10, 14)	# 2nd column
        quarter_round!(s, 3, 7, 11, 15)	# 3rd column
        quarter_round!(s, 4, 8, 12, 16)	# 4th column

        quarter_round!(s, 1, 6, 11, 16)	# diagonal 1 (main diagonal)
        quarter_round!(s, 2, 7, 12, 13)	# diagonal 2
        quarter_round!(s, 3, 8,  9, 14)	# diagonal 3
        quarter_round!(s, 4, 5, 10, 15) # diagonal 4
    end

    k = cipher.key
    n = cipher.nonce

    result = UInt8[]
    for counter in 0x00000001:cld(length(plaintext), 64)
        s0 = [
            0x61707865, 0x3320646e, 0x79622d32, 0x6b206574,
            k[1], k[2], k[3], k[4],
            k[5], k[6], k[7], k[8],
            counter, n[1], n[2], n[3],
        ]
        s = copy(s0)

        for _ in 1:2:rounds
            double_round!(s)
        end

        for i in 1:length(s)
            append!(result, reinterpret(UInt8, [s[i] + s0[i]]))
        end
    end

    ciphertext = xor.(plaintext, result[1:length(plaintext)])
    return ciphertext
 end

 function decrypt(cipher::ChaCha20, ciphertext::Vector{UInt8}, rounds::Int64=10)
    return encrypt(cipher, ciphertext, 20-rounds)
 end

ext: jl


text:
# Errors and exceptions

export DropboxError
struct DropboxError
    dict::Dict
end

# Metadata

export MediaInfo, SymlinkInfo, FileSharingInfo, PropertyGroup
struct MediaInfo end            # TODO
struct SymlinkInfo end          # TODO
struct FileSharingInfo end      # TODO
struct PropertyGroup end        # TODO

export FolderSharingInfo
struct FolderSharingInfo end    # TODO

export Metadata
abstract type Metadata end
Metadata(d::Dict) = Dict(
    "file" => FileMetadata,
    "folder" => FolderMetadata,
    "deleted" => DeletedMetadata,
)[d[".tag"]](d)

export FileMetadata
struct FileMetadata <: Metadata
    name::String
    id::String
    client_modified::String
    server_modified::String
    rev::String
    size::Int64
    path_lower::Union{Nothing, String}
    path_display::Union{Nothing, String}
    media_info::Union{Nothing, MediaInfo}
    symlink_info::Union{Nothing, SymlinkInfo}
    sharing_info::Union{Nothing, FileSharingInfo}
    property_groups::Union{Nothing, Vector{PropertyGroup}}
    has_explicit_shared_members::Union{Nothing, Bool}
    content_hash::Union{Nothing, String}
end
FileMetadata(d::Dict) = FileMetadata(
    d["name"],
    d["id"],
    d["client_modified"],
    d["server_modified"],
    d["rev"],
    d["size"],
    get(d, "path_lower", nothing),
    get(d, "path_display", nothing),
    nothing,                    # TODO
    nothing,                    # TODO
    nothing,                    # TODO
    nothing,                    # TODO
    get(d, "has_explicit_shared_members", nothing),
    get(d, "content_hash", nothing)
)

export FolderMetadata
struct FolderMetadata <: Metadata
    name::String
    id::String
    path_lower::Union{Nothing, String}
    path_display::Union{Nothing, String}
    sharing_info::Union{Nothing, FolderSharingInfo}
    property_groups::Union{Nothing, Vector{PropertyGroup}}
end
FolderMetadata(d::Dict) = FolderMetadata(
    d["name"],
    d["id"],
    get(d, "path_lower", nothing),
    get(d, "path_display", nothing),
    nothing,                    # TODO
    nothing                     # TODO
)

export DeletedMetadata
struct DeletedMetadata <: Metadata
    name::String
    path_lower::Union{Nothing, String}
    path_display::Union{Nothing, String}
end
DeletedMetadata(d::Dict) = DeletedMetadata(
    d["name"],
    get(d, "path_lower", nothing),
    get(d, "path_display", nothing)
)


text:
import Flugbuch

using Dates

exampleDir = @__DIR__

# in case last tear down did not work
if isfile(joinpath(exampleDir, "runtest.db"))
    rm(joinpath(exampleDir, "runtest.db"))
end

Flugbuch.loadConfig(joinpath(exampleDir, "examples", "flugbuchrc"))
Flugbuch.create(joinpath(exampleDir, "runtest.db"))
Flugbuch.csv(joinpath(exampleDir, "examples", "import.csv"), autoconfirm = true)
Flugbuch.prt()
Flugbuch.prt(beginNr=2, endNr=2)
Flugbuch.prt(beginDate=Date(2018,1,1), endDate=Date(2019,1,1))
Flugbuch.unload()

rm(joinpath(exampleDir, "runtest.db"))

max_stars_repo_name: julianschick/flugbuch

max_issues_repo_name: julianschick/flugbuch

max_forks_repo_name: julianschick/flugbuch


text:
#!/usr/bin/env julia

using Documenter, Logging, NarrativeTest, Test
using HypertextLiteral

subs = NarrativeTest.common_subs()

# Ignore the difference in the output of `print(Int)` between 32-bit and 64-bit platforms.
push!(subs, r"Int64" => s"Int(32|64)")

# Normalize printing of vector types.
if VERSION < v"1.6.0-DEV"
    Base.show_datatype(io::IO, x::Type{Vector{T}}) where {T} = print(io, "Vector{$T}")
end

# Normalize printing of type parameters.
if VERSION < v"1.6.0-DEV"
    function Base.show_datatype(io::IO, x::DataType)
        istuple = x.name === Tuple.name
        if (!isempty(x.parameters) || istuple) && x !== Tuple
            n = length(x.parameters)::Int
            if istuple && n > 3 && all(i -> (x.parameters[1] === i), x.parameters)
                print(io, "NTuple{", n, ", ", x.parameters[1], "}")
            else
                Base.show_type_name(io, x.name)
                print(io, '{')
                for (i, p) in enumerate(x.parameters)
                    show(io, p)
                    i < n && print(io, ", ")
                end
                print(io, '}')
            end
        else
            Base.show_type_name(io, x.name)
        end
    end
end

# Ignore line ending differences for Windows targets.
push!(subs, r"\r\n" => "\n")

# Set the width to 72 so that MD->PDF via pandoc fits the page.
ENV["COLUMNS"] = "72"

package_path(x) = relpath(joinpath(dirname(abspath(PROGRAM_FILE)), "..", x))
default = package_path.(["README.md", "docs/src"])

@testset "HypertextLiteral" begin

    if isempty(ARGS) || "doctest" in ARGS
        @info "Running doctest..."
        DocMeta.setdocmeta!(
            HypertextLiteral,
            :DocTestSetup,
            quote
                using HypertextLiteral:
                    HypertextLiteral, @htl, @htl_str,
                    Reprint, Render, Bypass, EscapeProxy,
                    attribute_value, content, attribute_pair,
                    inside_tag
                using Dates
            end)
        with_logger(Logging.ConsoleLogger(stderr, Logging.Warn)) do
            doctest(HypertextLiteral)
        end
    end

    if isempty(ARGS)
        @info "Running narrative tests..."
        NarrativeTest.testset(; default=default, subs=subs)
    else
        filter!(!=("doctest"), ARGS)
        if !isempty(ARGS)
            @info "Running narrative tests..."
            NarrativeTest.testset(ARGS; subs=subs)
        end
    end
end

@static if VERSION >= v"1.3"

    struct Foo end

    @testset "invalidation" begin
        @test HypertextLiteral.content(Foo()) isa HypertextLiteral.Reprint # fallback

        # Now define a html printing type
        @eval Base.show(io::IO, ::MIME"text/html", ::Foo) = "Foo"

        # Previously this would not have worked because content is a generated function depending on hasmethod in the generator
        @test repr(Base.invokelatest(HypertextLiteral.content, Foo())) == "HypertextLiteral.Render{Foo}(Foo())"
    end
end


text:
export treeverse, treeverse!, TreeverseLog

struct TreeverseAction
    action::Symbol
    τ::Int
    δ::Int
    step::Int
    depth::Int
end

struct TreeverseLog
    actions::Vector{TreeverseAction}
    depth::Base.RefValue{Int}
    peak_mem::Base.RefValue{Int}  # should be `n*(k-1)+2`
end
TreeverseLog() = TreeverseLog(TreeverseAction[], Ref(0), Ref(0))
Base.push!(tlog::TreeverseLog, args...) = push!(tlog.actions, TreeverseAction(args..., tlog.depth[]))

Base.show(io::IO, ::MIME"text/plain", logger::TreeverseLog) = Base.show(io, logger)
function Base.show(io::IO, logger::TreeverseLog)
    print(io, """Treeverse log
| peak memory usage = $(logger.peak_mem[])
| number of function calls = $(count(x->x.action==:call, logger.actions))
| number of gradient calls = $(count(x->x.action==:grad, logger.actions))
| number of stack push/pop = $(count(x->x.action==:store, logger.actions))/$(count(x->x.action==:fetch, logger.actions))""")
end

function binomial_fit(N::Int, δ::Int)
    τ = 1
    while N > binomial(τ+δ, τ)
        τ += 1
    end
    return τ
end

function mid(δ, τ, σ, ϕ)
    κ = ceil(Int, (δ*σ + τ*ϕ)/(τ+δ))
    if κ >= ϕ && δ > 0
        κ = max(σ+1, ϕ-1)
    end
    return κ
end

"""
    treeverse(f, gf, s; δ, N, τ=binomial_fit(N,δ), f_inplace=true, logger = TreeverseLog())

Treeverse algorithm for back-propagating a program memory efficiently.

Positional arguments
* `f`, the step function that ``s_{i+1} = f(s_i)``,
* `gf`, the single step gradient function that ``g_i = gf(s_i, g_{i+1})``
   !!! When ``g_{i+1}`` is `nothing`, it should return the gradient of the loss function,
* `s`, the initial state ``s_0``,

Keyword arguments
* `δ`, the number of checkpoints,
* `N`, the number of time steps,
* `τ`, the number of sweeps, it is chosen as the smallest integer that `binomial(τ+δ, τ) >= N` by default,
* `f_inplace = false`, whether `f` is inplace,
* `logger = TreeverseLog()`, the logger.

Ref: https://www.tandfonline.com/doi/abs/10.1080/10556789208805505
"""
function treeverse(f, gf, s::T; δ, N, τ=binomial_fit(N,δ), f_inplace=false, logger = TreeverseLog()) where T
    state = Dict(0=>s)
    if N > binomial(τ+δ, τ)
        error("please input a larger `τ` and `δ` so that `binomial(τ+δ, τ) >= N`!")
    end
    g = treeverse!(f, gf, state, nothing, δ, τ, 0, 0, N, logger, f_inplace)
    return g
end

function treeverse!(f, gf, state::Dict{Int,T}, g, δ, τ, β, σ, ϕ, logger, f_inplace) where T
    logger.depth[] += 1
    # cache sσ
    if σ > β
        δ -= 1
        s = state[β]
        for j=β:σ-1
            s = getf(f, j)(s)
            push!(logger, :call, τ, δ, j)
        end
        store_state!(state,σ, f_inplace ? copy(s) : s)
        s = nothing
        push!(logger, :store, τ, δ, σ)
        logger.peak_mem[] = max(logger.peak_mem[], length(state))
    elseif σ < β
        error("treeverse fails! σ < β")
    end

    κ = mid(δ, τ, σ, ϕ)
    while τ>0 && κ < ϕ
        g = treeverse!(f, gf, state, g, δ, τ, σ, κ, ϕ, logger, f_inplace)
        τ -= 1
        ϕ = κ
        κ = mid(δ, τ, σ, ϕ)
    end

    g = getf(gf, σ)(state[σ], g)
    push!(logger, :grad, τ, δ, σ)
    if σ>β
        # remove state[σ]
        delete_state!(state, σ)
        push!(logger, :fetch, τ, δ, σ)
    end
    return g
end

@inline function store_state!(state::Dict, i::Int, x)
    state[i] = x
end

@inline function delete_state!(state::Dict, i::Int)
    pop!(state, i)
end

getf(f, i::Int) = f
getf(f::AbstractArray, i::Int) = f[i]


text:
module TestPkg

using FilePathsBase
using FilePathsBase: /
using FilePaths

# Test for basic definitions taking abstract paths
FilePaths.@compat function typical(x::AbstractPath, y::AbstractPath)
    return relative(x, y)
end

# Test for definition taking parameterized paths
FilePaths.@compat function parameterized(x::P, y::P) where P <: AbstractPath
    return relative(x, y)
end

# Test for definition with a path return type
# (indicating that we should convert it back to a string)
FilePaths.@compat function rtype(x::P, y::P)::P where P <: AbstractPath
    return relative(x, y)
end

# Test for a mixed input definition
# z is an unused variable to make sure untyped variables don't break the macro
FilePaths.@compat function mixed(x::P, y::String, z)::P where P <: AbstractPath
    return x / y
end

# Test kwargs
FilePaths.@compat function kwargs(x::AbstractPath; y::AbstractPath=cwd())::AbstractPath
    return relative(x, y)
end

# Test optional args
FilePaths.@compat function optargs(x::AbstractPath, y::AbstractPath=cwd())::AbstractPath
    return relative(x, y)
end

# Test for inline definitions
FilePaths.@compat inline(x::AbstractPath, y::AbstractPath) = relative(x, y)

# Test for custom path type
struct MyPath <: AbstractPath
    x::String
end

__init__() = FilePathsBase.register(MyPath)
function Base.tryparse(::Type{MyPath}, str::AbstractString)
    startswith(str, "mypath://") ? MyPath(str) : nothing
end

FilePaths.@compat mypath_testem(path::MyPath) = "**"*path.x

end  # TestPkg module

@testset "@compat" begin
    cd(absolute(parent(Path(@__FILE__)))) do
        reg = Sys.iswindows() ? "..\\src\\FilePaths.jl" : "../src/FilePaths.jl"
        @test ispath(reg)
        p = Path(reg)

        @testset "typical" begin
            # Passing paths should should work just like calling relative
            @test TestPkg.typical(p, home()) == relative(p, home())

            # Passing strings should now also work just like calling relative
            @test TestPkg.typical(string(p), string(home())) == relative(p, home())

            # Mixing strings and paths should also work here
            @test TestPkg.typical(string(p), home()) == relative(p, home())
        end

        @testset "parameterized" begin
            # Passing paths should should work just like calling relative
            @test TestPkg.parameterized(p, home()) == relative(p, home())

            # Passing strings should now also work just like calling relative
            @test TestPkg.parameterized(string(p), string(home())) == relative(p, home())

            # Mixing strings and paths should error because x and y need to be the same type
            @test_throws MethodError TestPkg.parameterized(string(p), home())
        end

        @testset "rtype" begin
            # Passing paths should should work just like calling relative
            @test TestPkg.rtype(p, home()) == relative(p, home())

            # Passing strings should also convert the output to a string
            @test TestPkg.rtype(string(p), string(home())) == string(relative(p, home()))

            # Mixing strings and paths should error because x and y need to be the same type
            @test_throws MethodError TestPkg.rtype(string(p), home())
        end

        @testset "mixed" begin
            # Passing in a path as the first argument should return the same output as p / home()
            prefix = parent(p)
            suffix = basename(p)
            @test TestPkg.mixed(prefix, suffix, 2) == p

            # Passing 2 strings should also convert the output to a string
            @test TestPkg.mixed(string(prefix), suffix, 2) == string(p)

            # Passing a path where we explicitly want a string will error
            @test_throws MethodError TestPkg.mixed(prefix, Path(suffix), 2)
        end

        @testset "optargs" begin
            # Passing paths should should work just like calling relative
            @test TestPkg.optargs(p, home()) == relative(p, home())

            # Passing strings should also return the relative path as a string
            @test TestPkg.optargs(string(p), string(home())) == string(relative(p, home()))

            # With optional arguments we should be able to mix strings on either argument
            @test TestPkg.optargs(p, string(home())) == string(relative(p, home()))
            @test TestPkg.optargs(string(p), home()) == string(relative(p, home()))

            # Use default
            @test TestPkg.optargs(p) == relative(p, cwd())
            @test TestPkg.optargs(string(p)) == string(relative(p, cwd()))
        end

        @testset "inline" begin
            # Passing paths should should work just like calling relative
            @test TestPkg.inline(p, home()) == relative(p, home())

            # Passing strings should now also work just like calling relative
            @test TestPkg.inline(string(p), string(home())) == relative(p, home())

            # Mixing strings and paths should also work here
            @test TestPkg.inline(string(p), home()) == relative(p, home())
        end

        @testset "Custom path type" begin
            @test TestPkg.mypath_testem(TestPkg.MyPath("mypath://foo")) == "**mypath://foo"
            @test TestPkg.mypath_testem("mypath://foo") == "**mypath://foo"
        end
    end
end


text:
# Operations with periodic arrays
function pmshift(A::PeriodicArray, k::Int = 1)
    return PeriodicArray(A.M[:,:,mod.(k:A.dperiod+k-1,A.dperiod).+1], A.period; nperiod = A.nperiod)
end
function LinearAlgebra.transpose(A::PeriodicArray)
    return PeriodicArray(permutedims(A.M,(2,1,3)), A.period; nperiod = A.nperiod)
end
function LinearAlgebra.adjoint(A::PeriodicArray)
    return PeriodicArray(permutedims(A.M,(2,1,3)), A.period; nperiod = A.nperiod)
end
function LinearAlgebra.norm(A::PeriodicArray, p::Real = 2)
    return norm([norm(A.M[:,:,i],p) for i in 1:A.dperiod],p)
end
function +(A::PeriodicArray, B::PeriodicArray)
    A.Ts ≈ B.Ts || error("A and B must have the same sampling time")
    period = promote_period(A, B)
    m, n, pa = size(A.M)
    mb, nb, pb = size(B.M)
    (m, n) == (mb, nb) || throw(DimensionMismatch("A and B must have the same dimensions"))
    nta = numerator(rationalize(period/A.period))
    ntb = numerator(rationalize(period/B.period))
    pa == pb && nta == 1 && ntb == 1 && (return PeriodicArray(A.M+B.M, A.period; nperiod = A.nperiod))
    K = nta*A.nperiod*pa
    p = lcm(pa,pb)
    T = promote_type(eltype(A),eltype(B))
    X = Array{T,3}(undef, m, n, p)
    for i = 1:p
        ia = mod(i-1,pa)+1
        ib = mod(i-1,pb)+1
        X[:,:,i] = A.M[:,:,ia]+B.M[:,:,ib]
    end
    return PeriodicArray(X, period; nperiod = div(K,p))
end
+(A::PeriodicArray, C::AbstractMatrix) = +(A, PeriodicArray(C, A.Ts; nperiod = 1))
+(A::AbstractMatrix, C::PeriodicArray) = +(PeriodicArray(A, C.Ts; nperiod = 1), C)
-(A::PeriodicArray) = PeriodicArray(-A.M, A.period; nperiod = A.nperiod)
-(A::PeriodicArray, B::PeriodicArray) = +(A,-B)
function *(A::PeriodicArray, B::PeriodicArray)
    A.Ts ≈ B.Ts || error("A and B must have the same sampling time")
    period = promote_period(A, B)
    m, na, pa = size(A.M)
    mb, n, pb = size(B.M)
    na == mb || throw(DimensionMismatch("A and B have incompatible dimension"))
    nta = numerator(rationalize(period/A.period))
    K = nta*A.nperiod*pa
    #K == B.nperiod*pb || error("A and B must have the same sampling time")
    #pa == pb && (return PeriodicArray(A.M*B.M, A.period; nperiod = A.nperiod))
    p = lcm(pa,pb)
    T = promote_type(eltype(A),eltype(B))
    X = Array{T,3}(undef, m, n, p)
    for i = 1:p
        ia = mod(i-1,pa)+1
        ib = mod(i-1,pb)+1
        mul!(view(X,:,:,i), view(A.M,:,:,ia), view(B.M,:,:,ib))
    end
    return PeriodicArray(X, period; nperiod = div(K,p))
end
*(A::PeriodicArray, C::AbstractMatrix) = *(A, PeriodicArray(C, A.Ts; nperiod = 1))
*(A::AbstractMatrix, C::PeriodicArray) = *(PeriodicArray(A, C.Ts; nperiod = 1), C)
*(A::PeriodicArray, C::Real) = PeriodicArray(C*A.M, A.period; nperiod = A.nperiod)
*(A::Real, C::PeriodicArray) = PeriodicArray(A*C.M, C.period; nperiod = C.nperiod)

# function -(A::PeriodicArray, B::PeriodicArray)
#     A.period ≈ B.period || error("A and B must have the same period")
#     m, n, pa = size(A.M)
#     mb, nb, pb = size(B.M)
#     (m, n) == (mb, nb) || throw(DimensionMismatch("A and B must have the same dimensions"))
#     K = A.nperiod*pa
#     K == B.nperiod*pb || error("A and B must have the same sampling time")
#     pa == pb && (return PeriodicArray(A.M-B.M, A.period; nperiod = A.nperiod))
#     p = lcm(pa,pb)
#     T = promote_type(eltype(A),eltype(B))
#     X = Array{T,3}(undef, m, n, p)
#     for i = 1:p
#         ia = mod(i-1,pa)+1
#         ib = mod(i-1,pb)+1
#         X[:,:,i] = A.M[:,:,ia]-B.M[:,:,ib]
#     end
#     return PeriodicArray(X, A.period; nperiod = div(K,p))
# end


text:
# bias/unbias
biased_exponent(x::T) where T<:Union{SysFloat, Unsigned} = Int(biased_exponent_field(x))
unbias_exponent(x::T) where T<:Union{SysFloat, Unsigned} = Int(x - exponent_bias(T))
bias_exponent(x::T) where T<:Union{SysFloat, Unsigned} = Int(x + exponent_bias(T))

# fetch the field[s] into the low order bits of an Unsigned

@inline sign_field(x::T) where T<:Unsigned = isolate_sign_field(x) >> sign_field_offset(T)
@inline exponent_field(x::T) where T<:Unsigned = isolate_exponent_field(x) >> exponent_field_offset(T)
@inline significand_field(x::T) where T<:Unsigned = isolate_significand_field(x) >> significand_field_offset(T)
@inline sign_and_exponent_fields(x::T) where T<:Unsigned = isolate_sign_and_exponent_fields(x) >> exponent_field_offset(T)
@inline exponent_and_significand_fields(x::T) where T<:Unsigned = isolate_exponent_and_significand_fields(x) >> significand_field_offset(T)

@inline biased_exponent_field(x::T) where T<:Unsigned = exponent_field(x)
@inline unbiased_exponent_field(x::T) where T<:Unsigned = exponent_field(x) - exponent_bias(T)

for F in (:sign_field, :exponent_field, :significand_field,
          :biased_exponent_field,
          :sign_and_exponent_fields, :exponent_and_significand_fields)
  @eval begin
    @inline $F(x::T) where T<:SysFloat = $F(convert(Unsigned,x))
  end
end

@inline unbiased_exponent_field(x::Float64) = convert(UInt64, exponent_field(x) - exponent_bias(Float64))
@inline unbiased_exponent_field(x::Float32) = convert(UInt32, exponent_field(x) - exponent_bias(Float32))
@inline unbiased_exponent_field(x::Float16) = convert(UInt16, exponent_field(x) - exponent_bias(Float16))

# set field[s]: sign_field(1.0, 1%UInt64) == -1.0

for (F,C,P) in ((:sign_field, :clear_sign_field, :prepare_sign_field),
                (:exponent_field, :clear_exponent_field, :prepare_exponent_field),
                (:significand_field, :clear_significand_field, :prepare_significand_field),
                (:sign_and_exponent_fields, :clear_sign_and_exponent_fields, :prepare_sign_and_exponent_fields),
                (:exponent_and_significand_fields, :clear_exponent_and_significand_fields, :prepare_exponent_and_significand_fields),
                (:biased_exponent_field, :clear_exponent_field, :prepare_exponent_field))
  for (T,U,S) in ((:Float64, :UInt64, :Int64), (:Float32, :UInt32, :Int32), (:Float16, :UInt16, :Int32))
    @eval begin
        @inline $F(x::$T, y::$U) = reinterpret($T, $C(x) | $P(y))
        @inline $F(x::$T, y::U) where U<:Unsigned = convert($T, $C(x) | $P(convert($U, y)))
    end
  end
end

@inline unbiased_exponent_field(x::U, y::U) where U<:Unsigned = clear_exponent_field(x) | prepare_unbiased_exponent_field(y)
@inline unbiased_exponent_field(x::Float64, y::UInt64) = reinterpret(Float64, unbiased_exponent_field(reinterpret(UInt64, x), y))
@inline unbiased_exponent_field(x::Float32, y::UInt32) = reinterpret(Float32, unbiased_exponent_field(reinterpret(UInt32, x), y))
@inline unbiased_exponent_field(x::Float16, y::UInt16) = reinterpret(Float16, unbiased_exponent_field(reinterpret(UInt16, x), y))

# clear the field[s] and yield the value, as Unsigned bits in place

@inline clear_sign_field(x::T) where T<:Unsigned = x & sign_field_filter(T)
@inline clear_exponent_field(x::T) where T<:Unsigned = x & exponent_field_filter(T)
@inline clear_significand_field(x::T) where T<:Unsigned = x & significand_field_filter(T)
@inline clear_sign_and_exponent_fields(x::T) where T<:Unsigned = x & sign_and_exponent_fields_filter(T)
@inline clear_exponent_and_significand_fields(x::T) where T<:Unsigned = x & exponent_and_significand_fields_filter(T)

for F in (:clear_sign_field, :clear_exponent_field, :clear_significand_field,
          :clear_sign_and_exponent_fields, :clear_exponent_and_significand_fields)
  @eval begin
    @inline $F(x::T) where T<:SysFloat = $F(convert(Unsigned,x))
  end
end

# prepare Unsigned low order bits to occupy field[s]

@inline prepare_sign_field(x::T) where T<:Unsigned = (x & sign_field_mask_lsbs(T)) << sign_field_offset(T)
@inline prepare_exponent_field(x::T) where T<:Unsigned = (x & exponent_field_mask_lsbs(T)) << exponent_field_offset(T)
@inline prepare_significand_field(x::T) where T<:Unsigned = (x & significand_field_mask_lsbs(T)) << significand_field_offset(T)
@inline prepare_sign_and_exponent_fields(x::T) where T<:Unsigned = (x & sign_and_exponent_fields_mask_lsbs(T)) << exponent_field_offset(T)
@inline prepare_exponent_and_significand_fields(x::T) where T<:Unsigned = (x & exponent_and_significand_fields_mask_lsbs(T)) << exponent_and_significand_fields_offset(T)
@inline prepare_biased_exponent_field(x::T) where T<:Unsigned = prepare_exponent_field(x)
@inline prepare_unbiased_exponent_field(x::T) where T<:Unsigned = prepare_exponent_field(x + exponent_bias(T))