text:
struct ElementFEAInfo{dim, T, TK<:AbstractMatrix{T}, Tf<:AbstractVector{T}, TKes<:AbstractVector{TK}, Tfes<:AbstractVector{Tf}, Tcload<:AbstractVector{T}, refshape, TCV<:CellValues{dim, T, refshape}, dimless1, TFV<:FaceValues{dimless1, T, refshape}}
    Kes::TKes
    fes::Tfes
    fixedload::Tcload
    cellvolumes::Tcload
    cellvalues::TCV
    facevalues::TFV
    metadata::Metadata
end
function ElementFEAInfo(sp::RectilinearPointLoad{dim, T}, quad_order=2, ::Type{Val{mat_type}}=Val{:Static}) where {dim, T, mat_type}
    Kes, cellvalues, facevalues = make_Kes(sp, quad_order, Val{mat_type})
    fes = [zeros(T, ndofs_per_cell(sp.ch.dh)) for i in 1:getncells(sp.ch.dh.grid)]
    fixedload = make_f(sp)
    cellvolumes = get_cell_volumes(sp, cellvalues)
    return ElementFEAInfo(Kes, fes, fixedload, cellvolumes, cellvalues, facevalues, sp.metadata)
end
function ElementFEAInfo(sp::InpLinearElasticity, quad_order=2, ::Type{Val{mat_type}}=Val{:Static}) where {mat_type}
    Kes, weights, dloads, cellvalues, facevalues = make_Kes_and_fes(sp, quad_order, Val{mat_type})
    fixedload = full(make_cload(sp))
    assemble_f!(fixedload, sp, dloads)
    cellvolumes = get_cell_volumes(sp, cellvalues)
    ElementFEAInfo(Kes, weights, fixedload, cellvolumes, cellvalues, facevalues, sp.metadata)
end
function ElementFEAInfo(Kes::AbstractVector{TK}, fes::AbstractVector{Tf},fixedload::AbstractVector{T}, cellvolumes::AbstractVector{T}, cellvalues::TCV, facevalues::TFV, metadata::Metadata) where {T, TK<:AbstractMatrix{T}, Tf<:AbstractVector{T}, dim, dimless1, refshape, TCV<:CellValues{dim, T, refshape}, TFV<:FaceValues{dimless1, T, refshape}}
    ElementFEAInfo{dim, T, TK, Tf, typeof(Kes), typeof(fes), typeof(fixedload), refshape, TCV, dimless1, TFV}(Kes, fes, fixedload, cellvolumes, cellvalues, facevalues, metadata)
end

function get_cell_volumes(sp::LinearElasticityProblem{dim, T}, cellvalues) where {dim, T}
    dh = sp.ch.dh
    cellvolumes = zeros(T, getncells(dh.grid))
    for (i, cell) in enumerate(CellIterator(dh))
        reinit!(cellvalues, cell)
        cellvolumes[i] = sum(JuAFEM.getdetJdV(cellvalues, q_point) for q_point in 1:JuAFEM.getnquadpoints(cellvalues))
    end
    return cellvolumes
end

mutable struct GlobalFEAInfo{T, TK<:AbstractMatrix{T}, Tf<:AbstractVector{T}}
    K::TK
    f::Tf
end
GlobalFEAInfo(K::AbstractMatrix{T}, f::AbstractVector{T}) where T = GlobalFEAInfo{T, typeof(K), typeof(f)}(K, f)
GlobalFEAInfo(::Type{T}) where T = GlobalFEAInfo{T}()
GlobalFEAInfo() = GlobalFEAInfo{Float64}()
GlobalFEAInfo{T}() where T = GlobalFEAInfo{T, SparseMatrixCSC{T, Int}, Vector{T}}(sparse(zeros(T, 0, 0)), zeros(T, 0))
GlobalFEAInfo(sp::LinearElasticityProblem) = GlobalFEAInfo(make_empty_K(sp), make_empty_f(sp))
make_empty_K(sp::LinearElasticityProblem) = Symmetric(create_sparsity_pattern(sp.ch.dh))
make_empty_f(sp::LinearElasticityProblem{dim, T}) where {dim, T} = zeros(T, ndofs(sp.ch.dh))

function make_Kes(sp::RectilinearPointLoad, quad_order=2)
    make_Kes(sp, quad_order, Val{:Static})
end
function make_Kes(sp::RectilinearPointLoad, ::Type{Val{mat_type}}) where mat_type
    make_Kes(sp, 2, Val{mat_type})
end
function make_Kes(sp::RectilinearPointLoad{dim, T}, quad_order, ::Type{Val{mat_type}}) where {dim, T, mat_type}
    E = sp.E
    ν = sp.ν

    λ = E*ν / ((1 + ν) * (1 - 2*ν))
    μ = E / (2*(1 + ν))
    δ(i,j) = i == j ? T(1) : T(0)
    g(i,j,k,l) = λ*δ(i,j)*δ(k,l) + μ*(δ(i,k)*δ(j,l) + δ(i,l)*δ(j,k))
    C = SymmetricTensor{4, dim}(g)

    dh = sp.ch.dh

    # Shape functions and quadrature rule
    interpolation_space = Lagrange{dim, RefCube, 1}()
    quadrature_rule = QuadratureRule{dim, RefCube}(quad_order)
    cellvalues = CellScalarValues(quadrature_rule, interpolation_space)
    facevalues = FaceScalarValues(QuadratureRule{dim-1, RefCube}(quad_order), interpolation_space)

    # Calculate element stiffness matrices
    n_basefuncs = getnbasefunctions(cellvalues)

    return _make_Kes(sp, Val{mat_type}, Val{n_basefuncs}, Val{dim*n_basefuncs}, C, interpolation_space, quadrature_rule, cellvalues, facevalues), cellvalues, facevalues
end
function _make_Kes(sp::RectilinearPointLoad{dim, T}, ::Type{Val{mat_type}}, ::Type{Val{n_basefuncs}}, ::Type{Val{ndofs_per_cell}}, C, interpolation_space, quadrature_rule, cellvalues, facevalues) where {dim, T, mat_type, n_basefuncs, ndofs_per_cell}
    dh = sp.ch.dh

    # Calculate element stiffness matrices
    Kesize = ndofs_per_cell
    nel = getncells(dh.grid)

    if mat_type === :Static || mat_type === :SMatrix
        if !(T === BigFloat)
            Kes = Symmetric{T, SMatrix{Kesize, Kesize, T, Kesize^2}}[]
        else
            Kes = Symmetric{T, SizedMatrix{Kesize, Kesize, T, Kesize^2}}[]
        end
        sizehint!(Kes, nel)

        Ke_e = zeros(T, dim, dim)
        Ke_0 = Matrix{T}(Kesize, Kesize)
        celliteratortype = CellIterator{typeof(dh).parameters...}
        _celliterator::celliteratortype = CellIterator(dh)
        for (k, cell) in enumerate(_celliterator)
            Ke_0 .= 0
            reinit!(cellvalues, cell)
            for q_point in 1:getnquadpoints(cellvalues)
                dΩ = getdetJdV(cellvalues, q_point)
                for b in 1:n_basefuncs
                    ∇ϕb = shape_gradient(cellvalues, q_point, b)
                    for d2 in 1:dim
                        for a in 1:n_basefuncs
                            ∇ϕa = shape_gradient(cellvalues, q_point, a)
                            Ke_e .= dotdot(∇ϕa, C, ∇ϕb) * dΩ
                            for d1 in 1:dim
                                #if dim*(b-1) + d2 >= dim*(a-1) + d1
                                Ke_0[dim*(a-1) + d1, dim*(b-1) + d2] += Ke_e[d1,d2]
                                #end
                            end
                        end
                    end
                end
            end
            if !(T === BigFloat)
                push!(Kes, Symmetric(SMatrix{Kesize, Kesize, T, Kesize*Kesize}(Ke_0)))
            else
                push!(Kes, Symmetric(SizedMatrix{Kesize, Kesize, T}(Ke_0)))
            end
        end
    else
        Kes = let Kesize=Kesize, nel=nel
            [Symmetric(zeros(T, Kesize, Kesize), :U) for i = 1:nel]
        end

        Ke_e = zeros(T, dim, dim)

        celliteratortype = CellIterator{typeof(dh).parameters...}
        _celliterator = CellIterator(dh)
        for (k, cell) in enumerate(_celliterator)
            reinit!(cellvalues, cell)
            for q_point in 1:getnquadpoints(cellvalues)
                dΩ = getdetJdV(cellvalues, q_point)
                for b in 1:n_basefuncs
                    ∇ϕb = shape_gradient(cellvalues, q_point, b)
                    for d2 in 1:dim
                        for a in 1:n_basefuncs
                            ∇ϕa = shape_gradient(cellvalues, q_point, a)
                            Ke_e .= dotdot(∇ϕa, C, ∇ϕb) * dΩ
                            for d1 in 1:dim
                                #if dim*(b-1) + d2 >= dim*(a-1) + d1
                                Kes[k].data[dim*(a-1) + d1, dim*(b-1) + d2] += Ke_e[d1,d2]
                                #end
                            end
                        end
                    end
                end
            end
        end
    end
    return Kes
end

function make_f(sp::RectilinearPointLoad)
    dof_cells = sp.metadata.dof_cells
    fdof = sp.force_dof
    force = sp.force
    if fdof isa Vector
        f = sparsevec(fdof, [-force for i in 1:length(fdof)], ndofs(sp.ch.dh))
    else
        f = sparsevec([fdof], [-force], ndofs(sp.ch.dh))
    end
    return f
end

function make_Kes_and_fes(problem::InpLinearElasticity, quad_order=2)
    make_Kes_and_fes(problem, quad_order, Val{:Static})
end
function make_Kes_and_fes(problem::InpLinearElasticity, ::Type{Val{mat_type}}) where mat_type
    make_Kes_and_fes(problem, 2, Val{mat_type})
end
function make_Kes_and_fes(inp_problem::InpLinearElasticity{dim, N, T, M, TI, GO}, quad_order, ::Type{Val{mat_type}}) where {dim, N, T, M, TI, GO, mat_type}
    problem = inp_problem.inp_content
    dh = inp_problem.ch.dh
    E = problem.E
    ν = problem.ν
    ρ = problem.density

    refshape = JuAFEM.getrefshape(dh.field_interpolations[1])
    geom_order = GO

    λ = E*ν / ((1 + ν) * (1 - 2*ν))
    μ = E / (2*(1 + ν))
    δ(i,j) = i == j ? T(1) : T(0)
    g(i,j,k,l) = λ*δ(i,j)*δ(k,l) + μ*(δ(i,k)*δ(j,l) + δ(i,l)*δ(j,k))
    C = SymmetricTensor{4, dim}(g)

    # Shape functions and quadrature rule
    interpolation_space = Lagrange{dim, refshape, geom_order}()
    quadrature_rule = QuadratureRule{dim, refshape}(quad_order)
    cellvalues = CellScalarValues(quadrature_rule, interpolation_space)
    facevalues = FaceScalarValues(QuadratureRule{dim-1, refshape}(quad_order), interpolation_space)

    # Calculate element stiffness matrices
    n_basefuncs = getnbasefunctions(cellvalues)

    Kes, weights = _make_Kes_and_weights(dh, Val{mat_type}, Val{n_basefuncs}, Val{dim*n_basefuncs}, C, ρ, quadrature_rule, cellvalues)
    dloads = _make_dloads(weights, inp_problem, facevalues)

    return Kes, weights, dloads, cellvalues, facevalues
end

const g = [0., 9.81, 0.] # N/kg or m/s^2

function _make_Kes_and_weights(dh::DofHandler{dim, N, T}, ::Type{Val{mat_type}}, ::Type{Val{n_basefuncs}}, ::Type{Val{ndofs_per_cell}}, C, ρ, quadrature_rule, cellvalues) where {dim, N, T, mat_type, n_basefuncs, ndofs_per_cell}
    # Calculate element stiffness matrices
    Kesize = ndofs_per_cell
    nel = getncells(dh.grid)
    body_force = ρ .* g # Force per unit volume

    if !(T === BigFloat)
        if mat_type === :Static || mat_type === :SMatrix
            MatrixType = SMatrix{Kesize, Kesize, T, Kesize^2}
            VectorType = MVector{Kesize, T}
        elseif mat_type === :MMatrix
            MatrixType = MMatrix{Kesize, Kesize, T, Kesize^2}
            VectorType = MVector{Kesize, T}
        else
            MatrixType = Matrix{T}
            VectorType = Vector{T}
        end
    else
        if mat_type === :Static || mat_type === :SMatrix  || mat_type === :MMatrix
            MatrixType = SizedMatrix{Kesize, Kesize, T, Kesize^2}
            VectorType = SizedVector{Kesize, T}
        else
            MatrixType = Matrix{T}
            VectorType = Vector{T}
        end
    end

    if MatrixType <: StaticArray
        Kes = Symmetric{T, MatrixType}[]
        sizehint!(Kes, nel)
        weights = [zeros(VectorType) for i in 1:nel]

        Ke_e = zeros(T, dim, dim)
        fe = zeros(T, Kesize)
        Ke_0 = Matrix{T}(Kesize, Kesize)
        celliteratortype = CellIterator{typeof(dh).parameters...}
        _celliterator::celliteratortype = CellIterator(dh)
        for (k, cell) in enumerate(_celliterator)
            Ke_0 .= 0
            try
                reinit!(cellvalues, cell)
            catch
                @show cell.coords, cell.nodes, cell.current_cellid[]
            end
            fe = weights[k]
            for q_point in 1:getnquadpoints(cellvalues)
                dΩ = getdetJdV(cellvalues, q_point)
                for b in 1:n_basefuncs
                    ∇ϕb = shape_gradient(cellvalues, q_point, b)
                    ϕb = shape_value(cellvalues, q_point, b)
                    for d2 in 1:dim
                        fe[(b-1)*dim + d2] += ϕb * body_force[d2] * dΩ
                        for a in 1:n_basefuncs
                            ∇ϕa = shape_gradient(cellvalues, q_point, a)
                            Ke_e .= dotdot(∇ϕa, C, ∇ϕb) * dΩ
                            for d1 in 1:dim
                                #if dim*(b-1) + d2 >= dim*(a-1) + d1
                                Ke_0[dim*(a-1) + d1, dim*(b-1) + d2] += Ke_e[d1,d2]
                                #end
                            end
                        end
                    end
                end
            end
            if MatrixType <: SizedMatrix # Work around because full constructor errors
                push!(Kes, Symmetric(SizedMatrix{Kesize,Kesize,T}(Ke_0)))
            else
                push!(Kes, Symmetric(MatrixType(Ke_0)))
            end
        end
    else
        Kes = let Kesize=Kesize, nel=nel
            [Symmetric(zeros(T, Kesize, Kesize), :U) for i = 1:nel]
        end
        weights = let Kesize=Kesize, nel=nel
            [zeros(T, Kesize) for i = 1:nel]
        end

        Ke_e = zeros(T, dim, dim)

        celliteratortype = CellIterator{typeof(dh).parameters...}
        _celliterator = CellIterator(dh)
        for (k, cell) in enumerate(_celliterator)
            reinit!(cellvalues, cell)
            fe = weights[k]
            for q_point in 1:getnquadpoints(cellvalues)
                dΩ = getdetJdV(cellvalues, q_point)
                for b in 1:n_basefuncs
                    ∇ϕb = shape_gradient(cellvalues, q_point, b)
                    ϕb = shape_value(cellvalues, q_point, b)
                    for d2 in 1:dim
                        # Force business
                        fe[(b-1)*dim + d2] += ϕb * body_force[d2] * dΩ
                        for a in 1:n_basefuncs
                            ∇ϕa = shape_gradient(cellvalues, q_point, a)
                            Ke_e .= dotdot(∇ϕa, C, ∇ϕb) * dΩ
                            for d1 in 1:dim
                                #if dim*(b-1) + d2 >= dim*(a-1) + d1
                                Kes[k].data[dim*(a-1) + d1, dim*(b-1) + d2] += Ke_e[d1,d2]
                                #end
                            end
                        end
                    end
                end
            end
        end
    end
    return Kes, weights
end

function _make_dloads(fes, inp_problem::InpLinearElasticity{dim, N, T}, facevalues) where {dim, N, T}
    dloads = deepcopy(fes)
    for i in 1:length(dloads)
        dloads[i] .= 0
    end
    problem = inp_problem.inp_content
    grid = inp_problem.ch.dh.grid
    boundary_matrix = grid.boundary_matrix
    cell_coords = zeros(JuAFEM.Vec{dim, T}, N)
    n_basefuncs = getnbasefunctions(facevalues)
    for k in keys(problem.dloads)
        t = -problem.dloads[k] # traction = negative the pressure
        faceset = problem.facesets[k]
        for (cellid, faceid) in faceset
            boundary_matrix[faceid, cellid] || throw("Face $((cellid, faceid)) not on boundary.")
            fe = dloads[cellid]
            getcoordinates!(cell_coords, grid, cellid)
            reinit!(facevalues, cell_coords, faceid)
            for q_point in 1:getnquadpoints(facevalues)
                dΓ = getdetJdV(facevalues, q_point) # Face area
                normal = getnormal(facevalues, q_point) # Nomral vector at quad point
                for i in 1:n_basefuncs
                    ϕ = shape_value(facevalues, q_point, i) # Shape function value
                    for d = 1:dim
                        fe[(i-1)*dim + d] += ϕ * t * normal[d] * dΓ
                    end
                end
            end
        end
    end

    return dloads
end

function make_cload(inp_problem::InpLinearElasticity{dim, N, T}) where {dim, N, T}
    cloads = inp_problem.inp_content.cloads
    dh = inp_problem.ch.dh
    node_dofs = inp_problem.metadata.node_dofs
    inds = Int[]
    vals = T[]
    for nodeidx in keys(cloads)
        for (dofidx, force) in enumerate(cloads[nodeidx])
            if force != 0
                dof = node_dofs[(nodeidx-1)*dim+dofidx]
                push!(inds, dof)
                push!(vals, force)
            end
        end
    end
    return sparsevec(inds, vals, ndofs(dh))
end


text:
using Test
using Statistics
using LinearAlgebra
using Random

using EnhancedGJK
using EnhancedGJK: projection_weights, projection_weights_reference, reset!, normal
using CoordinateTransformations: IdentityTransformation, Translation
using StaticArrays: SVector
import GeometryTypes
const gt = GeometryTypes
using FileIO

const mesh_dir = joinpath(dirname(@__FILE__), "meshes")

@testset "Issue #17" begin
    # Note: when this test was written, it mattered whether it was the first
    # test or not.
    geomA = gt.FlexibleConvexHull(gt.Point{2,Float64}[
        gt.Point(1.125, 2.0),
        gt.Point(0.125, 2.0),
        gt.Point(0.125, 2.8),
        gt.Point(0.325, 3.0),
        gt.Point(0.925, 3.0),
        gt.Point(1.125, 2.8)]
    )
    geomB = gt.FlexibleConvexHull(gt.Point{2,Float64}[
        gt.Point(-0.025, 1.2321067811865476),
        gt.Point(-0.025, 1.2821067811865474),
        gt.Point(0.3267766952966369, 1.2821067811865474),
        gt.Point(0.3267766952966369, 1.2321067811865476)]
    )

    cache = CollisionCache(geomA, geomB)
    poseA = poseB = IdentityTransformation()
    simplex = EnhancedGJK.transform_simplex(cache, poseA, poseB)
    @test isapprox(projection_weights(simplex), projection_weights_reference(simplex))

    result = gjk(geomA, geomB)
    @test !result.in_collision
end

@testset "reference distance" begin
    mesh = load(joinpath(mesh_dir, "base_link.obj"))
    for x in range(0.05, stop=1, length=10)
        for y in range(-1, stop=1, length=10)
            for z in range(-1, stop=1, length=10)
                point = SVector(x, y, z)
                result = gjk(mesh, point)
                simple_dist = ReferenceDistance.signed_distance(mesh, point)
                if result.in_collision
                    @test isapprox(simplex_penetration_distance(result), -simple_dist, atol=1e-3)
                else
                    @test isapprox(separation_distance(result), simple_dist, atol=1e-3)
                end
            end
        end
    end
end

@testset "reference interior distance" begin
    verts = [gt.Point(1., -1, 1), gt.Point(1., 1, 1), gt.Point(-1., 1, 1), gt.Point(-1., -1, 1),
             gt.Point(1., -1, -1), gt.Point(1., 1, -1), gt.Point(-1., 1, -1), gt.Point(-1., -1, -1)]
    ft = gt.Face{3, Int}
    faces = [ft(1, 2, 5), ft(6, 5, 2),
             ft(2, 3, 6), ft(7, 6, 3),
             ft(3, 4, 7), ft(8, 7, 4),
             ft(4, 1, 8), ft(5, 8, 1),
             ft(2, 1, 3), ft(4, 3, 1),
             ft(5, 6, 7), ft(7, 8, 5)]
    mesh = gt.HomogenousMesh(verts, faces)
    ReferenceDistance.signed_distance(mesh, SVector(0, 0, 0))

    for x in range(-1, stop=1, length=10)
        for y in range(-1, stop=1, length=10)
            for z in range(-1, stop=1, length=10)
                point = SVector(x, y, z)
                expected = -min(1 - abs(x), 1 - abs(y), 1 - abs(z))
                actual = ReferenceDistance.signed_distance(mesh, point)
                @test isapprox(expected, actual, atol=1e-15)
            end
        end
    end
end

@testset "table" begin
    width = 0.5
    thickness = 0.05
    surface_points = Vector{SVector{3, Float64}}()
    for z in [-thickness, thickness]
    for x in [-width, width]
        for y in [-width, width]
            push!(surface_points, SVector(x, y, z))
        end
    end
    end
    geometry = SVector{length(surface_points)}(surface_points)
    point = zeros(SVector{3, Float64})
    for x in range(-width, stop=width, length=21)
        for y in range(-width, stop=width, length=21)
            z = 0.1
            let result = gjk(geometry, point, IdentityTransformation(), Translation(SVector(x, y, z)))
                @test separation_distance(result) ≈ 0.05
            end
            z = 0.06
            let result = gjk(geometry, point, IdentityTransformation(), Translation(SVector(x, y, z)))
                @test separation_distance(result) ≈ 0.01 atol=1e-12
            end
            z = 0.05
            let result = gjk(geometry, point, IdentityTransformation(), Translation(SVector(x, y, z)))
                if result.in_collision
                    @test simplex_penetration_distance(result) ≈ 0.0 atol=1e-12
                else
                    @test separation_distance(result) ≈ 0.0 atol=1e-12
                end
            end
        end
    end
end

@testset "johnson distance subalgorithm" begin
    include("johnson_distance.jl")
end

@testset "simplex distance" begin
    simplex = SVector{3}(SVector{2, Float64}[[1., 0], [2., 0], [1., 1]])
    pt = SVector(0., 0)
    cache = CollisionCache(simplex, pt);
    result = gjk!(cache, IdentityTransformation(), IdentityTransformation())
    @test isapprox(separation_distance(result), 1.0)
    @test isapprox(result.closest_point_in_body.a, [1.0, 0.0])
    @test isapprox(result.closest_point_in_body.b, [0.0, 0.0])
end

@testset "mesh to mesh" begin
    mesh = load(joinpath(mesh_dir, "r_foot_chull.obj"))
    dx = 1.0
    foot_length = 0.172786 + 0.090933

    cache = CollisionCache(mesh, mesh)
    result = gjk!(cache, IdentityTransformation(), Translation(SVector(dx, 0, 0)))
    @test isapprox(separation_distance(result), dx - foot_length, atol=1e-3)

    cache = CollisionCache(mesh, mesh)
    result = gjk!(cache, Translation(SVector(dx, 0, 0)), IdentityTransformation())
    @test isapprox(separation_distance(result), dx - foot_length, atol=1e-3)

    cache = CollisionCache(mesh, mesh)
    expected_penetration = 0.01
    result = gjk!(cache, IdentityTransformation(), Translation(foot_length - expected_penetration, 0, 0))
    # TODO: penetration distance is inconsistent and inaccurate
    @test simplex_penetration_distance(result) > 0
end

@testset "neighbor mesh to mesh" begin
    mesh = NeighborMesh(load(joinpath(mesh_dir, "r_foot_chull.obj")))
    dx = 1.0
    foot_length = 0.172786 + 0.090933

    cache = CollisionCache(mesh, mesh)
    result = gjk!(cache, IdentityTransformation(), Translation(SVector(dx, 0, 0)))
    @test isapprox(separation_distance(result), dx - foot_length, atol=1e-3)

    cache = CollisionCache(mesh, mesh)
    result = gjk!(cache, Translation(SVector(dx, 0, 0)), IdentityTransformation())
    @test isapprox(separation_distance(result), dx - foot_length, atol=1e-3)

    cache = CollisionCache(mesh, mesh)
    expected_penetration = 0.01
    result = gjk!(cache, IdentityTransformation(), Translation(foot_length - expected_penetration, 0, 0))
    # TODO: penetration distance is inconsistent and inaccurate
    @test simplex_penetration_distance(result) > 0
end

@testset "geometry types" begin
    # Adapted from
    # https://github.com/JuliaGeometry/GeometryTypes.jl/blob/master/test/gjk.jl
    @testset "gjk examples" begin
        c1 = gt.Simplex(gt.Vec(-1.))
        c2 = gt.Simplex(gt.Vec(4.))
        cache = CollisionCache(c1, c2)
        result = gjk!(cache, IdentityTransformation(), IdentityTransformation())
        @test isapprox(separation_distance(result), 5.0)

        c1 = gt.Simplex(gt.Vec(-1.,0,0))
        c2 = gt.Simplex(gt.Vec(4.,0,0))
        cache = CollisionCache(c1, c2)
        result = gjk!(cache, IdentityTransformation(), IdentityTransformation())
        @test isapprox(separation_distance(result), 5.0)

        c1 = gt.FlexibleConvexHull([gt.Vec(0.,0), gt.Vec(0.,1), gt.Vec(1.,0),gt.Vec(1.,1)])
        c2 = gt.Simplex(gt.Vec(4.,0.5))
        cache = CollisionCache(c1, c2)
        result = gjk!(cache, IdentityTransformation(), IdentityTransformation())
        @test isapprox(separation_distance(result), 3.0)

        pt1 = gt.Vec(1,2,3.)
        pt2 = gt.Vec(3,4,5.)
        cache = CollisionCache(pt1, pt2)
        result = gjk!(cache, IdentityTransformation(), IdentityTransformation())
        @test isapprox(separation_distance(result), norm(pt1 - pt2))

        pt1 = gt.Point(1,2,3.)
        pt2 = gt.Point(3,4,5.)
        cache = CollisionCache(pt1, pt2)
        result = gjk!(cache, IdentityTransformation(), IdentityTransformation())
        @test isapprox(separation_distance(result), norm(pt1 - pt2))
    end

    @testset "gjk intersecting lines" begin
        c1 = gt.Simplex(gt.Vec(1,1.), gt.Vec(1, 2.))
        @test simplex_penetration_distance(gjk(c1, c1)) == 0.

        c2 = gt.Simplex(gt.Vec(1,1.), gt.Vec(10, 2.))
        @test simplex_penetration_distance(gjk(c1, c2)) == 0.
    end
end

@testset "reset!" begin
    c1 = gt.FlexibleConvexHull([gt.Vec(0.,0), gt.Vec(0.,1), gt.Vec(1.,0),gt.Vec(1.,1)])
    c2 = gt.Simplex(gt.Vec(4.,0.5))
    cache = CollisionCache(c1, c2)
    initial_simplex = deepcopy(cache.simplex_points)
    gjk!(cache, IdentityTransformation(), IdentityTransformation())
    @test cache.simplex_points != initial_simplex
    reset!(cache)
    @test cache.simplex_points == initial_simplex
    allocs = @allocated reset!(cache)
    @test allocs == 0
end

@testset "normal, N = $N" for N = 2 : 4
    rng = MersenneTwister(1)
    for j = 1 : 100
        face = SVector{N}(ntuple(i -> rand(rng, SVector{N}), Val(N)))
        n = normal(face)
        dot_products = [point ⋅ n for point in face]
        @test all(x -> isapprox(x, 0; atol=1e-10), dot_products .- mean(dot_products))
    end
end

@testset "benchmarks" begin
    include("../perf/runbenchmarks.jl")
end

ext: jl


text:
import Base: show
import Compat.Distributed: AbstractWorkerPool

# abstracts
export AbstractSampler, SamplerState

export SamplePlan

# Interface
export state, update!, sample!

"""
Abstract type for all samplers

## Interface

- `show(io::IO, sampler)`: print verbose information for a sampler
- `algorithm(sampler) -> String`: print algorithm name
- `state(sampler) -> SamplerState`: get sampler state
"""
abstract type AbstractSampler end

"""
Abstract type for all sampler state
"""
abstract type SamplerState end

"""
Sample Plan contains a distribution and a sampler.

## Interface

- `state(plan::SamplePlan) -> SamplerState`: get sampler's state
- `update!(plan::SamplePlan) -> SamplePlan`: update sample plan according to the sampler
- `show(plan::SamplePlan)`: print verbose information

## Note

using a list of distribution (probability) is also allowed. This means you have
to define your list type as a subtype of `AbstractVector` and it has to support
`getindex` method at least.

using a `Function` as your distribution is more common. And you can define your
custom type of distribution function by define a subtype of `Function` with its
`call` method (AKA `function (f::YourFunctionType)(x)`).
"""
struct SamplePlan{S<:AbstractSampler, D}
    distribution::D
    sampler::S
end

"""
    state(plan)

get sampler's state in `SamplePlan`
"""
state(plan::SamplePlan) = state(plan.sampler)

"""
    update!(plan)

update plan accroding to certain sampler
"""
function update!(plan::SamplePlan) end

"""
    sample!(plan) -> data

execute your sample plan
"""
function sample!(plan::SamplePlan) end

"""
    sample!(pool, plan; nchain=4) -> data

execute your sample plan according to a `WorkPool`
"""
sample!(pool::AbstractWorkerPool, plan::SamplePlan; nchain=4) =
    pmap(pool, x->sample!(plan), 1:nchain)

"""
    sample!(nchain, plan) -> data

execute your sample plan in parallel, with `n` chains.
"""
sample!(nchain::Int, plan::SamplePlan) =
    pmap(x->sample!(plan), 1:nchain)

function show(io::IO, plan::SamplePlan{S, D}) where {S<:AbstractSampler, D<:Function}
    println(io, "Sample Plan:")
    println("Algorithm: $(algorithm(plan.sampler))")
    print("Distribution Type: Function Distribution")
end

function show(io::IO, plan::SamplePlan{S, D}) where {S<:AbstractSampler, D<:AbstractVector}
    println(io, "Sample Plan:")
    println("Algorithm: $(algorithm(plan.sampler))")
    print("Distribution Type: Probability List")
end

### sugars
export @sample

"""
    sample(expr)

a sugar macro to construct a `SamplePlan`
"""
macro sample(expr)
    warn("Not Implemented")
end


text:
cd(@__DIR__)
using Pkg; Pkg.activate("."); Pkg.instantiate()

using Printf

using Oceananigans
using Oceananigans.Diagnostics
using Oceananigans.OutputWriters

# Model parameters
N = 128  # Number of grid points in each dimension.
L = 1   # Physical length of each dimension.

κ = 0.05  # Diffusivity
ν = κ     # Assuming Prandtl number Pr = 1

end_time = 1.5

# Boundary conditions
c_bcs = HorizontallyPeriodicBCs(   top = BoundaryCondition(Neumann, 0),
                                bottom = BoundaryCondition(Neumann, 0))
model_bcs = HorizontallyPeriodicSolutionBCs(c=c_bcs)

# Set up forcing
@inline Fc(i, j, k, grid, t, U, C, params) = @inbounds cos(sin(C.c[i, j, k]^3)) + sin(cos(C.c[i, j, k]^2))
forcing = ModelForcing(c=Fc)

# Set up model
model = Model(
           architecture = CPU(),
             float_type = Float64,
                   grid = RegularCartesianGrid(size=(N, N, N), x=(0, L), y=(0, L), z=(0, L)),
                tracers = (:c,),
               coriolis = nothing,
               buoyancy = nothing,
                closure = ConstantIsotropicDiffusivity(ν=ν, κ=κ),
    boundary_conditions = model_bcs,
                forcing = forcing
)

# Setting initial conditions
c₀(x, y, z) = exp(-200(z-0.75)^2) + 1e-8 * randn()
set!(model, c=c₀)

# Add JLD2 output writer for 3D fields.
fields = Dict(
    :u => model -> Array(interior(model.velocities.u)),
    :v => model -> Array(interior(model.velocities.v)),
    :w => model -> Array(interior(model.velocities.w)),
    :c => model -> Array(interior(model.tracers.c))
)

model.output_writers[:fields] =
    JLD2OutputWriter(model, fields, dir=".", prefix="advection_diffusion_3d_fields",
                     interval=0.1, force=true, verbose=true)

# Add JLD2 output writer for horizontal averages.
 u̅ = HorizontalAverage(model.velocities.u; return_type=Array)
 v̅ = HorizontalAverage(model.velocities.v; return_type=Array)
 w̅ = HorizontalAverage(model.velocities.w; return_type=Array)
 c̅ = HorizontalAverage(model.tracers.c; return_type=Array)

 profiles = Dict(
    :u => model -> u̅(model),
    :v => model -> v̅(model),
    :w => model -> w̅(model),
    :c => model -> c̅(model)
)

model.output_writers[:horizontal_averages] =
    JLD2OutputWriter(model, profiles, dir=".", prefix="advection_diffusion_3d_horizontal_averages",
                     interval=0.01, force=true, verbose=true)

# Set up adaptive time stepping.
wizard = TimeStepWizard(cfl=0.1, Δt=1e-6, max_change=1.2, max_Δt=1e-1)

# Set up CFL diagnostics.
 cfl = AdvectiveCFL(wizard)

# Number of time steps to perform at a time before printing a progress
# statement and updating the adaptive time step.
Ni = 20

while model.clock.time < end_time
    walltime = @elapsed time_step!(model; Nt=Ni, Δt=wizard.Δt)

    # Calculate simulation progress in %.
    progress = 100 * (model.clock.time / end_time)

    # Calculate maximum velocities.
    umax = maximum(abs, model.velocities.u.data.parent)
    vmax = maximum(abs, model.velocities.v.data.parent)
    wmax = maximum(abs, model.velocities.w.data.parent)

    # Calculate a new adaptive time step.
    update_Δt!(wizard, model)

    # Print progress statement.
    i, t = model.clock.iteration, model.clock.time
    @printf("[%06.2f%%] i: %d, t: %5.2e, umax: (%6.3e, %6.3e, %6.3e), CFL: %6.4e, next Δt: %8.5e, ⟨wall time⟩: %s\n",
            progress, i, t, umax, vmax, wmax, cfl(model), wizard.Δt, prettytime(walltime / Ni))
end


text:
@testset "FieldFactory" begin

  println("Quadratic Fields")
  @time begin
    @testset "Quadratic Fields" begin
      lq = fields(2, 1, fmpz(10)^4)
      @test length(lq) == 6086

      lqreal = fields(2, 1, fmpz(10)^4, only_real = true)
      @test length(lqreal) == 3043

      ind = 0
      for x in lq
        if istotally_real(x.field)
          ind += 1
        end
      end
      @test ind == 3043
    end
  end

  println("Biquadratic Fields")
  @time begin
    @testset "Biquadratic Fields" begin
      l = fields(4, 2, fmpz(10)^6)
      @test length(l) ==  1014

      l1 = fields(4, 2, fmpz(10)^6, only_real = true)
      @test length(l1) == 196

      ind = 0
      for x in l
        if istotally_real(x.field)
          ind += 1
        end
      end
      @test ind == 196
    end
  end

  println("Abelian Fields")
  @time begin
    @testset "Abelian fields" begin
      l = fields(6, 2, fmpz(10)^6)
      @test length(l) ==  10

      l1 = fields(9, 2, fmpz(10)^13)
      @test length(l1) == 4
    end
  end

  println("Direct product decomposition")
  @time begin
    @testset "Direct product decomposition" begin
      @time l_direct_product = fields(12, 4, fmpz(10)^13)
      @test length(l_direct_product) == 25

      @time l_without = fields(12, 4, fmpz(10)^13, using_direct_product = false)
      @test length(l_without) == 25

      d1 = length(fields(24, 13, fmpz(10)^24))
      d2 = length(fields(24, 13, fmpz(10)^24, using_direct_product = false))
      @test d1 == d2

      @time l_direct_product = fields(18, 3, fmpz(10)^18)
      @time l_without = fields(18, 3, fmpz(10)^18, using_direct_product = false)
      @test length(l_direct_product) == length(l_without)
      for x in l_direct_product
        @test GAP.gap_to_julia(Vector{Int}, Hecke.IdGroup(closure(x.generators_of_automorphisms))) == [19, 3]
      end
      for x in l_without
        @test GAP.gap_to_julia(Vector{Int}, Hecke.IdGroup(closure(x.generators_of_automorphisms))) == [18, 3]
      end
    end
  end

  println("Split extension")
  @time begin
    @testset "Split extension" begin
      x = PolynomialRing(FlintZZ, "x", cached = false)[2]
      K =  number_field(x^2+1)[1]
      l = Hecke.FieldsTower[Hecke.field_context(K)]
      l1 = fields(42, 5, l, fmpz(10)^90)
      @test length(l1) == 1
      @test GAP.gap_to_julia(Vector{Int}, Hecke.IdGroup(closure(l1[1].generators_of_automorphisms))) == [42, 5]
    end
  end

  println("Obstructions: degree 2")
  @time begin
    @testset "Obstructions: degree 2" begin
      G = GAP.Globals.SmallGroup(8, 4)
      L = GAP.Globals.DerivedSeries(G)
      l = fields(4, 2, fmpz(10)^6)
      forQ8 = Hecke.check_obstruction(l, L, 2, [2])
      @test length(forQ8) == 53

      lQ8 = fields(8, 4, fmpz(10)^12)
      @test length(lQ8) == 2

      lQ8_2 = fields(8, 4, forQ8, fmpz(10)^12)
      @test length(lQ8_2) == 2

      lQ8real = fields(8, 4, forQ8, fmpz(10)^12, only_real = true)
      @test length(lQ8real) == 1
    end
  end

  println("Obstructions: prime degree")
  @time begin
    @testset "Obstructions: prime degree" begin
      l = fields(6, 2, fmpz(10)^8)
      G = GAP.Globals.SmallGroup(54, 5)
      L = GAP.Globals.DerivedSeries(G)
      lsieved = Hecke.check_obstruction(l, L, 2, [3, 3])
      @test length(lsieved) == length(l)
      l = fields(9, 2, fmpz(10)^20)
      G = GAP.Globals.SmallGroup(27, 3)
      L = GAP.Globals.DerivedSeries(G)
      lsieved = Hecke.check_obstruction(l, L, 2, [3])
      @test length(lsieved) == 24
    end
  end

  println("Obstructions: prime power case")
  @time begin
    @testset "Obstructions: prime_power_case" begin
      l = fields(4, 2, fmpz(10)^5);
      G = GAP.Globals.SmallGroup(16, 9)
      L = GAP.Globals.DerivedSeries(G)
      lsieved = Hecke.check_obstruction(l, L, 2, [4])
      @test length(lsieved) == 78
      l = fields(4, 2, fmpz(10)^6);
      G = GAP.Globals.SmallGroup(32, 20)
      L = GAP.Globals.DerivedSeries(G)
      lsieved = Hecke.check_obstruction(l, L, 2, [8])
      @test length(lsieved) == 232
    end
  end

  println("Obstructions: reduction to cyclic case")
  @time begin
    @testset "Obstructions: reduction to cyclic case" begin
      l = fields(8, 2, fmpz(10)^10)
      G = GAP.Globals.SmallGroup(96, 13)
      L = GAP.Globals.DerivedSeries(G)
      lsieved = Hecke.check_obstruction(l, L, 2, [2, 6])
      @test length(l) == length(lsieved)
      @test length(fields(32, 30, fmpz(10)^39)) == 2
    end
  end

  println("Some examples")
  @time begin
    @testset "Some examples" begin
      lQ = fields(1, 1, fmpz(1))
      @test length(lQ) == 1
      @test degree(number_field(lQ[1])) == 1
      l1 = fields(6, 1, fmpz(10)^8)
      @test length(l1) == 107
      for x in l1
        @test GAP.gap_to_julia(Vector{Int}, Hecke.IdGroup(closure(x.generators_of_automorphisms))) == [6, 1]
      end
      l2 = fields(24, 8, fmpz(10)^30)
      @test length(l2) == 15
      for x in l2
        @test GAP.gap_to_julia(Vector{Int}, Hecke.IdGroup(closure(x.generators_of_automorphisms))) == [24, 8]
      end
      l3 = fields(30, 3, fmpz(10)^40)
      @test length(l3) == 5
      for x in l3
        @test GAP.gap_to_julia(Vector{Int}, Hecke.IdGroup(closure(x.generators_of_automorphisms))) == [30, 3]
      end
      Qx, x = PolynomialRing(FlintQQ, "x", cached = false)
      f = x^36 - x^33 + x^27 - x^24 + x^18 - x^12 + x^9 - x^3 + 1
      K, a = number_field(f, cached = false, check = false)
      d = fmpz(1252291517600545939502745293690906945712691746311040212121628914687318440182651069503694911322360563684969)
      flds = [Hecke.field_context(K)]
      @test isone(length(fields(72, 37, flds, d)))

      f = x^24 + 12*x^23 + 12*x^22 - 372*x^21 - 1224*x^20 + 4192*x^19 + 22370*x^18 - 17518*x^17 - 198383*x^16 - 30608*x^15 + 1017900*x^14 + 582696*x^13 - 3257224*x^12 - 2385744*x^11 + 6751956*x^10 + 4794056*x^9 - 9162405*x^8 - 4964804*x^7 + 7887120*x^6 + 2258008*x^5 - 3853388*x^4 - 85944*x^3 + 797010*x^2 - 160358*x + 1483
      K, a = number_field(f, check = false, cached = false)
      c = Hecke.field_context(K)
      @test length(c.subfields) == 3
      @test degree(domain(c.subfields[1])) == 1
      @test degree(domain(c.subfields[2])) == 2
      @test degree(domain(c.subfields[3])) == 6
      @test domain(c.subfields[2]) == codomain(c.subfields[1])
      @test domain(c.subfields[3]) == codomain(c.subfields[2])

      f = x^16 + 2*x^15 - 9*x^14 + 157*x^12 + 56*x^11 + 505*x^10 - 2324*x^9 + 1608*x^8 - 4252*x^7 + 19587*x^6 - 26800*x^5 - 52305*x^4 - 97006*x^3 + 164822*x^2 + 444940*x + 245969
      K, a = number_field(f, check = false, cached = false)
      l = [Hecke.field_context(K)];
      @test length(fields(64, 6, l, fmpz(80)^64)) == 1

    end
  end

end


text:
using SymbolicRegression, Test
include("test_params.jl")

_inv(x::Float32)::Float32 = 1.0f0 / x
X = rand(Float32, 5, 100) .+ 1
y = 1.2f0 .+ 2 ./ X[3, :]

options = SymbolicRegression.Options(;
    default_params..., binary_operators=(+, *), unary_operators=(_inv,), npopulations=8
)
hallOfFame = EquationSearch(X, y; niterations=8, options=options, numprocs=4)

dominating = calculate_pareto_frontier(X, y, hallOfFame, options)

best = dominating[end]
# Test the score
@test best.loss < maximum_residual / 10


text:
using BinaryProvider

# Parse some basic command-line arguments
const verbose = "--verbose" in ARGS
const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))

# These are the two binary objects we care about
products = Product[
    LibraryProduct(prefix, "libmdal", :libmdal),
    ExecutableProduct(prefix, "mdalinfo", :mdalinfo),
]

# Download binaries from hosted location
bin_prefix = "https://github.com/staticfloat/small_bin/raw/51b13b44feb2a262e2e04690bfa54d03167533f2/libfoo"
download_info = Dict(
    Linux(:aarch64, :glibc) => ("$bin_prefix/libfoo.aarch64-linux-gnu.tar.gz", "36886ac25cf5678c01fe20630b413f9354b7a3721c6a2c2043162f7ebd147ff5"),
    Linux(:armv7l, :glibc)  => ("$bin_prefix/libfoo.arm-linux-gnueabihf.tar.gz", "147ebaeb1a722da43ee08705689aed71ac87c3c2c907af047c6721c0025ba383"),
    Linux(:powerpc64le, :glibc) => ("$bin_prefix/libfoo.powerpc64le-linux-gnu.tar.gz", "5c35295ac161272ada9a77d1f6b770e30ea864e521e31853258cbc36ad4c4468"),
    Linux(:i686, :glibc)    => ("$bin_prefix/libfoo.i686-linux-gnu.tar.gz", "97655b6a218d61284723b6923d7c96e6a256fa68b9419d723c588aa24404b102"),
    Linux(:x86_64, :glibc)  => ("$bin_prefix/libfoo.x86_64-linux-gnu.tar.gz", "5208c63a9d07e592c78f541fc13caa8cd191b11e7e77b31d407237c2b13ec391"),

    Linux(:aarch64, :musl)  => ("$bin_prefix/libfoo.aarch64-linux-musl.tar.gz", "81751477c1e3ee6c93e1c28ee7db2b99d1eed0d6ce86dc30d64c2e5dd4dfe88d"),
    Linux(:armv7l, :musl)   => ("$bin_prefix/libfoo.arm-linux-musleabihf.tar.gz", "bb65aad58f2e6fc39dc9688da1bca5e8103a3a3fa67dc589debbd2e98176f0e1"),
    Linux(:i686, :musl)     => ("$bin_prefix/libfoo.i686-linux-musl.tar.gz", "5f02fd1fe19f3a565fb128d3673b35c7b3214a101cef9dcbb202c0092438a87b"),
    Linux(:x86_64, :musl)   => ("$bin_prefix/libfoo.x86_64-linux-musl.tar.gz", "ea630600a12d2c1846bc93bcc8d9638a4991f63329205c534d93e0a3de5f641d"),

    FreeBSD(:x86_64)        => ("$bin_prefix/libfoo.x86_64-unknown-freebsd11.1.tar.gz", "5f6edd6247b3685fa5c42c98a53d2a3e1eef6242c2bb3cdbb5fe23f538703fe4"),
    MacOS(:x86_64)          => ("$bin_prefix/libfoo.x86_64-apple-darwin14.tar.gz", "fcc268772d6f21d65b45fcf3854a3142679b78e53c7673dac26c95d6ccc89a24"),

    Windows(:i686)          => ("$bin_prefix/libfoo.i686-w64-mingw32.tar.gz", "79181cf62ca8e0b2e0851fa0ace52f4ab335d0cad26fb7f9cd4ff356a9a96e70"),
    Windows(:x86_64)        => ("$bin_prefix/libfoo.x86_64-w64-mingw32.tar.gz", "7f8939e9529835b83810d3ae7e2556f6e002d571f619894e54ece42ea5262b7f"),
)
# First, check to see if we're all satisfied
if any(!satisfied(p; verbose=verbose) for p in products)
    try
        # Download and install binaries
        url, tarball_hash = choose_download(download_info)
        install(url, tarball_hash; prefix=prefix, force=true, verbose=true)
    catch e
        if typeof(e) <: ArgumentError
            error("Your platform $(Sys.MACHINE) is not supported by this package!")
        else
            rethrow(e)
        end
    end

    # Finally, write out a deps.jl file
    write_deps_file(joinpath(@__DIR__, "deps.jl"), products)
end
write_deps_file(joinpath(@__DIR__, "deps.jl"), products)


text:
using QuadraticRelaxations

using Mosek, MosekTools

using Dates

const ROOT_DIR = joinpath(@__DIR__, "..")
include(joinpath(ROOT_DIR, "local_config.jl"))
include(joinpath(ROOT_DIR, "solver_factories.jl"))
include(joinpath(ROOT_DIR, "box_qp/util.jl"))
include(joinpath(ROOT_DIR, "box_qp/justitia_interface.jl"))
include("justitia_addendum.jl")

const TIME_LIMIT = 10 * 60.0
const NUM_THREADS = 1

const NUM_LAYERS = 3

const SHIFT = QuadraticRelaxations.SemidefiniteShift(Mosek.Optimizer)
const SOC_LOWER_BOUND = false

const FAMILIES = ("basic", "extended", "extended2")

const INSTANCE_DIR = joinpath(ROOT_DIR, "box_qp", "data")
const OUTPUT_ROOT_DIR = joinpath(@__DIR__, "output/")

const CONFIG = FactoryConfig(TIME_LIMIT, NUM_THREADS, CPLEX_PATH)

INSTANCES = Dict{String,NearestGoodBoxQPInstance}()
for family in FAMILIES, instance_name in readdir(joinpath(INSTANCE_DIR, family))
    INSTANCES["$family-$instance_name"] =
        _create_nearest_good_instance(family, instance_name, primal_factory(CONFIG))
end

APPROACHES = Dict{String,BoxQPApproach}()
for (s_name, solver) in [
    "brn" => baron_direct_factory,
    "grb" => gurobi_direct_factory,
    "grb3" => gurobi_direct_mipfocus_3_factory,
]
    for (t_name, transformation) in [
        "direct" => identity,
        "direct-shift" => reformulate_with(
            QuadraticRelaxations.ExactFormulation(),
            SHIFT,
            SOC_LOWER_BOUND,
        ),
    ]
        APPROACHES[string(s_name, "-", t_name)] = BoxQPApproach(solver, transformation)
    end
end
for (s_name, solver) in
    ["grb" => gurobi_mip_factory, "grb3" => gurobi_mip_mipfocus_3_factory]
    for (r_name, relaxation) in [
        "nn" => QuadraticRelaxations.NeuralNetRelaxation(NUM_LAYERS),
        "cda" => QuadraticRelaxations.CDARelaxation(NUM_LAYERS),
        "standard-nmdt" => QuadraticRelaxations.StandardNMDT(NUM_LAYERS),
        "tightened-nmdt" => QuadraticRelaxations.TightenedNMDT(NUM_LAYERS),
    ]
        APPROACHES[string(s_name, "-", r_name)] =
            BoxQPApproach(solver, reformulate_with(relaxation, SHIFT, SOC_LOWER_BOUND))
    end
end


text:
using Gen
import Random

#########
# model #
#########

struct Params
    prob_outlier::Float64
    inlier_std::Float64
    outlier_std::Float64
    slope::Float64
    intercept::Float64
end

@gen function datum(x::Float64, params::Params)
    is_outlier = @addr(bernoulli(params.prob_outlier), :z)
    std = is_outlier ? params.inlier_std : params.outlier_std
    y = @addr(normal(x * params.slope + params.intercept, std), :y)
    return y
end

data = plate(datum)

@gen function model(xs::Vector{Float64})
    inlier_std = @addr(Gen.gamma(1, 1), :inlier_std)
    outlier_std = @addr(Gen.gamma(1, 1), :outlier_std)
    slope = @addr(normal(0, 2), :slope)
    intercept = @addr(normal(0, 2), :intercept)
    params = Params(0.5, inlier_std, outlier_std, slope, intercept)
    if all([@change(addr) == NoChange() for addr in [:slope, :intercept, :inlier_std, :outlier_std]])
        change = NoChange()
    else
        change = nothing
    end
    ys = @addr(data(xs, fill(params, length(xs))), :data, change)
    return ys
end

#######################
# inference operators #
#######################

@gen function slope_proposal(prev)
    slope = get_assignment(prev)[:slope]
    @addr(normal(slope, 0.5), :slope)
end

@gen function intercept_proposal(prev)
    intercept = get_assignment(prev)[:intercept]
    @addr(normal(intercept, 0.5), :intercept)
end

@gen function inlier_std_proposal(prev)
    inlier_std = get_assignment(prev)[:inlier_std]
    @addr(normal(inlier_std, 0.5), :inlier_std)
end

@gen function outlier_std_proposal(prev)
    outlier_std = get_assignment(prev)[:outlier_std]
    @addr(normal(outlier_std, 0.5), :outlier_std)
end

@gen function is_outlier_proposal(prev, i::Int)
    prev = get_assignment(prev)[:data => i => :z]
    @addr(bernoulli(prev ? 0.0 : 1.0), :data => i => :z)
end

@gen function observer(ys::Vector{Float64})
    for (i, y) in enumerate(ys)
        @addr(dirac(y), :data => i => :y)
    end
end

Gen.load_generated_functions()

#####################
# generate data set #
#####################

Random.seed!(1)

prob_outlier = 0.5
true_inlier_noise = 0.5
true_outlier_noise = 5.0
true_slope = -1
true_intercept = 2
xs = collect(range(-5, stop=5, length=200))
ys = Float64[]
for (i, x) in enumerate(xs)
    if rand() < prob_outlier
        y = true_slope * x + true_intercept + randn() * true_inlier_noise
    else
        y = true_slope * x + true_intercept + randn() * true_outlier_noise
    end
    push!(ys, y)
end

##################
# run experiment #
##################

function do_inference(n)
    observations = get_assignment(simulate(observer, (ys,)))

    # initial trace
    (trace, _) = generate(model, (xs,), observations)

    for i=1:n

        # steps on the parameters
        for j=1:5
            trace = mh(model, slope_proposal, (), trace)
            trace = mh(model, intercept_proposal, (), trace)
            trace = mh(model, inlier_std_proposal, (), trace)
            trace = mh(model, outlier_std_proposal, (), trace)
        end

        # step on the outliers
        for j=1:length(xs)
            trace = mh(model, is_outlier_proposal, (j,), trace)
        end

        score = get_call_record(trace).score

        # print
        assignment = get_assignment(trace)
        slope = assignment[:slope]
        intercept = assignment[:intercept]
        inlier_std = assignment[:inlier_std]
        outlier_std = assignment[:outlier_std]
        println("score: $score, slope: $slope, intercept: $intercept, inlier_std: $inlier_std, outlier_std: $outlier_std")
    end
end

@time do_inference(100)
@time do_inference(100)

max_stars_repo_name: BoltzmannBrain/Gen

max_issues_repo_name: BoltzmannBrain/Gen

max_forks_repo_name: BoltzmannBrain/Gen


text:
# # Ice Flow over a Bumpy Bed

# This is the first two dimensional example: flow down a plane with bumps on different length scales.
# This example is similar in spirit to Experiment A in the Ice Sheet Model Intercomparison Exercise - Higher Order Model (ISMIP-HOM): doi: 10.5194/tc-2-95-2008
# This example demonstrates: How to use WAVI.jl in plan view (two horizontal spatial dimensions) with arbitrary bed shapes, and how to interact with Grid objects.

# ## Install dependencies
#
# First let's make sure we have all required packages installed.

# using Pkg
# Pkg.add("https://github.com/RJArthern/WAVI.jl"), Pkg.add(Plots)
using WAVI, Plots

# ## Basal Topography
#
# Following ISMIP-HOM, let's take a bed with a series of sinusoidal oscillations with an amplitude of 500m:
# $$
# z_b(x,y) = -x \tan \alpha + 500 \sin (\omega x) \sin(\omega y)
# $$
# Here $L$ is the lengthscale of the domain, $\alpha$ is the net slope of the plane, and $\omega = 2\pi / L$ is the frequency of the bumps.
# Let's define this bed as a function and take a look at it:

z_b(x,y; α, ω) = -x * tand(α) + 500sin(ω*x)*sin(ω*y);

# It's useful to define our model grid here. We can do this for an arbitary domain length L.
grid(; L, nx = 80,ny = 80) = Grid(nx = nx, ny = ny, dx = L/nx, dy = L/ny, y0 = 0.0, x0 = 0.0);

# The final two arguments simply set the origin of the co-ordinate system. This grid has 80 grid points in each direction by default

# Let's choose a domain of 80km with 80 grid points
L = 80000.;
grid80 = grid(L = L);

# This grid object contains information about the location of grid points. We use this to construct and array defining the bed:
z_b80 = z_b.(grid80.xxh,grid80.yyh; α = 0.5, ω = 2π/L );
plt = Plots.heatmap(grid80.xxh[:,1]/1e3, grid80.yyh[1,:]/1e3, z_b80,
                        xlabel = "x (km)",
                        ylabel = "y (km)",
                        colorbar_title = "bed depth (m)")
plot!(size = (800,400))
#display(plt)

# ## Model Instantiation and Initial Conditions
# In the ISMIP-HOM comparison, the main test is velocity along the line $y = L/4$, for various different values of $L$.
# Before we do that, lets look at the velocity for the example we started above with $L = 80$km.

# To begin, we create an `InitialConditions` object to prescribe the ice thickness of 1000m everywhere
initial_conditions = InitialConditions(initial_thickness = 1000. .* ones(grid80.nx, grid80.ny));

# Now we can build our model
model80 = Model(grid = grid80,
            bed_elevation = z_b80,
            initial_conditions = initial_conditions);

# ## Determining the velocity
# To bring the velocity in line with the ice thickness, we have to use the `update_state!` function
update_state!(model80);

# Let's have a look at the velocity component in the x-direction
Plots.heatmap(model80.grid.xxh[:,1]/1e3, model80.grid.yyh[1,:]/1e3, model80.fields.gh.u',
                        xlabel = "x (km)",
                        ylabel = "y (km)",
                        colorbar_title = "ice velocity in x-direction (m/yr)")
plot!(size = (1000,550))

# ## Different Lengthscales
# Now let's look at how the velocity along a flowline changes with L.
L_values =  [160, 80, 40, 20, 10, 5]*1.0e3;

# We loop over these values and store the info:
U_flowline = zeros(80, length(L_values));
grid_flowline = zeros(80, length(L_values));
for (count,L) in enumerate(L_values) ;
    gridL = grid(L = L);
    z_bL = z_b.(gridL.xxh,gridL.yyh; α = 0.5, ω = 2π/L );
    initial_conditions = InitialConditions(initial_thickness = 1000. .* ones(80, 80));
    model = Model(grid = gridL,
                bed_elevation = z_bL,
                initial_conditions = initial_conditions);
    update_state!(model);
    grid_flowline[:, count] .= model.grid.xxh[:, round(Int, gridL.nx/4)];
    U_flowline[:,count] .= model.fields.gh.u[:, round(Int, gridL.nx/4)];
end

#And make the plot in a single command
plot(grid_flowline/1e3,
    U_flowline,
    layout = (2,3),
    framestyle = :box,
    xlabel = "x (km)",
    ylabel = "horizontal velocity (m/yr)",
    label = :none,
    title =  ["160km" "80km" "40km" "20km" "10km" "5km"])
plot!(size = (1000,550))


text:
module quantum

using ..interaction, ..system
using QuantumOptics, LinearAlgebra

try
    eval(Expr(:toplevel,:(import Optim)))
    global optimize = Optim.optimize
catch e
    if typeof(e) == ArgumentError
        println("Optim package not available. (Needed for calculation of squeezing parameter)")
    else
        rethrow(e)
    end
end

export Hamiltonian, JumpOperators

# Define Spin 1/2 operators
const spinbasis = SpinBasis(1//2)
const sigmax_ = sigmax(spinbasis)
const sigmay_ = sigmay(spinbasis)
const sigmaz_ = sigmaz(spinbasis)
const sigmap_ = sigmap(spinbasis)
const sigmam_ = sigmam(spinbasis)
const I_spin = identityoperator(spinbasis)

"""
    quantum.basis(x)

Get basis of the given System.
"""
basis(x::Spin) = spinbasis
basis(x::SpinCollection) = CompositeBasis([basis(s) for s=x.spins]...)
basis(N::Int) = CompositeBasis([spinbasis for s=1:N]...)
basis(x::CavityMode) = FockBasis(x.cutoff)
basis(x::CavitySpinCollection) = compose(basis(x.cavity), basis(x.spincollection))

"""
    quantum.blochstate(phi, theta[, N=1])

Product state of `N` single spin Bloch states.

All spins have the same azimuthal angle `phi` and polar angle `theta`.
"""
function blochstate(phi::Vector{T1}, theta::Vector{T2}) where {T1<:Real, T2<:Real}
    N = length(phi)
    @assert length(theta)==N
    state_g = basisstate(spinbasis, 1)
    state_e = basisstate(spinbasis, 2)

    states = [cos(theta[k]/2)*state_g + exp(1im*phi[k])*sin(theta[k]/2)*state_e for k=1:N]
    return reduce(tensor, states)
end

function blochstate(phi::Real, theta::Real, N::Int=1)
    state_g = basisstate(spinbasis, 1)
    state_e = basisstate(spinbasis, 2)
    state = cos(theta/2)*state_g + exp(1im*phi)*sin(theta/2)*state_e
    if N>1
        return reduce(tensor, [state for i=1:N])
    else
        return state
    end
end

"""
    quantum.dim(state)

Number of spins described by this state.
"""
function dim(ρ::AbstractOperator)
    return length(ρ.basis_l.bases)
end

"""
    quantum.Hamiltonian(S)

Hamiltonian of the given System.
"""
function Hamiltonian(S::system.SpinCollection)
    spins = S.spins
    N = length(spins)
    b = basis(S)
    H = SparseOperator(b)
    for i=1:N
        if S.spins[i].delta != 0.
            H += 0.5*S.spins[i].delta * embed(b, i, sigmaz_)
        end
    end
    for i=1:N, j=1:N
        if i==j
            continue
        end
        sigmap_i = embed(b, i, sigmap_)
        sigmam_j = embed(b, j, sigmam_)
        H += interaction.Omega(spins[i].position, spins[j].position, S.polarizations[i], S.polarizations[j], S.gammas[i], S.gammas[j])*sigmap_i*sigmam_j
    end
    return H
end

function Hamiltonian(S::system.CavityMode)
    b = basis(S)
    H = SparseOperator(S.delta*number(b) + S.eta*(create(b) + destroy(b)))
    return H
end

function Hamiltonian(S::system.CavitySpinCollection)
    b = basis(S)
    bs = basis(S.spincollection)
    bc = basis(S.cavity)
    Hs = Hamiltonian(S.spincollection)
    Hc = Hamiltonian(S.cavity)
    Ic = identityoperator(basis(S.cavity))
    H = embed(b, 1, Hc) + tensor(Ic, Hs)
    a = SparseOperator(destroy(bc))
    at = SparseOperator(create(bc))
    for i=1:length(S.spincollection.spins)
        if S.g[i] != 0.
            H += S.g[i]*(tensor(a, embed(bs, i, sigmap_)) + tensor(at, embed(bs, i, sigmam_)))
        end
    end
    return H
end

"""
    quantum.JumpOperators(S)

Jump operators of the given system.
"""
function JumpOperators(S::system.SpinCollection)
    J = SparseOpType[embed(basis(S), i, sigmam_) for i=1:length(S.spins)]
    Γ = interaction.GammaMatrix(S)
    return Γ, J
end

JumpOperators(S::system.CavityMode) = (Float64[2*S.kappa], SparseOpType[SparseOperator(destroy(basis(S)))])

function JumpOperators(S::system.CavitySpinCollection)
    Γs, Js = JumpOperators(S.spincollection)
    Γc, Jc = JumpOperators(S.cavity)
    Ns = length(Js)
    Nc = length(Jc)
    N = Ns + Nc
    T = promote_type(eltype(Γs),eltype(Γc))
    Γ = zeros(T, N, N)
    Γ[1:Nc, 1:Nc] = Γc
    Γ[Nc+1:end, Nc+1:end] = Γs
    Ic = identityoperator(basis(S.cavity))
    J = SparseOpType[embed(basis(S), 1, Jc[1])]
    for j in Js
        push!(J, tensor(Ic, j))
    end
    return Γ, J
end

"""
    quantum.JumpOperators_diagonal(S)

Jump operators of the given system. (diagonalized)

Diagonalized means that the Gamma matrix is diagonalized and
the jump operators are changed accordingly.
"""
function JumpOperators_diagonal(S::system.SpinCollection)
    spins = S.spins
    N = length(spins)
    b = basis(S)
    Γ = [interaction.Gamma(spins[i].position, spins[j].position, S.polarizations[i], S.polarizations[j], S.gammas[i], S.gammas[j]) for i=1:N, j=1:N]
    λ, M = eig(Γ)
    J = Any[]
    for i=1:N
        op = Operator(b)
        for j=1:N
            op += M[j,i]*embed(b, j, sigmam_)
        end
        push!(J, sqrt(λ[i])*op)
    end
    return J
end

"""
    quantum.timeevolution_diagonal(T, S, state0[; fout])

Master equation time evolution. (diagonalized)

Diagonalized means that the Gamma matrix is diagonalized and
the jump operators are changed accordingly.

# Arguments
* `T`: Points of time for which output will be generated.
* `S`: System
* `ρ₀`: Initial density operator.
* `fout` (optional): Function with signature fout(t, state) that is called
    whenever output should be generated.
"""
function timeevolution_diagonal(T, S::system.System, ρ₀::Union{StateVector, DenseOpType}; fout=nothing, kwargs...)
    H = Hamiltonian(S)
    J = JumpOperators_diagonal(S)
    Hnh = H - 0.5im*sum([dagger(J[i])*J[i] for i=1:length(J)])
    return QuantumOptics.timeevolution.master_nh(T, ρ₀, Hnh, J; fout=fout, kwargs...)
end

"""
    quantum.timeevolution(T, S, state0[; fout])

Master equation time evolution.

Diagonalized means that the Gamma matrix is diagonalized and
the jump operators are changed accordingly.

# Arguments
* `T`: Points of time for which output will be generated.
* `S`: System
* `ρ₀`: Initial density operator.
* `fout` (optional): Function with signature fout(t, state) that is called
    whenever output should be generated.
"""
function timeevolution(T, S::system.System, ρ₀::Union{StateVector, DenseOpType}; fout=nothing, kwargs...)
    b = basis(S)
    H = Hamiltonian(S)
    Γ, J = JumpOperators(S)
    return QuantumOptics.timeevolution.master_h(T, ρ₀, H, J; fout=fout, rates=Γ, kwargs...)
end

"""
    meanfield.rotate(axis, angles, state)

Rotations on the Bloch sphere for the given density operator.

# Arguments
* `axis`: Rotation axis.
* `angles`: Rotation angle(s).
* `ρ`: Density operator that should be rotated.
"""
function rotate(axis::Vector{T1}, angles::Vector{T2}, ρ::DenseOpType) where {T1<:Real, T2<:Real}
    N = dim(ρ)
    @assert length(axis)==3
    @assert length(angles)==N
    basis = ρ.basis_l
    n = axis/norm(axis)
    for i=1:N
        nσ = n[1]*sigmax_ + n[2]*sigmay_ + n[3]*sigmaz_
        α = angles[i]
        R = I_spin*cos(α/2) - 1im*nσ*sin(α/2)
        R_ = embed(basis, i, R)
        ρ = R_*ρ*dagger(R_)
    end
    return ρ
end

rotate(axis::Vector{T}, angle::Real, ρ::DenseOpType) where {T<:Real} = rotate(axis, ones(float(T), dim(ρ))*angle, ρ)

"""
    quantum.squeeze_sx(χT, ρ₀)

Spin squeezing along sx.

# Arguments
* `χT`: Squeezing strength.
* `ρ₀`: Operator that should be squeezed.
"""
function squeeze_sx(χT::Real, ρ₀::DenseOpType)
    N = dim(ρ₀)
    basis = ρ₀.basis_l
    totaloperator(op::SparseOpType) = sum([embed(basis, i, op) for i=1:N])/N
    sigmax_total = totaloperator(sigmax_)
    H = χT*sigmax_total^2
    T = [0.,1.]
    t, states = QuantumOptics.timeevolution.master(T, ρ₀, H, [])
    return states[end]
end

"""
    quantum.squeeze(axis, χT, ρ₀)

Spin squeezing along an arbitrary axis.

# Arguments
* `axis`: Squeezing axis.
* `χT`: Squeezing strength.
* `ρ₀`: Operator that should be squeezed.
"""
function squeeze(axis::Vector{T}, χT::Real, ρ₀::DenseOpType) where T<:Real
    @assert length(axis)==3
    axis = axis/norm(axis)
    N = dim(ρ₀)
    basis = ρ₀.basis_l
    totaloperator(op::SparseOpType) = sum([embed(basis, i, op) for i=1:N])/N
    σ = map(totaloperator, [sigmax_, sigmay_, sigmaz_])
    σn = sum([axis[i]*σ[i] for i=1:3])
    H = χT*σn^2
    tout, states = QuantumOptics.timeevolution.master([0,1], ρ₀, H, [])
    return states[end]
end

"""
    quantum.orthogonal_vectors(n)

Create 3 orthonormal vectors where one is in the given direction `n`.
"""
function orthogonal_vectors(n::Vector{<:Real})
    @assert length(n)==3
    n = n/norm(n)
    v = (n[1]<n[2] ? [1.,0.,0.] : [0.,1.,0.])
    e1 = v - dot(n,v)*n
    e1 = e1/norm(e1)
    e2 = cross(n, e1)
    e2 = e2/norm(e2)
    return e1, e2
end

"""
    quantum.variance(op, state)

Variance of the operator for the given state.
"""
variance(op::AbstractOperator, state) = (expect(op^2, state) - expect(op, state)^2)

"""
    quantum.squeezingparameter(ρ)

Calculate squeezing parameter for the given state.
"""
function squeezingparameter(ρ::DenseOpType)
    N = dim(ρ)
    basis = ρ.basis_l
    totaloperator(op::SparseOpType) = sum([embed(basis, i, op) for i=1:N])/N
    S = map(totaloperator, [sigmax_, sigmay_, sigmaz_])
    n = real([expect(s, ρ) for s=S])
    e1, e2 = orthogonal_vectors(n)
    function f(phi)
        nphi = cos(phi)*e1 + sin(phi)*e2
        Sphi = sum([nphi[i]*S[i] for i=1:3])
        return real(variance(Sphi, ρ))
    end
    varSmin = Optim.optimize(f, 0., 2*pi).f_minimum
    return sqrt(N*varSmin)/norm(n)
end

end # module


text:
# MvNormalCDF
# Copyright © 2019-2021 Vladimir Arnautov aka PharmCat <mail@pharmcat.net>, Andrew Gough

module MvNormalCDF
    using Distributions, Primes, Random, LinearAlgebra, FillArrays, StatsBase

    export MvNormal, mvnormcdf

    const sqrt2π = sqrt(2π)

    include("functions.jl")
end # module

max_stars_repo_name: PharmCat/MvNormalCDF

max_issues_repo_name: PharmCat/MvNormalCDF

max_forks_repo_name: PharmCat/MvNormalCDF


text:
using Knet: @diff, Param, value, grad

sigm(x) = 1.0 / (1.0 + exp(-x))

mutable struct Layer
    wai::Param
    was::Param
    ba::Param
    wri::Param
    wrs::Param
    br::Param
    wsi::Param
    bs::Param
    state
end

Layer(in_size, out_size) =
begin
    sq = sqrt(in_size)
    wai = Param(randn(in_size, out_size)/sq)
    was = Param(randn(out_size, out_size)/sq)
    ba  = Param(zeros(1, out_size))
    wri = Param(randn(in_size, out_size)/sq)
    wrs = Param(randn(out_size, out_size)/sq)
    br  = Param(zeros(1, out_size))
    wsi = Param(randn(in_size, out_size)/sq)
    bs  = Param(zeros(1, out_size))
    state = zeros(1, out_size)
    layer = Layer(wai, was, ba, wri, wrs, br, wsi, bs, state)
layer
end

(layer::Layer)(in) =
begin
    attention = sigm.(in * layer.wai + layer.state * layer.was + layer.ba)
    remember  = sigm.(in * layer.wri + layer.state * layer.wrs + layer.br)
    short_mem = tanh.(in * layer.wsi + attention .* layer.state + layer.bs)
    layer.state = remember .* short_mem + (1 .- remember) .* layer.state
end

make(in_size, layers, out_size) =
begin
    model = []
    hm_layers = length(layers) +1
    for i in 1:hm_layers
        if     i == 1         layer = Layer(in_size, layers[i])
        elseif i == hm_layers layer = Layer(layers[end], out_size)
        else                  layer = Layer(layers[i-1], layers[i])
        end
        push!(model, layer)
    end
model
end

prop(model, in) =
begin
    for layer in model
        in = layer(in)
    end
    for layer in model
        setfield!(layer, :state, zeros(1, size(getfield(layer, :bs))[end]))
    end
in
end

cross_entropy(out, label) =
begin
    - sum(label .* log.(out))
end


text:
using Test
using JDF
using RDatasets
using DataFrames: DataFrame
using CategoricalArrays
using Tables

@testset "JDF.jl categorical arrays" begin
    df = DataFrame([collect(1:100) for i = 1:3], :auto)
    df[!, :x1] = categorical(df[!, :x1])
    df[!, :x2] = categorical(string.(df[!, :x2]))
    JDF.save(df, "a3cate.jdf")
    df_loaded_back = JDF.load("a3cate.jdf", cols = [:x2, :x1])

    df2 = DataFrame(df_loaded_back; copycols = true)
    @test size(df2, 2) == 2
    @test size(df2, 1) == 100
    @time df2[!, :x1] isa CategoricalVector{Int}
    @time df2[!, :x2] isa CategoricalVector{String}

    rm("a3cate.jdf", force = true, recursive = true)
end

@testset "Guard against Github #27" begin
    iris = dataset("datasets", "iris")
    JDF.save(iris, "iris.jdf")
    JDF.load("iris.jdf")

    rm("iris.jdf", force = true, recursive = true)
end

@testset "CategoricalArray{Union{Missing, String}}" begin
    # Guard against github 73
    df2 = DataFrame(sex = categorical(["Male", missing, "Female"]))
    JDF.save("df2.jdf", df2)

    b = JDF.load("df2.jdf") |> DataFrame

    @test any(ismissing, b.sex)

    rm("df2.jdf", force=true, recursive=true)
end

max_stars_repo_name: xiaodaigh/JDF

max_issues_repo_name: xiaodaigh/JDF

max_forks_repo_name: xiaodaigh/JDF


text: # European option pricing using simulation.


text:
# f_3p_1im_dep(x) = 5.0 + (x[1]-3.0)^2 + (x[1]-x[2]-1.0)^2 + 0*x[3]^2 # [3., 2., any]

### base example

@testset "default options" begin
    res0 = [get_endpoint(
        [3., 2., 2.1],
        i,
        (x::Vector{Float64}) -> f_3p_1im_dep(x),
        :CICO_ONE_PASS;
        loss_crit = 9.,
        silent = true
    ) for i in 1:3]

    @test isapprox(res0[1].value, 5.0, atol=1e-2)
    @test length(res0[1].profilePoints) > 0
    @test res0[1].status == :BORDER_FOUND_BY_SCAN_TOL
    @test res0[1].direction == :right
    @test isapprox(res0[2].value, 2.0+2.0*sqrt(2.), atol=1e-2)
    @test length(res0[2].profilePoints) > 0
    @test res0[2].status == :BORDER_FOUND_BY_SCAN_TOL
    @test res0[2].direction == :right
    @test length(res0[3].profilePoints) == 0
    @test res0[3].status == :SCAN_BOUND_REACHED
    @test res0[3].direction == :right
end

@testset ":left" begin
    res0 = [get_endpoint(
        [3., 2., 2.1],
        i,
        (x::Vector{Float64}) -> f_3p_1im_dep(x),
        :CICO_ONE_PASS,
        :left;
        loss_crit = 9.,
        silent = true
    ) for i in 1:3]

    @test isapprox(res0[1].value, 1.0, atol=1e-2)
    @test length(res0[1].profilePoints) > 0
    @test res0[1].status == :BORDER_FOUND_BY_SCAN_TOL
    @test res0[1].direction == :left
    @test isapprox(res0[2].value, 2.0-2.0*sqrt(2.), atol=1e-2)
    @test length(res0[2].profilePoints) > 0
    @test res0[2].status == :BORDER_FOUND_BY_SCAN_TOL
    @test res0[2].direction == :left
    @test length(res0[3].profilePoints) == 0
    @test res0[3].status == :SCAN_BOUND_REACHED
    @test res0[3].direction == :left
end

@testset ":log" begin
    res0 = [get_endpoint(
        [3., 2., 2.1],
        i,
        (x::Vector{Float64}) -> f_3p_1im_dep(x),
        :CICO_ONE_PASS,
        :right;
        loss_crit = 9.,
        scale = [:log, :direct, :log],
        silent = true
    ) for i in 1:3]

    @test isapprox(log10(res0[1].value), log10(5.), atol=1e-2)
    @test length(res0[1].profilePoints) > 0
    @test res0[1].status == :BORDER_FOUND_BY_SCAN_TOL
    @test res0[1].direction == :right
    @test isapprox(res0[2].value, 2.0+2.0*sqrt(2.), atol=1e-2)
    @test length(res0[2].profilePoints) > 0
    @test res0[2].status == :BORDER_FOUND_BY_SCAN_TOL
    @test res0[2].direction == :right
    @test length(res0[3].profilePoints) == 0
    @test res0[3].status == :SCAN_BOUND_REACHED
    @test res0[3].direction == :right
end

### derivation free

@testset ":LN_NELDERMEAD, no scale" begin
    res0 = get_endpoint(
        [3., 2., 2.1],
        1,            # parameter number
        f_3p_1im_dep, # loss_func
        :CICO_ONE_PASS,
        loss_crit = 9.,
        local_alg = :LN_NELDERMEAD,
        silent = true
    )
    @test res0.status == :BORDER_FOUND_BY_SCAN_TOL
    @test res0.counter < 1000
    @test isapprox(res0.value, 5.; atol=1e-2)
end

@testset ":LN_NELDERMEAD, no scale, !unpredictable" begin
    res0 = get_endpoint(
        [3., 2., 2.1],
        3,            # parameter number
        f_3p_1im_dep, # loss_func
        :CICO_ONE_PASS,
        loss_crit = 9.,
        local_alg = :LN_NELDERMEAD,
        silent = true
    )
    @test res0.status == :SCAN_BOUND_REACHED
    @test res0.counter < 1000
end

@testset ":LN_NELDERMEAD, log scale" begin
    res0 = get_endpoint(
        [3., 2., 2.1],
        1,            # parameter number
        f_3p_1im_dep, # loss_func
        :CICO_ONE_PASS,
        loss_crit = 9.,
        scale = [:log, :log, :log],
        local_alg = :LN_NELDERMEAD,
        silent = true
    )
    @test res0.status == :BORDER_FOUND_BY_SCAN_TOL
    @test res0.counter < 1000
    @test isapprox(res0.value, 5.; atol=1e-2)
end

@testset ":LN_NELDERMEAD, log scale, !unpredictable" begin
    res0 = get_endpoint(
        [3., 2., 2.1],
        3,            # parameter number
        f_3p_1im_dep, # loss_func
        :CICO_ONE_PASS,
        loss_crit = 9.,
        scale = [:log, :log, :log],
        local_alg = :LN_NELDERMEAD,
        silent = true
    )
    @test res0.status == :SCAN_BOUND_REACHED
    @test res0.counter < 1000
end

### gradient-based

@testset ":LD_SLSQP, no scale, no loss_grad => error" begin
    @test_throws ArgumentError get_endpoint(
        [3., 2., 2.1],
        1,            # parameter number
        f_3p_1im_dep, # loss_func
        :CICO_ONE_PASS,
        loss_crit = 9.,
        local_alg = :LD_SLSQP,
        silent = true
    )
end

@testset ":LD_SLSQP, log scale, no loss_grad => error" begin
    @test_throws ArgumentError get_endpoint(
        [3., 2., 2.1],
        1,            # parameter number
        f_3p_1im_dep, # loss_func
        :CICO_ONE_PASS,
        loss_crit = 9.,
        scale = [:log, :log, :log],
        local_alg = :LD_SLSQP,
        silent = true
    )
end

## autodiff

@testset ":LD_SLSQP, no scale, :AUTODIFF" begin
    res0 = get_endpoint(
        [3., 2., 2.1],
        1,            # parameter number
        f_3p_1im_dep, # loss_func
        :CICO_ONE_PASS,
        loss_crit = 9.,
        local_alg = :LD_SLSQP,
        loss_grad = :AUTODIFF,
        silent = true
    )
    @test res0.status == :BORDER_FOUND_BY_SCAN_TOL
    @test res0.counter < 1000
    @test isapprox(res0.value, 5.; atol=1e-2)
end

@testset ":LD_SLSQP, no scale, :AUTODIFF, !unpredictable" begin
    res0 = get_endpoint(
        [3., 2., 2.1],
        3,            # parameter number
        f_3p_1im_dep, # loss_func
        :CICO_ONE_PASS,
        loss_crit = 9.,
        local_alg = :LD_SLSQP,
        loss_grad = :AUTODIFF,
        silent = true
    )
    @test res0.status == :SCAN_BOUND_REACHED
    @test res0.counter < 1000
end

@testset ":LD_SLSQP, log scale, :AUTODIFF" begin # XXX: scan_tol should be updated
    res0 = get_endpoint(
        [3., 2., 2.1],
        1,            # parameter number
        f_3p_1im_dep, # loss_func
        :CICO_ONE_PASS,
        loss_crit = 9.,
        scale = [:log, :log, :log],
        local_alg = :LD_SLSQP,
        loss_grad = :AUTODIFF,
        scan_tol = 1e-4,
        silent = true
    )
    @test res0.status == :BORDER_FOUND_BY_SCAN_TOL
    @test res0.counter < 1000
    @test isapprox(res0.value, 5.; atol=1e-2)
end

@testset ":LD_SLSQP, log scale, :AUTODIFF, !unpredictable" begin
    res0 = get_endpoint(
        [3., 2., 2.1],
        3,            # parameter number
        f_3p_1im_dep, # loss_func
        :CICO_ONE_PASS,
        loss_crit = 9.,
        scale = [:log, :log, :log],
        local_alg = :LD_SLSQP,
        loss_grad = :AUTODIFF,
        silent = true
    )
    @test res0.status == :SCAN_BOUND_REACHED
    @test res0.counter < 1000
end

## FINITE

@testset ":LD_SLSQP, no scale, :FINITE" begin
    res0 = get_endpoint(
        [3., 2., 2.1],
        1,            # parameter number
        f_3p_1im_dep, # loss_func
        :CICO_ONE_PASS,
        loss_crit = 9.,
        local_alg = :LD_SLSQP,
        loss_grad = :FINITE,
        silent = true
    )
    @test res0.status == :BORDER_FOUND_BY_SCAN_TOL
    @test res0.counter < 1000
    @test isapprox(res0.value, 5.; atol=1e-2)
end

@testset ":LD_SLSQP, no scale, :FINITE, !unpredictable" begin
    res0 = get_endpoint(
        [3., 2., 2.1],
        3,            # parameter number
        f_3p_1im_dep, # loss_func
        :CICO_ONE_PASS,
        loss_crit = 9.,
        local_alg = :LD_SLSQP,
        loss_grad = :FINITE,
        silent = true
    )
    @test res0.status == :SCAN_BOUND_REACHED
    @test res0.counter < 1000
end

@testset ":LD_SLSQP, log scale, :FINITE" begin # XXX: scan_tol should be updated
    res0 = get_endpoint(
        [3., 2., 2.1],
        1,            # parameter number
        f_3p_1im_dep, # loss_func
        :CICO_ONE_PASS,
        loss_crit = 9.,
        scale = [:log, :log, :log],
        local_alg = :LD_SLSQP,
        loss_grad = :FINITE,
        scan_tol = 1e-4,
        silent = true
    )
    @test res0.status == :BORDER_FOUND_BY_SCAN_TOL
    @test res0.counter < 1000
    @test isapprox(log10(res0.value), log10(5.); atol=1e-2)
end

@testset ":LD_SLSQP, log scale, :FINITE, !unpredictable" begin
    res0 = get_endpoint(
        [3., 2., 2.1],
        3,            # parameter number
        f_3p_1im_dep, # loss_func
        :CICO_ONE_PASS,
        loss_crit = 9.,
        scale = [:log, :log, :log],
        local_alg = :LD_SLSQP,
        loss_grad = :FINITE,
        silent = true
    )
    @test res0.status == :SCAN_BOUND_REACHED
    @test res0.counter < 1000
end

## explicit gradient

@testset ":LD_SLSQP, no scale, Function" begin
    res0 = get_endpoint(
        [3., 2., 2.1],
        1,            # parameter number
        f_3p_1im_dep, # loss_func
        :CICO_ONE_PASS,
        loss_crit = 9.,
        local_alg = :LD_SLSQP,
        loss_grad = (x) -> [  # (x[1]-3.0)^2 + (x[1]-x[2]-1.0)^2 + 0*x[3]^2
            2.0 * (x[1]-3.0) + 2.0 * (x[1]-x[2]-1.0),
            -2.0 * (x[1]-x[2]-1.0),
            0.0
        ],
        silent = true
    )
    @test res0.status == :BORDER_FOUND_BY_SCAN_TOL
    @test res0.counter < 1000
    @test isapprox(res0.value, 5.; atol=1e-2)
end

@testset ":LD_SLSQP, no scale, Function, !unpredictable" begin
    res0 = get_endpoint(
        [3., 2., 2.1],
        3,            # parameter number
        f_3p_1im_dep, # loss_func
        :CICO_ONE_PASS,
        loss_crit = 9.,
        local_alg = :LD_SLSQP,
        loss_grad = (x) -> [  # (x[1]-3.0)^2 + (x[1]-x[2]-1.0)^2 + 0*x[3]^2
            2.0 * (x[1]-3.0) + 2.0 * (x[1]-x[2]-1.0),
            -2.0 * (x[1]-x[2]-1.0),
            0.0
        ],
        silent = true
    )
    @test res0.status == :SCAN_BOUND_REACHED
    @test res0.counter < 1000
end

@testset ":LD_SLSQP, log scale, Function" begin
    res0 = get_endpoint(
        [3., 2., 2.1],
        1,            # parameter number
        f_3p_1im_dep, # loss_func
        :CICO_ONE_PASS,
        loss_crit = 9.,
        scale = [:log, :log, :log],
        local_alg = :LD_SLSQP,
        loss_grad = (x) -> [  # (x[1]-3.0)^2 + (x[1]-x[2]-1.0)^2 + 0*x[3]^2
            2.0 * (x[1]-3.0) + 2.0 * (x[1]-x[2]-1.0),
            -2.0 * (x[1]-x[2]-1.0),
            0.0
        ],
        scan_tol = 1e-4,
        silent = true
    )
    @test res0.status == :BORDER_FOUND_BY_SCAN_TOL
    @test res0.counter < 1000
    @test isapprox(res0.value, 5.; atol=1e-2)
end

@testset ":LD_SLSQP, log scale, Function, !unpredictable" begin
    res0 = get_endpoint(
        [3., 2., 2.1],
        3,            # parameter number
        f_3p_1im_dep, # loss_func
        :CICO_ONE_PASS,
        loss_crit = 9.,
        scale = [:log, :log, :log],
        local_alg = :LD_SLSQP,
        loss_grad = (x) -> [  # (x[1]-3.0)^2 + (x[1]-x[2]-1.0)^2 + 0*x[3]^2
            2.0 * (x[1]-3.0) + 2.0 * (x[1]-x[2]-1.0),
            -2.0 * (x[1]-x[2]-1.0),
            0.0
        ],
        silent = true
    )
    @test res0.status == :SCAN_BOUND_REACHED
    @test res0.counter < 1000
end


text:
# ---
# title: 407. Trapping Rain Water II
# id: problem407
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Heap, Breadth-first Search
# link: <https://leetcode.com/problems/trapping-rain-water-ii/description/>
# hidden: true
# ---
#
# Given an `m x n` matrix of positive integers representing the height of each
# unit cell in a 2D elevation map, compute the volume of water it is able to
# trap after raining.
#
# **Example:**
#
#
#
#     Given the following 3x6 height map:
#     [
#       [1,4,3,1,3,2],
#       [3,2,1,3,2,4],
#       [2,3,3,2,3,1]
#     ]
#
#     Return 4.
#
#
# ![](https://assets.leetcode.com/uploads/2018/10/13/rainwater_empty.png)
#
# The above image represents the elevation map
# `[[1,4,3,1,3,2],[3,2,1,3,2,4],[2,3,3,2,3,1]]` before the rain.
#
#
#
# ![](https://assets.leetcode.com/uploads/2018/10/13/rainwater_fill.png)
#
# After the rain, water is trapped between the blocks. The total volume of water
# trapped is 4.
#
#
#
# **Constraints:**
#
#   * `1 <= m, n <= 110`
#   * `0 <= heightMap[i][j] <= 20000`
#
#
## @lc code=start
using LeetCode

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


text:
"
    Represents a periodic statmech mpo
"
struct InfiniteMPO{O<:MPOTensor} <: Operator
    opp::PeriodicArray{O,1}
end

InfiniteMPO(t::AbstractTensorMap) = InfiniteMPO(fill(t,1));
InfiniteMPO(t::AbstractArray{T,1}) where T<:MPOTensor = InfiniteMPO(PeriodicArray(t));
Base.length(t::InfiniteMPO) = length(t.opp);
Base.size(t::InfiniteMPO) = (length(t),)

Base.getindex(t::InfiniteMPO,i) = getindex(t.opp,i);
Base.eltype(t::InfiniteMPO{O}) where O = O;
Base.iterate(t::InfiniteMPO,itstate = 1) = iterate(t.opp,itstate);

function Base.convert(::Type{InfiniteMPS},mpo::InfiniteMPO)
    InfiniteMPS(map(mpo.opp) do t
        permute(t,(1,2,4),(3,))
    end)
end

function Base.convert(::Type{InfiniteMPO},mps::InfiniteMPS)
    InfiniteMPO(map(mps.AL) do t
        permute(t,(1,2),(4,3))
    end)
end

#naively apply the mpo to the mps
function Base.:*(mpo::InfiniteMPO,st::InfiniteMPS)
    length(st) == length(mpo) || throw(ArgumentError("dimension mismatch"))

    fusers = PeriodicArray(map(zip(st.AL,mpo)) do (al,mp)
        isometry(fuse(_firstspace(al),_firstspace(mp)),_firstspace(al)*_firstspace(mp))
    end)

    InfiniteMPS(map(1:length(st)) do i
        @tensor t[-1 -2;-3] := st.AL[i][1,2,3]*mpo[i][4,-2,5,2]*fusers[i][-1,1,4]*conj(fusers[i+1][-3,3,5])
    end)
end
function Base.:*(mpo::InfiniteMPO,st::FiniteMPS)
    mod(length(mpo),length(st)) == 0 || throw(ArgumentError("dimension mismatch"))

    tensors = [st.AC[1];st.AR[2:end]];
    mpot = mpo[1:length(st)];

    fusers = PeriodicArray(map(zip(tensors,mpot)) do (al,mp)
        isometry(fuse(_firstspace(al),_firstspace(mp)),_firstspace(al)*_firstspace(mp))
    end)

    (_firstspace(mpot[1]) == oneunit(_firstspace(mpot[1])) && space(mpot[end],3)' == _firstspace(mpot[1])) ||
        @warn "mpo does not start/end with a trivial leg"

    FiniteMPS(map(1:length(st)) do i
        @tensor t[-1 -2;-3] := tensors[i][1,2,3]*mpot[i][4,-2,5,2]*fusers[i][-1,1,4]*conj(fusers[i+1][-3,3,5])
    end)
end

function Base.:*(mpo1::InfiniteMPO,mpo2::InfiniteMPO)
    length(mpo1) == length(mpo2) || throw(ArgumentError("dimension mismatch"))

    fusers = PeriodicArray(map(zip(mpo2.opp,mpo1.opp)) do (mp1,mp2)
        isometry(fuse(_firstspace(mp1),_firstspace(mp2)),_firstspace(mp1)*_firstspace(mp2))
    end)

    InfiniteMPO(map(1:length(mpo1)) do i
        @tensor t[-1 -2;-3 -4] := mpo2[i][1,2,3,-4]*mpo1[i][4,-2,5,2]*fusers[i][-1,1,4]*conj(fusers[i+1][-3,3,5])
    end)
end

#---
const MPOMultiline = Multiline{<:InfiniteMPO}

Base.CartesianIndices(t::Union{InfiniteMPO,MPOMultiline}) = CartesianIndices(size(t))
Base.eltype(t::MPOMultiline) = eltype(t[1]);

MPOMultiline(t::AbstractTensorMap) = MPOMultiline(fill(t,1,1));
function MPOMultiline(t::AbstractArray{T,2}) where T<: MPOTensor
    Multiline(PeriodicArray(map(1:size(t,1)) do row
        InfiniteMPO(t[row,:]);
    end))
end

Base.lastindex(t::MPOMultiline) = prod(size(t));
Base.iterate(t::MPOMultiline,i=1) = i <= lastindex(t) ? (t[(i÷end)+1][mod1(i,size(t,1))],i+1) : nothing;
Base.getindex(t::MPOMultiline,i,j) = t[i][j];
Base.getindex(t::MPOMultiline,i::CartesianIndex{2}) = t[i[1],i[2]];
Base.convert(::Type{MPOMultiline},t::InfiniteMPO) = Multiline([t]);
Base.convert(::Type{InfiniteMPO},t::MPOMultiline) = t[1];

function Base.:*(mpo::MPOMultiline,st::MPSMultiline)
    size(mpo) == size(st) || throw(ArgumentError("dimension mismatch"))
    Multiline(map(*,zip(mpo,st)))
end

function Base.:*(mpo1::MPOMultiline,mpo2::MPOMultiline)
    size(mpo1) == size(mpo2) || throw(ArgumentError("dimension mismatch"))
    Multiline(map(*,zip(mpo1,mpo2)));
end


text:
using ArrayViewsAPL, Base.Test

## Higher dimensional fuzz testing
A = reshape(1:625, 5, 5, 5, 5)
function test0(A)
## Indexing to extract a scalar
# Linear indexing
for i = 5:7:400
    S = subview(A, i)
    @test S[1] == A[i]
    S = sliceview(A, i)
    @test S[1] == A[i]
end
# Linear indexing in trailing dimensions
for j = 10:13:125, i = 1:5
    S = subview(A, i, j)
    @test S[1] == A[i, j]
    S = sliceview(A, i, j)
    @test S[1] == A[i, j]
end
for k = 3:2:25, j = 1:5, i = 1:5
    S = subview(A, i, j, k)
    @test S[1] == A[i, j, k]
    S = sliceview(A, i, j, k)
    @test S[1] == A[i, j, k]
end
# Scalar indexing
for l = 1:5, k = 1:5, j = 1:5, i = 1:5
    S = subview(A, i, j, k, l)
    @test S[1] == A[i,j,k,l]
    S = sliceview(A, i, j, k, l)
    @test S[1] == A[i,j,k,l]
end
end
test0(A)

if !isdefined(:eq_arrays_linear) # to make subsequent runs much faster
function eq_arrays_linear(S, B)
    length(S) == length(B) || error("length mismatch")
    alleq = true
    for j = 1:length(B)
        alleq &= S[j] == B[j]
    end
    @test alleq
    nothing
end
end

## Generic indexing
index_choices = (1,3,5,2:3,3:3,1:2:5,4:-1:2,1:5,[1,4,3])
function test1(n)
for i = 1:n
    ii = rand(1:length(index_choices), 4)
    indexes = index_choices[ii]
    B = A[indexes...]
    S = subview(A, indexes...)
    @test size(S) == size(B)
    eq_arrays_linear(S, B)
    S = sliceview(A, indexes...)
    @test size(S) == size(B)[[map(x->!isa(x,Int), indexes)...]]
    eq_arrays_linear(S, B)
end
end
test1(10^4)

index3_choices = (2,7,13,22,5:17,2:3:20,[18,14,9,11])
function test2(n)
indexes = Array(Any, 3)
for i = 1:n
    ii = rand(1:length(index_choices), 2)
    indexes[1] = index_choices[ii[1]]
    indexes[2] = index_choices[ii[2]]
    indexes[3] = index3_choices[rand(1:length(index3_choices))]
    B = A[indexes...]
    S = subview(A, indexes...)
    @test size(S) == size(B)
    eq_arrays_linear(S, B)
    S = sliceview(A, indexes...)
    @test size(S) == size(B)[[map(x->!isa(x,Int), indexes)...]]
    eq_arrays_linear(S, B)
end
end
test2(10^4)

index2_choices = (13,42,55,99,111,88:99,2:5:54,[72,38,37,101])
function test3(n)
indexes = Array(Any, 2)
for i = 1:n
    indexes[1] = index_choices[rand(1:length(index_choices))]
    indexes[2] = index2_choices[rand(1:length(index2_choices))]
    B = A[indexes...]
    S = subview(A, indexes...)
    @test size(S) == size(B)
    eq_arrays_linear(S, B)
    S = sliceview(A, indexes...)
    @test size(S) == size(B)[[map(x->!isa(x,Int), indexes)...]]
    eq_arrays_linear(S, B)
end
end
test3(10^4)

index1_choices = (25,63,128,344,599,121:315,43:51:600,[19,1,603,623,555,229])
function test4()
for indexes in index1_choices
    B = A[indexes]
    S = subview(A, indexes)
    @test size(S) == size(B)
    eq_arrays_linear(S, B)
    S = sliceview(A, indexes)
    @test size(S) == size(B)[[!isa(indexes,Int)]]
    eq_arrays_linear(S, B)
end
end
test4()

function randsubset(indexes1)
    n = length(indexes1)
    indexesnew = Array(Any, n)
    for i = 1:n
        I = indexes1[i]
        if isa(I, Int)
            indexesnew[i] = 1
        else
            k = length(I)
            r = rand()
            if r < 0.3
                l = rand(1:k)
                indexesnew[i] = rand(1:k, l)
            elseif r < 0.6
                indexesnew[i] = rand(1:k)
            else
                i1 = rand(1:k)
                i2 = rand(1:k)
                indexesnew[i] = i1 >= i2 ? (i1:-1:i2) : (i1:1:i2)
            end
        end
    end
    indexes2 = [indexes1[k][indexesnew[k]] for k = 1:n]
    indexesnew, indexes2
end

## Views of views
function test5(A, n)
for i = 1:n
    ii = rand(1:length(index_choices), ndims(A))
    indexes = index_choices[ii]
    indexesnew, indexescomposed = randsubset(indexes)
    B = A[indexescomposed...]
    S1 = subview(A, indexes...)
    S = subview(S1, indexesnew...)
    @test size(S) == size(B)
    eq_arrays_linear(S, B)
    S1 = sliceview(A, indexes...)
    T, N, P, IV = typeof(S1).parameters
    keepdim = [map(x->!(x<:Int), IV)...]
    S = sliceview(S1, indexesnew[keepdim]...)
    @test size(S) == size(B)[[map(x->!isa(x,Int), indexescomposed)...]]
    eq_arrays_linear(S, B)
end
end
# Test this one in just 3d because of the combinatorial explosion from views-of-views
A = reshape(1:125, 5, 5, 5)
test5(A, 10^5)  # peaks out at more than 2000 separate stagedfunction generation events


text:
@testset "1089.duplicate-zeros.jl" begin
    arr = [1,0,2,3,0,4,5,0]
    duplicate_zeros!(arr)
    @test arr == [1,0,0,2,3,0,0,4]

    arr2 = [1,2,3]
    duplicate_zeros!(arr2)
    @test arr2 == [1,2,3]
end


text:
using Documenter, LinearSVM

makedocs(modules=[LinearSVM], doctest=true)
deploydocs(deps = Deps.pip("mkdocs", "python-markdown-math"),
repo   = "https://github.com/rlefmann/LinearSVM.jl.git",
julia  = "0.4.5",
osname = "linux")


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMQuad.jl/blob/master/LICENSE

### Gauss quadrature rules for tetrahedrons

""" Gauss-Legendre quadrature, 1 point rule on tetrahedron. """
function get_quadrature_points(::Type{Val{:GLTET1}})
    weights = (1.0/6.0, )
    points = ((1.0/4.0, 1.0/4.0, 1.0/4.0), )
    return zip(weights, points)
end

function get_order(::Type{Val{:GLTET1}})
    return 1
end

""" Gauss-Legendre quadrature, 4 point rule on tetrahedron. """
function get_quadrature_points(::Type{Val{:GLTET4}})
    a = (5.0+3.0*sqrt(5.0))/20.0
    b = (5.0-sqrt(5.0))/20.0
    w = 1.0/24.0
    weights = (w, w, w, w)
    points = ((a, b, b), (b, a, b), (b, b, a), (b, b, b))
    return zip(weights, points)
end

function get_order(::Type{Val{:GLTET4}})
    return 2
end

""" Gauss-Legendre quadrature, 5 point rule on tetrahedron. """
function get_quadrature_points(::Type{Val{:GLTET5}})
    a = 1.0/4.0
    b = 1.0/6.0
    c = 1.0/2.0
    weights = (-2.0/15.0, 3.0/40.0, 3.0/40.0, 3.0/40.0, 3.0/40.0)
    points = ((a, a, a), (b, b, b), (b, b, c), (b, c, b), (c, b, b))
    return zip(weights, points)
end

function get_order(::Type{Val{:GLTET5}})
    return 3
end

""" Gauss-Legendre quadrature, 15 point rule on tetrahedron. """
function get_quadrature_points(::Type{Val{:GLTET15}})
    a = 1.0/4.0
    b1 = 1.0/34.0*(7.0 + sqrt(15.0))
    b2 = 1.0/34.0*(7.0 - sqrt(15.0))
    c1 = 1.0/34.0*(13.0 - 3.0*sqrt(15.0))
    c2 = 1.0/34.0*(13.0 + 3.0*sqrt(15.0))
    d = 1.0/20.0*(5.0 - sqrt(15.0))
    f = 1.0/20.0*(5.0 + sqrt(15.0))
    w1 = 8.0/405.0
    w2 = (2665.0 - 14.0*sqrt(15.0))/226800.0
    w3 = (2665.0 + 14.0*sqrt(15.0))/226800.0
    w4 = 5.0/567.0
    weights = (w1, w2, w2, w2, w2, w3, w3, w3, w3, w4, w4, w4, w4, w4, w4)
    points = ((a, a, a), (b1, b1, b1), (b1, b1, c1), (b1, c1, b1), (c1, b1, b1),
              (b2, b2, b2), (b2, b2, c2), (b2, c2, b2), (c2, b2, b2), (d, d, f),
              (d, f, d), (f, d, d), (d, f, f), (f, d, f), (f, f, d))
    return zip(weights, points)
end

function get_order(::Type{Val{:GLTET15}})
    return 4
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using JuliaFEM
using JuliaFEM.Preprocess
using JuliaFEM.Postprocess
using JuliaFEM.Testing

@testset "3d upper side curved contact" begin

    # TODO: accurate solution is not known, verify using another fem software
    # however results look very meaningful and probably this is right.

    meshfile = Pkg.dir("JuliaFEM") * "/test/testdata/block_3d_curved.med"
    mesh = aster_read_mesh(meshfile)

    upper = Problem(Elasticity, "upper", 3)
    upper.elements = create_elements(mesh, "UPPER")
    update!(upper, "youngs modulus", 96.0)
    update!(upper, "poissons ratio", 1/3)

    lower = Problem(Elasticity, "lower", 3)
    lower.elements = create_elements(mesh, "LOWER")
    update!(lower, "youngs modulus", 96.0)
    update!(lower, "poissons ratio", 1/3)

    bc_upper = Problem(Dirichlet, "upper boundary", 3, "displacement")
    bc_upper.elements = create_elements(mesh, "UPPER_TOP")
    update!(bc_upper, "displacement 1",  0.0)
    update!(bc_upper, "displacement 2",  0.0)
    update!(bc_upper, "displacement 3", -0.1)

    bc_lower = Problem(Dirichlet, "lower boundary", 3, "displacement")
    bc_lower.elements = create_elements(mesh, "LOWER_BOTTOM")
    update!(bc_lower, "displacement 1", 0.0)
    update!(bc_lower, "displacement 2", 0.0)
    update!(bc_lower, "displacement 3", 0.0)

    contact = Problem(Contact, "contact between upper and lower block", 3, "displacement")
    contact_slave_elements = create_elements(mesh, "LOWER_TOP")
    contact_master_elements = create_elements(mesh, "UPPER_BOTTOM")
    update!(contact_slave_elements, "master elements", contact_master_elements)
    contact.elements = [contact_master_elements; contact_slave_elements]

    solver = Solver(Nonlinear)
    push!(solver, upper, lower, bc_upper, bc_lower, contact)

    solver()
    for element in get_slave_elements(contact)
        normal = element("normal", [1/3, 1/3], solver.time)
        @test isapprox(normal, [0.0, 0.0, 1.0])
        pres = dot(normal, element("reaction force", [1/3, 1/3], solver.time))
        info("pressure = $pres")
        #info(element("displacement", [1/3, 1/3], solver.time))
    end
    normu = norm(contact.assembly.u)
    info("displacement field norm = $normu")
    @test isapprox(normu, 0.7417557629004985)

end