text:
const VI = MOI.VariableIndex

struct HS071 <: MOI.AbstractNLPEvaluator
    enable_hessian::Bool
end

# hs071
# min x1 * x4 * (x1 + x2 + x3) + x3
# st  x1 * x2 * x3 * x4 >= 25
#     x1^2 + x2^2 + x3^2 + x4^2 = 40
#     1 <= x1, x2, x3, x4 <= 5
# Start at (1,5,5,1)
# End at (1.000..., 4.743..., 3.821..., 1.379...)

function MOI.initialize(d::HS071, requested_features::Vector{Symbol})
    for feat in requested_features
        if !(feat in MOI.features_available(d))
            error("Unsupported feature $feat")
            # TODO: implement Jac-vec and Hess-vec products
            # for solvers that need them
        end
    end
end

function MOI.features_available(d::HS071)
    if d.enable_hessian
        return [:Grad, :Jac, :Hess, :ExprGraph]
    else
        return [:Grad, :Jac, :ExprGraph]
    end
end

MOI.objective_expr(d::HS071) = :(x[$(VI(1))] * x[$(VI(4))] * (x[$(VI(1))] +
                                 x[$(VI(2))] + x[$(VI(3))]) + x[$(VI(3))])

function MOI.constraint_expr(d::HS071, i::Int)
    if i == 1
        return :(x[$(VI(1))] * x[$(VI(2))] * x[$(VI(3))] * x[$(VI(4))] >= 25.0)
    elseif i == 2
        return :(x[$(VI(1))]^2 + x[$(VI(2))]^2 +
                 x[$(VI(3))]^2 + x[$(VI(4))]^2 == 40.0)
    else
        error("Out of bounds constraint.")
    end
end

MOI.eval_objective(d::HS071, x) = x[1] * x[4] * (x[1] + x[2] + x[3]) + x[3]

function MOI.eval_constraint(d::HS071, g, x)
    g[1] = x[1]   * x[2]   * x[3]   * x[4]
    g[2] = x[1]^2 + x[2]^2 + x[3]^2 + x[4]^2
end

function MOI.eval_objective_gradient(d::HS071, grad_f, x)
    grad_f[1] = x[1] * x[4] + x[4] * (x[1] + x[2] + x[3])
    grad_f[2] = x[1] * x[4]
    grad_f[3] = x[1] * x[4] + 1
    grad_f[4] = x[1] * (x[1] + x[2] + x[3])
end

function MOI.jacobian_structure(d::HS071)
    return Tuple{Int64,Int64}[(1,1), (1,2), (1,3), (1,4), (2,1), (2,2),
                              (2,3), (2,4)]
end
# lower triangle only
function MOI.hessian_lagrangian_structure(d::HS071)
    @assert d.enable_hessian
    return Tuple{Int64,Int64}[(1,1), (2,1), (2,2), (3,1), (3,2), (3,3),
                              (4,1), (4,2), (4,3), (4,4)]
end

function MOI.eval_constraint_jacobian(d::HS071, J, x)
    # Constraint (row) 1
    J[1] = x[2]*x[3]*x[4]  # 1,1
    J[2] = x[1]*x[3]*x[4]  # 1,2
    J[3] = x[1]*x[2]*x[4]  # 1,3
    J[4] = x[1]*x[2]*x[3]  # 1,4
    # Constraint (row) 2
    J[5] = 2*x[1]  # 2,1
    J[6] = 2*x[2]  # 2,2
    J[7] = 2*x[3]  # 2,3
    J[8] = 2*x[4]  # 2,4
end

function MOI.eval_hessian_lagrangian(d::HS071, H, x, σ, μ)
    @assert d.enable_hessian
    # Again, only lower left triangle
    # Objective
    H[1] = σ * (2*x[4])               # 1,1
    H[2] = σ * (  x[4])               # 2,1
    H[3] = 0                          # 2,2
    H[4] = σ * (  x[4])               # 3,1
    H[5] = 0                          # 3,2
    H[6] = 0                          # 3,3
    H[7] = σ* (2*x[1] + x[2] + x[3])  # 4,1
    H[8] = σ * (  x[1])               # 4,2
    H[9] = σ * (  x[1])               # 4,3
    H[10] = 0                         # 4,4

    # First constraint
    H[2] += μ[1] * (x[3] * x[4])  # 2,1
    H[4] += μ[1] * (x[2] * x[4])  # 3,1
    H[5] += μ[1] * (x[1] * x[4])  # 3,2
    H[7] += μ[1] * (x[2] * x[3])  # 4,1
    H[8] += μ[1] * (x[1] * x[3])  # 4,2
    H[9] += μ[1] * (x[1] * x[2])  # 4,3

    # Second constraint
    H[1]  += μ[2] * 2  # 1,1
    H[3]  += μ[2] * 2  # 2,2
    H[6]  += μ[2] * 2  # 3,3
    H[10] += μ[2] * 2  # 4,4

end

function hs071test_template(model::MOI.ModelLike, config::TestConfig, evaluator::HS071)
    atol = config.atol
    rtol = config.rtol

    @test MOI.supports(model, MOI.NLPBlock())
    @test MOI.supports_constraint(model, MOI.SingleVariable, MOI.LessThan{Float64})
    @test MOI.supports_constraint(model, MOI.SingleVariable, MOI.GreaterThan{Float64})
    @test MOI.supports(model, MOI.VariablePrimalStart(), MOI.VariableIndex)

    MOI.empty!(model)
    @test MOI.is_empty(model)

    lb = [25.0, 40.0]
    ub = [Inf, 40.0]

    block_data = MOI.NLPBlockData(MOI.NLPBoundsPair.(lb, ub), evaluator, true)

    v = MOI.add_variables(model, 4)
    @test MOI.get(model, MOI.NumberOfVariables()) == 4

    l = [1.0,1.0,1.0,1.0]
    u = [5.0,5.0,5.0,5.0]
    start = [1,5,5,1]

    for i in 1:4
        cub = MOI.add_constraint(model, MOI.SingleVariable(v[i]), MOI.LessThan(u[i]))
        # We test this after the creation of every `SingleVariable` constraint
        # to ensure a good coverage of corner cases.
        @test cub.value == v[i].value
        clb = MOI.add_constraint(model, MOI.SingleVariable(v[i]), MOI.GreaterThan(l[i]))
        @test clb.value == v[i].value
        MOI.set(model, MOI.VariablePrimalStart(), v[i], start[i])
    end

    MOI.set(model, MOI.NLPBlock(), block_data)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    # TODO: config.query tests
    if config.solve
        MOI.optimize!(model)

        @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status

        @test MOI.get(model, MOI.ResultCount()) >= 1

        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 17.014017145179164 atol=atol rtol=rtol

        optimal_v = [1.0, 4.7429996418092970, 3.8211499817883077, 1.3794082897556983]

        @test MOI.get(model, MOI.VariablePrimal(), v) ≈ optimal_v atol=atol rtol=rtol

        # TODO: Duals? Maybe better to test on a convex instance.
    end
end

hs071_test(model, config) = hs071test_template(model, config, HS071(true))
hs071_no_hessian_test(model, config) = hs071test_template(model, config, HS071(false))

# Test for FEASIBILITY_SENSE.
# Find x satisfying x^2 == 1.
struct FeasibilitySenseEvaluator <: MOI.AbstractNLPEvaluator
    enable_hessian::Bool
end

function MOI.initialize(d::FeasibilitySenseEvaluator,
                         requested_features::Vector{Symbol})
    for feat in requested_features
        if !(feat in MOI.features_available(d))
            error("Unsupported feature $feat")
            # TODO: implement Jac-vec and Hess-vec products
            # for solvers that need them
        end
    end
end

function MOI.features_available(d::FeasibilitySenseEvaluator)
    if d.enable_hessian
        return [:Grad, :Jac, :Hess, :ExprGraph]
    else
        return [:Grad, :Jac, :ExprGraph]
    end
end

MOI.objective_expr(d::FeasibilitySenseEvaluator) = :()

function MOI.constraint_expr(d::FeasibilitySenseEvaluator, i::Int)
    if i == 1
        return :(x[$(VI(1))]^2 == 1)
    else
        error("Out of bounds constraint.")
    end
end

MOI.eval_objective(d::FeasibilitySenseEvaluator, x) = 0.0

function MOI.eval_constraint(d::FeasibilitySenseEvaluator, g, x)
    g[1] = x[1]^2
end

function MOI.eval_objective_gradient(d::FeasibilitySenseEvaluator, grad_f, x)
    grad_f[1] = 0.0
end

function MOI.jacobian_structure(d::FeasibilitySenseEvaluator)
    return Tuple{Int64,Int64}[(1,1)]
end

function MOI.hessian_lagrangian_structure(d::FeasibilitySenseEvaluator)
    @assert d.enable_hessian
    return Tuple{Int64,Int64}[(1,1)]
end

function MOI.eval_constraint_jacobian(d::FeasibilitySenseEvaluator, J, x)
    J[1] = 2x[1]
end

function MOI.eval_hessian_lagrangian(d::FeasibilitySenseEvaluator, H, x, σ, μ)
    @assert d.enable_hessian
    H[1] = 2μ[1] # 1,1
end

function feasibility_sense_test_template(model::MOI.ModelLike,
                                         config::TestConfig,
                                         set_has_objective::Bool,
                                         evaluator::FeasibilitySenseEvaluator)
    atol = config.atol
    rtol = config.rtol

    @test MOI.supports(model, MOI.NLPBlock())
    @test MOI.supports(model, MOI.VariablePrimalStart(), MOI.VariableIndex)

    MOI.empty!(model)
    @test MOI.is_empty(model)

    lb = [1.0]
    ub = [1.0]

    block_data = MOI.NLPBlockData(MOI.NLPBoundsPair.(lb, ub), evaluator,
                                  set_has_objective)

    x = MOI.add_variable(model)
    @test MOI.get(model, MOI.NumberOfVariables()) == 1

    # Avoid starting at zero because it's a critial point.
    MOI.set(model, MOI.VariablePrimalStart(), x, 1.5)

    MOI.set(model, MOI.NLPBlock(), block_data)
    MOI.set(model, MOI.ObjectiveSense(), MOI.FEASIBILITY_SENSE)

    # TODO: config.query tests
    if config.solve
        MOI.optimize!(model)

        @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status

        @test MOI.get(model, MOI.ResultCount()) >= 1

        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 0.0 atol=atol rtol=rtol

        @test abs(MOI.get(model, MOI.VariablePrimal(), x)) ≈ 1.0 atol=atol rtol=rtol
    end
end

function feasibility_sense_with_objective_and_hessian_test(model, config)
    feasibility_sense_test_template(model, config, true,
                                    FeasibilitySenseEvaluator(true))
end

function feasibility_sense_with_objective_and_no_hessian_test(model, config)
    feasibility_sense_test_template(model, config, true,
                                    FeasibilitySenseEvaluator(false))
end

function feasibility_sense_with_no_objective_and_with_hessian_test(model, config)
    feasibility_sense_test_template(model, config, false,
                                    FeasibilitySenseEvaluator(true))
end

function feasibility_sense_with_no_objective_and_no_hessian_test(model, config)
    feasibility_sense_test_template(model, config, false,
                                    FeasibilitySenseEvaluator(false))
end

function nlp_objective_and_moi_objective_test(model::MOI.ModelLike,
                                              config::TestConfig)
    atol = config.atol
    rtol = config.rtol

    @test MOI.supports(model, MOI.NLPBlock())
    @test MOI.supports(model, MOI.VariablePrimalStart(), MOI.VariableIndex)

    MOI.empty!(model)
    @test MOI.is_empty(model)

    lb = [1.0]
    ub = [2.0]

    block_data = MOI.NLPBlockData(MOI.NLPBoundsPair.(lb, ub),
                                  FeasibilitySenseEvaluator(true), true)

    x = MOI.add_variable(model)
    @test MOI.get(model, MOI.NumberOfVariables()) == 1

    # Avoid starting at zero because it's a critial point.
    MOI.set(model, MOI.VariablePrimalStart(), x, 1.5)

    MOI.set(model, MOI.NLPBlock(), block_data)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    f_x =  MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x)], 0.0)
    # This objective function should be ignored.
    MOI.set(model, MOI.ObjectiveFunction{typeof(f_x)}(), f_x)

    # TODO: config.query tests
    if config.solve
        MOI.optimize!(model)

        @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status

        @test MOI.get(model, MOI.ResultCount()) >= 1

        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 0.0 atol=atol rtol=rtol

        @test 1.0 ≤ MOI.get(model, MOI.VariablePrimal(), x) ≤ 2.0
    end
end

const nlptests = Dict("hs071" => hs071_test,
                      "hs071_no_hessian" => hs071_no_hessian_test,
                      "feasibility_sense_with_objective_and_hessian" =>
                      feasibility_sense_with_objective_and_hessian_test,
                      "feasibility_sense_with_objective_and_no_hessian" =>
                      feasibility_sense_with_objective_and_no_hessian_test,
                      "feasibility_sense_with_no_objective_and_with_hessian" =>
                      feasibility_sense_with_no_objective_and_with_hessian_test,
                      "feasibility_sense_with_no_objective_and_no_hessian" =>
                      feasibility_sense_with_no_objective_and_no_hessian_test,
                      "nlp_objective_and_moi_objective" =>
                      nlp_objective_and_moi_objective_test
                      )

@moitestset nlp

"""
    test_linear_mixed_complementarity(model::MOI.ModelLike, config::TestConfig)

Test the solution of the linear mixed-complementarity problem:
`F(x) complements x`, where `F(x) = M * x .+ q` and `0 <= x <= 10`.
"""
function test_linear_mixed_complementarity(model::MOI.ModelLike, config::TestConfig)
    MOI.empty!(model)
    x = MOI.add_variables(model, 4)
    MOI.add_constraint.(model, MOI.SingleVariable.(x), MOI.Interval(0.0, 10.0))
    MOI.set.(model, MOI.VariablePrimalStart(), x, 0.0)
    M = Float64[0  0 -1 -1; 0  0  1 -2; 1 -1  2 -2; 1  2 -2  4]
    q = [2; 2; -2; -6]
    terms = MOI.VectorAffineTerm{Float64}[]
    for i = 1:4
        push!(
            terms,
            MOI.VectorAffineTerm(4 + i, MOI.ScalarAffineTerm(1.0, x[i]))
        )
        for j = 1:4
            iszero(M[i, j]) && continue
            push!(
                terms,
                MOI.VectorAffineTerm(i, MOI.ScalarAffineTerm(M[i, j], x[j]))
            )
        end
    end
    MOI.add_constraint(
        model,
        MOI.VectorAffineFunction(terms, [q; 0.0; 0.0; 0.0; 0.0]),
        MOI.Complements(4)
    )
    @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED
    if config.solve
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status
        x_val = MOI.get.(model, MOI.VariablePrimal(), x)
        @test isapprox(
            x_val, [2.8, 0.0, 0.8, 1.2], atol = config.atol, rtol = config.rtol
        )
    end
end

const mixed_complementaritytests = Dict(
    "test_linear_mixed_complementarity" => test_linear_mixed_complementarity,
)

@moitestset mixed_complementarity

"""
    test_qp_complementarity_constraint(model::MOI.ModelLike, config::TestConfig)

Test the solution of the quadratic program with complementarity constraints:

```
    min  (x0 - 5)^2 +(2 x1 + 1)^2
    s.t.  -1.5 x0 + 2 x1 + x2 - 0.5 x3 + x4 = 2
        x2 complements(3 x0 - x1 - 3)
        x3 complements(-x0 + 0.5 x1 + 4)
        x4 complements(-x0 - x1 + 7)
        x0, x1, x2, x3, x4 >= 0

```
which rewrites, with auxiliary variables

```
    min  (x0 - 5)^2 +(2 x1 + 1)^2
    s.t.  -1.5 x0 + 2 x1 + x2 - 0.5 x3 + x4 = 2  (cf1)
        3 x0 - x1 - 3 - x5 = 0                 (cf2)
        -x0 + 0.5 x1 + 4 - x6 = 0              (cf3)
        -x0 - x1 + 7 - x7 = 0                  (cf4)
        x2 complements x5
        x3 complements x6
        x4 complements x7
        x0, x1, x2, x3, x4, x5, x6, x7 >= 0

```
"""
function test_qp_complementarity_constraint(
    model::MOI.ModelLike, config::TestConfig
)
    MOI.empty!(model)
    x = MOI.add_variables(model, 8)
    MOI.set.(model, MOI.VariablePrimalStart(), x, 0.0)
    MOI.add_constraint.(model, MOI.SingleVariable.(x), MOI.GreaterThan(0.0))
    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}(),
        MOI.ScalarQuadraticFunction(
            MOI.ScalarAffineTerm.([-10.0, 4.0], x[[1, 2]]),
            MOI.ScalarQuadraticTerm.([2.0, 8.0], x[1:2], x[1:2]),
            26.0
        )
    )
    MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([-1.5, 2.0, 1.0, 0.5, 1.0], x[1:5]), 0.0
        ),
        MOI.EqualTo(2.0)
   )
    MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([3.0, -1.0, -1.0], x[[1, 2, 6]]), 0.0
        ),
        MOI.EqualTo(3.0)
   )
    MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([-1.0, 0.5, -1.0], x[[1, 2, 7]]), 0.0
        ),
        MOI.EqualTo(-4.0)
   )
    MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction(
           MOI.ScalarAffineTerm.(-1.0, x[[1, 2, 8]]), 0.0
        ),
        MOI.EqualTo(-7.0)
   )
    MOI.add_constraint(
        model,
        MOI.VectorOfVariables([x[3], x[4], x[5], x[6], x[7], x[8]]),
        MOI.Complements(3)
    )
    @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED
    if config.solve
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status
        x_val = MOI.get.(model, MOI.VariablePrimal(), x)
        @test isapprox(
            x_val,
            [1.0, 0.0, 3.5, 0.0, 0.0, 0.0, 3.0, 6.0],
            atol = config.atol,
            rtol = config.rtol
        )
        @test isapprox(
            MOI.get(model, MOI.ObjectiveValue()),
            17.0,
            atol = config.atol,
            rtol = config.rtol
        )
    end
end

const math_program_complementarity_constraintstests = Dict(
    "test_qp_complementarity_constraint" => test_qp_complementarity_constraint,
)

@moitestset math_program_complementarity_constraints


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

mutable struct IoK8sApiStorageV1beta1CSINodeDriver <: SwaggerModel
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    nodeID::Any # spec type: Union{ Nothing, String } # spec name: nodeID
    topologyKeys::Any # spec type: Union{ Nothing, Vector{String} } # spec name: topologyKeys

    function IoK8sApiStorageV1beta1CSINodeDriver(;name=nothing, nodeID=nothing, topologyKeys=nothing)
        o = new()
        validate_property(IoK8sApiStorageV1beta1CSINodeDriver, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(IoK8sApiStorageV1beta1CSINodeDriver, Symbol("nodeID"), nodeID)
        setfield!(o, Symbol("nodeID"), nodeID)
        validate_property(IoK8sApiStorageV1beta1CSINodeDriver, Symbol("topologyKeys"), topologyKeys)
        setfield!(o, Symbol("topologyKeys"), topologyKeys)
        o
    end
end # type IoK8sApiStorageV1beta1CSINodeDriver

const _property_map_IoK8sApiStorageV1beta1CSINodeDriver = Dict{Symbol,Symbol}(Symbol("name")=>Symbol("name"), Symbol("nodeID")=>Symbol("nodeID"), Symbol("topologyKeys")=>Symbol("topologyKeys"))
const _property_types_IoK8sApiStorageV1beta1CSINodeDriver = Dict{Symbol,String}(Symbol("name")=>"String", Symbol("nodeID")=>"String", Symbol("topologyKeys")=>"Vector{String}")
Base.propertynames(::Type{ IoK8sApiStorageV1beta1CSINodeDriver }) = collect(keys(_property_map_IoK8sApiStorageV1beta1CSINodeDriver))
Swagger.property_type(::Type{ IoK8sApiStorageV1beta1CSINodeDriver }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiStorageV1beta1CSINodeDriver[name]))}
Swagger.field_name(::Type{ IoK8sApiStorageV1beta1CSINodeDriver }, property_name::Symbol) =  _property_map_IoK8sApiStorageV1beta1CSINodeDriver[property_name]

function check_required(o::IoK8sApiStorageV1beta1CSINodeDriver)
    (getproperty(o, Symbol("name")) === nothing) && (return false)
    (getproperty(o, Symbol("nodeID")) === nothing) && (return false)
    true
end

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

ext: jl


text:
"""
Computes convex hull of 2D points. Assumes non-degenerate inputs.
"""
module ConvexHull

const Point{T} = Tuple{T, T}

raw"""
Classifies whether a point $r$ lies to the left or right of the oriented line $pq$.

Returns:

* -1 if $r$ lies to the left
* 0 if $r$ lies on the line
* 1 if $r$ lies to the right

**Warning**: This computation is not robust. See [Robust Predicates for Computational
Geometry][robust] for details. (If `BigFloat` is used, there is no rounding problems.)

[robust]: http://www.cs.cmu.edu/%7Equake/robust.html
"""
function orient(p::Point{T},
                q::Point{T},
                r::Point{T})::T where T<:AbstractFloat
    pr_x = r[1] - p[1]
    pr_y = r[2] - p[2]
    rq_x = q[1] - r[1]
    rq_y = q[2] - r[2]
    sign(pr_x * rq_y - pr_y * rq_x)
end

# Graham scan
#
# Computational Geometry: Algorithms and Applications, Chapter 1.
function convex_hull_graham_scan!(points::Vector{Point{T}})::Vector{Point{T}} where T<:AbstractFloat
    @assert length(points) >= 3 "Too few points"
    sort!(points)  # in lexical order

    function build_convex_hull(points::AbstractVector{Point{T}})
        n = length(points)
        hull = points[1:2]
        sizehint!(hull, n)

        for i in 3:n
            while true
                if orient(hull[end - 1], hull[end], points[i]) > 0
                    break
                else
                    pop!(hull)
                    if length(hull) < 2
                        break
                    end
                end
            end
            push!(hull, points[i])
        end

        hull
    end

    upper_hull = build_convex_hull(points)
    pop!(upper_hull)
    lower_hull = build_convex_hull(@view points[end:-1:1])
    pop!(lower_hull)

    [upper_hull; lower_hull]
end

# Unsorted version, complexity: O(n^2)
function convex_hull_incremental(points::Vector{Point{T}})::Vector{Point{T}} where T<:AbstractFloat
    n = length(points)
    @assert n >= 3 "Too few points"
    hull = points[1:3]
    if orient(hull[1], hull[2], hull[3]) < 0
        hull[[1, 2]] .= hull[[2, 1]]
    end
    sizehint!(hull, n)

    for i in 4:n
        p = points[i]
        # Vertices inside `[start, stop)` will be removed.
        start::Union{Nothing, Int} = nothing
        stop::Union{Nothing, Int} = nothing
        h_len = length(hull)
        prev_ori = orient(hull[end - 1], hull[end], p) >= 0

        # trying to find two tangent lines
        for (from, to) in ((mod1(i - 1, h_len), i) for i in 1:h_len)
            curr_ori = orient(hull[from], hull[to], p) >= 0
            if !curr_ori && prev_ori
                start = to  # inclusive
                if stop != nothing
                    break
                end
            elseif curr_ori && !prev_ori
                stop = from  # exclusive
                if start != nothing
                    break
                end
            end
            prev_ori = curr_ori
        end

        # Always being on the right side, `p` must be inside the hull.
        if start == nothing
            continue
        end
        # If `start != nothing`, then `stop != nothing`.

        # Note: `start == stop` means no vertices are removed.
        if start <= stop
            splice!(hull, start:(stop - 1), [p])
        else
            # The same as
            #
            #     hull = [hull[stop:(start - 1)]; p]
            #
            # but without allocating extra memory.
            hull[start] = p
            resize!(hull, start)
            deleteat!(hull, 1:(stop - 1))
        end
    end
    hull
end

#===== Unit Tests =====#
using Test
using Random: seed!

function convex_hull_length(hull::Vector{Point{T}})::T where T<:AbstractFloat
    hull2 = circshift(hull, -1)
    map((x, y) -> hypot((x .- y)...), hull2, hull) |> sum
end

@testset "Convex Hull" begin
    seed!()
    n = 100
    a = rand(n, 2)
    points = [tuple(a[i, :]...) for i in 1:n]

    cvh_increm = convex_hull_incremental(points)
    cvh_graham = convex_hull_graham_scan!(points)
    @test convex_hull_length(cvh_increm) ≈ convex_hull_length(cvh_graham)
end
end # module

max_stars_repo_name: gyk/TrivialSolutions

max_issues_repo_name: gyk/TrivialSolutions

max_forks_repo_name: gyk/TrivialSolutions


text:
# TODO: fix and test smoothing
module Simulations2DSymmetric

using ..Curlcurls
using ..Dispersions
using ..Domains
using ..Laplacians
using ..Lattices
using ..Points
using ..SelfEnergies
using ..Shapes
using ..VectorFields
# using LinearAlgebra
using SparseArrays
using Statistics

import ..AbstractDomain
import ..Symmetric, ..Unsymmetric
import ..Simulation
import ..smooth_dielectric!
import ..smooth_pump!
import ..which_domains
import ..simulation_dielectric!
import ..simulation_pump!

function Simulation(
			ω₀::Real,
			lattice_domains::Tuple{LatticeDomain{2,Symmetric,Cartesian}},
			nondispersive_domains::NTuple{L,NondispersiveDomain{2}},
			dispersive_domains::NTuple{M,DispersiveDomain{2}};
			k₃₀::Real=0
			) where {L,M}

	lattice_domain = lattice_domains[1] # there is only one lattice by construction, and it's Cartesian
	smoothed = Ref(false)

	lattice = lattice_domain.lattice
	indices = lattice_domain.indices
	imin, imax = extrema(map(ld->ld[1],indices))
	jmin, jmax = extrema(map(ld->ld[2],indices))
	x = [lattice[i,0] for i ∈ imin:imax]
	y = [lattice[0,j] for j ∈ jmin:jmax]

	lattice_domain_indices = fill(1,length(lattice_domain.x))
	nondispersive_domain_indices = which_domains(nondispersive_domains,lattice_domain.x)
	dispersive_domain_indices = which_domains(dispersive_domains,lattice_domain.x)

	# populate dielectric
	ε = Vector{ComplexF64}(undef,length(lattice_domain.x))
	if isempty(nondispersive_domains)
		for i ∈ eachindex(ε) ε[i] = lattice_domain.ε end
	else
		dielectrics = map(n -> getfield(n,:dielectric), nondispersive_domains)
		for i ∈ eachindex(ε)
			d = nondispersive_domain_indices[i]
			ε[i] = iszero(d) ? lattice_domain.ε : dielectrics[d](lattice_domain.x[i])
		end
	end

	# populate pump and dispersive susceptability
	F = Vector{Float64}(undef,length(lattice_domain.x))
	χ = Vector{AbstractDispersion}(undef,length(lattice_domain.x))
	Fs = Vector{Vector{Float64}}(undef,length(dispersive_domains))
	for i ∈ eachindex(Fs) Fs[i] = zeros(Float64,length(lattice_domain.x)) end
	if isempty(dispersive_domains)
		for i ∈ eachindex(F) F[i] = 0 end
		for i ∈ eachindex(F) χ[i] = NoDispersion() end
	else
		pumps = map(n -> n.pump, dispersive_domains)
		χs = map(n -> n.χ, dispersive_domains)
		for i ∈ eachindex(F)
			d = dispersive_domain_indices[i]
			if iszero(d)
				F[i] = 0
				χ[i] = NoDispersion()
			else
				F[i] = pumps[d](lattice_domain.x[i])
				χ[i] = χs[d]
				Fs[d][i] = pumps[d](lattice_domain.x[i])
			end
		end
	end

	# generate boundary layers
	σx, _ = boundary_layer(lattice_domain, x)
	_, σy = boundary_layer(lattice_domain, y)
	αx, αy = 1 .+ 1im*σx/sqrt(ω₀^2-k₃₀^2), 1 .+ 1im*σy/sqrt(ω₀^2-k₃₀^2)

	# generate half sites and boundary layers
	x_half, y_half = _generate_half_xy(lattice_domain, x)
	σx_half, _ = boundary_layer(lattice_domain, x_half)
	_, σy_half = boundary_layer(lattice_domain, y_half)
	αx_half, αy_half = 1 .+ 1im*σx_half/sqrt(ω₀^2-k₃₀^2), 1 .+ 1im*σy_half/sqrt(ω₀^2-k₃₀^2)

	laplacian = Laplacian{Symmetric}(lattice_domain.lattice, αx, αy, αx_half, αy_half)
	# curlcurl = Curlcurl{Symmetric}(lattice_domain.lattice,α[1],α_half[1],nnm,nnp,indices,interior,surface)

	Σ = SelfEnergy{Symmetric}(lattice_domain, αx_half, αy_half)

	return Simulation{2,Symmetric,ComplexF64,typeof(lattice_domains),typeof(nondispersive_domains),typeof(dispersive_domains),typeof(Σ)}(
		lattice_domains,
		nondispersive_domains,
		dispersive_domains,
		lattice_domain.x,
		lattice_domain_indices,
		nondispersive_domain_indices,
		dispersive_domain_indices,
		ε,
		F,
		χ,
		Fs,
		laplacian,
		# curlcurl,
		Σ,
		(αx, αy),
		(σx, σy),
		(x_half, y_half),
		(σx_half, σy_half),
		ω₀,
		NaN,
		NaN,
		k₃₀,
		smoothed)
end

################################################################################
# SIMULATION{2} building utilities

# used in building Simulation{2}
function _generate_half_xy(lattice_domain::LatticeDomain{2,Symmetric,Cartesian}, x::Vector{TP}) where TP<:Point{2}
	lattice = lattice_domain.lattice
	indices = lattice_domain.indices
	imin, imax = extrema(map(ld->ld[1],indices))
	jmin, jmax = extrema(map(ld->ld[2],indices))
	half_x = [lattice[i,0] for i ∈ ((imin:(imax+1)).-1/2)]
	half_y = [lattice[0,j] for j ∈ ((jmin:(jmax+1)).-1/2)]
	return half_x, half_y
end

@inline function boundary_layer(domain::LatticeDomain{2,Symmetric,Cartesian}, x::Vector{Point{2,Cartesian}})
	σx = domain.boundary.bls[1].(x) .+ domain.boundary.bls[2].(x)
	σy = domain.boundary.bls[3].(x) .+ domain.boundary.bls[4].(x)
	return σx, σy
end

################################################################################
# 2-d smoothing functions

function smooth_dielectric!(sim::Simulation{2,Symmetric}, num_sub_pixel::Integer=NUM_SUBPIXELS)
    indices = sim.nondispersive_domain_indices
    lattice = sim.lattice
    X = Matrix{Point{2,Cartesian}}(undef,num_sub_pixel,num_sub_pixel)
    E = Vector{ComplexF64}(undef,num_sub_pixel)
	r = LinRange(-.5,.5,num_sub_pixel)
    for i ∈ eachindex(indices)
		if 1<i<length(indices)
			idx, = Tuple(sim.lattice_domain.indices[i])
			if !(indices[i] == indices[i-1] == indices[i+1])
				for j ∈ eachindex(r) X[j] = lattice[idx+r[j]] end
		        simulation_dielectric!(E,sim,X)
				sim.ε[i] = mean(E)
			end
		else
			# placeholder for periodic lattice domains
		end
    end
    return nothing
end

function smooth_pump!(sim::Simulation{2}, num_sub_pixel::Integer = NUM_SUBPIXELS)
	indices = sim.dispersive_domain_indices
    lattice = sim.lattice
    X = Vector{Point{1}}(undef,num_sub_pixel)
    F = Vector{ComplexF64}(undef,num_sub_pixel)
	r = LinRange(-.5,.5,num_sub_pixel)
    for i ∈ eachindex(indices)
		if 1<i<length(indices)
			idx, = Tuple(sim.lattice_domain.indices[i])
			if !(indices[i] == indices[i-1] == indices[i+1])
				for j ∈ eachindex(r) X[j] = lattice[idx+r[j]] end
		        simulation_pump!(F,sim,X)
				sim.F[i] = mean(F)
			end
		end
    end
    return nothing
end

################################################################################
# extras

function Base.getproperty(sim::Simulation{2}, sym::Symbol)
	if sym==:lattice_domain
		return getfield(sim,:lattice_domains)[1]
	elseif sym==:latticesize
		indices = getfield(sim,:lattice_domains)[1].indices
		imin, imax = extrema(map(ld->ld[1],indices))
		jmin, jmax = extrema(map(ld->ld[2],indices))
		return imax-imin+1, jmax-jmin+1
	elseif Base.sym_in(sym,(:k10,:k1₀))
		return getfield(sim,:k₁₀)
	elseif Base.sym_in(sym,(:k20,:k2₀))
		return getfield(sim,:k₂₀)
	elseif Base.sym_in(sym,(:k30,:k3₀))
		return getfield(sim,:k₃₀)
	elseif Base.sym_in(sym,(:shape,:boundary,:lattice))
		return getproperty(getfield(sim,:lattice_domains)[1],sym)
	elseif Base.sym_in(sym,propertynames(getproperty(getfield(sim,:lattice_domains)[1],:lattice)))
		return getproperty(getfield(sim,:lattice_domains)[1],sym)
	elseif Base.sym_in(sym,(:lat,:lattice,:Lat,:Lattice))
		return getfield(sim,:lattice_domains)[1].lattice
	else
		return getfield(sim,sym)
	end
end

function Base.propertynames(sim::Simulation{2,Symmetric}, private=false)
	if private
		return fieldnames(Simulation)
	else
		return (:shape, :boundary, :lattice, :lattice_domain, :lattice_domains, :nondispersive_domains, :dispersive_domains, :ε, :F, :x, :χ, :ω₀, :k1₀, :k2₀, :k3₀, :lattice, propertynames(sim.lattice)...)
	end
end

################################################################################
# Pretty Printing

import ...PRINTED_COLOR_GOOD
import ...PRINTED_COLOR_WARN
import ...PRINTED_COLOR_DARK

function Base.show(io::IO,sim::Simulation{2,CLASS}) where CLASS
	if sim.smoothed[]
		printstyled(io,"Smoothed ",color=PRINTED_COLOR_GOOD)
	else
		printstyled(io,"Unsmoothed ",color=PRINTED_COLOR_WARN)
	end
	print(io,"2D ")
	CLASS<:Symmetric ? print(io,"Symmetric") : nothing
	CLASS<:Unsymmetric ? print(io,"Unsymmetric") : nothing
	printstyled(io," Simulation\n",color=PRINTED_COLOR_DARK)
    print(io,"\tn sites: ")
	printstyled(io,length(sim),"\n",color=:light_cyan)
	println(io,"\t====================")
	println(IOContext(io,:tabbed2=>true),sim.boundary)
	println(io)
	println(IOContext(io,:tabbed2=>true),sim.lattice)
	println(io)
	domains = (sim.nondispersive_domains...,sim.dispersive_domains...)
	for d ∈ eachindex(domains)
		print(io,"\t\tDomain ", d, " ")
		print(IOContext(io,:tabbed2=>true),domains[d])
		d < length(domains) ? println(io) : nothing
	end
end

end

using .Simulations2DSymmetric


text:
# Export
export slew_time_single_axis
export compute_los_vector
export constraint_agility_single_axis

##########################
# Spacecraft Slew Models #
##########################

"""
Compute the required slew time to maneuver from pointing aligned with the start
vector to the end vector.

Arguments:
- `z_start::Array{<:Real, 1}` Initial pointing vector
- `z_end::Array{<:Real, 1}` Final pointing vector
- `slew_rate::Real` Slew rate [deg/s]

Returns:
- `slew_time::Float64` Time required to slew from start to end axis
"""
function slew_time_single_axis(z_start::Array{<:Real, 1}, z_end::Array{<:Real, 1}; slew_rate=1.0::Real)

    if length(z_start) != 3
        throw(ErrorException("Length of line-of-sight vector z_start must be 3 not $(length(z_start))."))
    end

    if length(z_end) != 3
        throw(ErrorException("Length of line-of-sight vector z_end must be 3 not $(length(z_end))."))
    end

    # Use the dot product to get the angle between two vectors
    proj = dot(z_start, z_end)/(norm(z_start)*norm(z_end))

    # Normalize argument if greather than one
    if proj > 1.0
        proj = proj/norm(proj)
    end

    # Compute slew angle
    slew_angle = acos(proj)*RAD2DEG

    # Compute slew time
    slew_time = slew_angle/slew_rate

    return slew_time
end

"""
Compute line of sight vector
"""
function compute_los_vector(opp::Opportunity, epc::Epoch)
    # Compute starting and ending attitude
    r_eci_location = rECEFtoECI(epc) * opp.location.ecef
    r_eci_sat = state(opp.spacecraft.tle, epc)[1:3]

    # Compute initial look angle
    z_los = r_eci_location - r_eci_sat # Compute look angle

    # Normalize vector
    z_los = z_los/norm(z_los)

    return z_los
end

"""
Computes whether it is feasible for a spacecraft to slew from the start collect
to the end collect. The transition is based on the spacecraft ability to satisfy
single_axis pointing constraints.

Arguments:
- `start_collect::Opportunity` Initial collect. Spacecraft is assumed pointing here.
- `end_collect::Opportunity` End collect. Spacecraft is assumed pointing here.

Returns:
- `feasible::Bool` `true` if the transition from start to end is feasible. `false` otherwise
"""
function constraint_agility_single_axis(start_collect::Opportunity, end_collect::Opportunity)
    # Check inputs
    if start_collect.spacecraft.id != end_collect.spacecraft.id
        throw(ErrorException("Collect opportunities $(start_collect.id) and $(end_collect.id) for different spacecraft."))
    end

    # Get spacecraft slew-rate
    slew_rate = start_collect.spacecraft.slew_rate

    # Maximum slew time based on worst-case reorigntation from slew-rate
    max_slew_time = 180.0/slew_rate

    # Can't go backwards in time
    if end_collect.t_start < start_collect.t_end
        return false
    end

    # Exit early if time separation is large enough to guarantee feasibility
    if (end_collect.t_start - start_collect.t_end) > max_slew_time
        return true
    end

    # Compute start and end line of sight vectors
    z_start = compute_los_vector(start_collect, start_collect.t_end)
    z_end   = compute_los_vector(end_collect, end_collect.t_start)

    # Compute slew time between orientations
    t_slew = slew_time_single_axis(z_start, z_end, slew_rate=slew_rate)

    # Check if slew time is less than the time available to complete the slew
    if t_slew > (end_collect.t_start - start_collect.t_end)
        return false
    else
        return true
    end
end


text:
# This simply holds the feasibility checks
function isFeasiblePolyHedron(A, b, x)
    #=
    Checks if all polyhedron checks are satisfied
    Namely: Ax ≤ b
    =#

    state = A * x - b

    if maximum(state) > 0
        return false
    end

    return true
end

max_stars_repo_name: danineamati/TrajOptSOCPs


text:
module TestVcat
    using Base.Test
    using DataFrames

    null_df = DataFrame(Int, 0, 0)
    df = DataFrame(Int, 4, 3)

    # Assignment of rows
    df[1, :] = df[1, :]
    df[1:2, :] = df[1:2, :]

    # Broadcasting assignment of rows
    df[1, :] = 1

    # Assignment of columns
    df[1] = zeros(4)

    # Broadcasting assignment of columns
    df[:, 1] = 1
    df[1] = 3
    df[:x3] = 2

    vcat(null_df)
    vcat(null_df, null_df)
    vcat(null_df, df)
    vcat(df, null_df)
    vcat(df, df)
    vcat(df, df, df)

    alt_df = deepcopy(df)
    vcat(df, alt_df)
    df[1] = zeros(Int, nrow(df))
    # Fail on non-matching types
    vcat(df, alt_df)

    alt_df = deepcopy(df)
    names!(alt_df, [:A, :B, :C])
    # Fail on non-matching names
    vcat(df, alt_df)

    # df[:, 1] = dvzeros(Int, nrow(df))
end


text:
# ---
# title: 934. Shortest Bridge
# id: problem934
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Depth-first Search, Breadth-first Search
# link: <https://leetcode.com/problems/shortest-bridge/description/>
# hidden: true
# ---
#
# In a given 2D binary array `A`, there are two islands.  (An island is a
# 4-directionally connected group of `1`s not connected to any other 1s.)
#
# Now, we may change `0`s to `1`s so as to connect the two islands together to
# form 1 island.
#
# Return the smallest number of `0`s that must be flipped.  (It is guaranteed
# that the answer is at least 1.)
#
#
#
# **Example 1:**
#
#
#
#     Input: A = [[0,1],[1,0]]
#     Output: 1
#
#
# **Example 2:**
#
#
#
#     Input: A = [[0,1,0],[0,0,0],[0,0,1]]
#     Output: 2
#
#
# **Example 3:**
#
#
#
#     Input: A = [[1,1,1,1,1],[1,0,0,0,1],[1,0,1,0,1],[1,0,0,0,1],[1,1,1,1,1]]
#     Output: 1
#
#
#
#
# **Constraints:**
#
#   * `2 <= A.length == A[0].length <= 100`
#   * `A[i][j] == 0` or `A[i][j] == 1`
#
#
## @lc code=start
using LeetCode

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


text:
using PEG
using Test

function parse_fails_at(rule, input)
  try
    parse_whole(rule, input)
    "parse succeeded!"
  catch err
    isa(err, Meta.ParseError) || rethrow()
    m = match(r"^On line \d+, at column \d+ \(byte (\d+)\):", err.msg)
    m == nothing && rethrow()
    parse(Int, m.captures[1])
  end
end

# test miscellaneous features not covered by calculator or meta

# repeat range
@rule a34 = "a"[3:4]
@test parse_fails_at(a34, "aa") == 1
@test parse_whole(a34, "aaa") == ["a", "a", "a"]
@test parse_whole(a34, "aaaa") == ["a", "a", "a", "a"]
@test parse_fails_at(a34, "aaaaa") == 1

# 1 or more times
@rule aplus = "a"[+]
@test parse_fails_at(aplus, "") == 1
@test parse_whole(aplus, "a") == ["a"]
@test parse_whole(aplus, "aaa") == ["a", "a", "a"]

# positive lookahead
@rule abc = +"abc" & r"\w+"
@test parse_fails_at(abc, "abdargh") == 1
@test parse_whole(abc, "abc") == [nothing, "abc"]
@test parse_whole(abc, "abcfoo") == [nothing, "abcfoo"]

# negative lookahead
@rule notabc = -"abc" & r"\w+"
@test parse_whole(notabc, "ab") == [nothing, "ab"]
@test parse_whole(notabc, "abd") == [nothing, "abd"]
@test parse_whole(notabc, "abdfoo") == [nothing, "abdfoo"]
@test parse_fails_at(notabc, "abcargh") == 1

# interpolating semantics
@rule interp = r"\d+" & r"[a-z]+" & r"\d+" > (a,b,c) -> b * Base.string(parse(Int, a) + parse(Int, c))
@test parse_whole(interp, "123abc456") == "abc579"

# multibyte characters in input and rule ("∈" is 3 bytes)
@rule inop = "∈"
@rule inexpr = r"[a-z]+" & inop & r"[a-z]+"
@test parse_whole(inexpr, "foo∈bar") == ["foo", "∈", "bar"]
@test parse_fails_at(inexpr, "foo∈42") == 7
@test parse_fails_at(inexpr, "foo∈") == 7
@test parse_fails_at(inexpr, "foo42") == 4
@test parse_fails_at(inexpr, "foo") == 4


text:
pushfirst!(LOAD_PATH, joinpath(@__DIR__, "..")) # add ParameterEstimocean to environment stack

using
  Documenter,
  Literate,
  CairoMakie,  # so that Literate.jl does not capture precompilation output or warnings
  Distributions,
  ParameterEstimocean

# Gotta set this environment variable when using the GR run-time on CI machines.
# This happens when examples, e.g., use Plots.jl to make plots and movies.
# See: https://github.com/jheinen/GR.jl/issues/278
ENV["GKSwstype"] = "100"
ENV["DATADEPS_ALWAYS_ACCEPT"] = "true"

#####
##### Generate examples
#####

const EXAMPLES_DIR = joinpath(@__DIR__, "..", "examples")
const OUTPUT_DIR   = joinpath(@__DIR__, "src/literated")

to_be_literated = [
  "intro_to_observations.jl",
  "intro_to_inverse_problems.jl",
  "exploring_priors.jl",
  "perfect_convective_adjustment_calibration.jl",
  "single_case_lesbrary_ri_based_calibration.jl",
  "multi_case_lesbrary_ri_based_calibration.jl",
  "perfect_baroclinic_adjustment_calibration.jl"
]

for file in to_be_literated
    filepath = joinpath(EXAMPLES_DIR, file)
    Literate.markdown(filepath, OUTPUT_DIR; flavor = Literate.DocumenterFlavor())
end

#####
##### Build and deploy docs
#####

# Set up a timer to print a space ' ' every 240 seconds. This is to avoid CI
# timing out when building demanding Literate.jl examples.
Timer(t -> println(" "), 0, interval=240)

format = Documenter.HTML(
  collapselevel = 2,
     prettyurls = get(ENV, "CI", nothing) == "true",
      canonical = "https://clima.github.io/ParameterEstimocean/dev/",
)

pages = [
    "Home" => "index.md",
    "Installation Instructions" => "installation_instructions.md",

    "Intro to observations" => "literated/intro_to_observations.md",
    "Intro to inverse problems" => "literated/intro_to_inverse_problems.md",
    "Exploring Prior distributions" => "literated/exploring_priors.md",

    "Examples" => [
        "literated/perfect_convective_adjustment_calibration.md",
        "literated/single_case_lesbrary_ri_based_calibration.md",
        "literated/multi_case_lesbrary_ri_based_calibration.md",
        "literated/perfect_baroclinic_adjustment_calibration.md"
        ],

    "Library" => [
        "Contents"       => "library/outline.md",
        "Public"         => "library/public.md",
        "Private"        => "library/internals.md",
        "Function index" => "library/function_index.md",
        ],
]

makedocs(
   sitename = "ParameterEstimocean.jl",
    modules = [ParameterEstimocean],
     format = format,
      pages = pages,
    doctest = true,
     strict = true,
      clean = true,
  checkdocs = :exports
)

deploydocs(        repo = "github.com/CliMA/ParameterEstimocean.jl",
               versions = ["stable" => "v^", "v#.#.#", "dev" => "dev"],
              forcepush = true,
              devbranch = "main",
           push_preview = true
)


text:
#=
Directory
---------
Akl Toussaint = line 13
Graham Scan = line 89
Monotone Chain = line 110
Quickhull
=#

#=
Akl Toussaint algorithm
=#
function _akltoussaint{T}(p::Vector{Point{2, T}})

    # Get number of points
    npts = length(p)
    if npts <= 4
        return p
    end

    # Get min and max values for x and y
    xmin, xminind, xmax, xmaxind = extremax(p)
    ymin, yminind, ymax, ymaxind = extremay(p)

    # Find unique set of points
    uptsind = unique([xminind, ymaxind, xmaxind, yminind])
    upts = p[uptsind]
    nu = length(uptsind)

    if nu == 4
        out = _akltoussaint4(p, upts, npts)
    elseif nu == 3
        out =  _akltoussaint3(p, upts, npts)
    else
        warn("Could not apply Akl Toussaint. Returning all points")
        out = p
    end

    return out
end

function _akltoussaint4{T<:Real}(p::Vector{Point{2, T}}, upts::Vector{Point{2, T}}, npts::Int)

    # Get the points from upts
    pt1, pt2, pt3, pt4 = upts

    # Allocate space for output of each triangle
    out = Array(Bool, npts)

    # Create a "T inverse" matrix for each triangle
    tinv_1 = create_t_inv(pt1, pt2, pt3)
    tinv_2 = create_t_inv(pt1, pt4, pt3)

    @inbounds for n=1:npts
        temp = (p[n] - pt3)

        λ1_1, λ2_1 = _unsafe_mult_2by2_point(tinv_1, temp)
        λ1_2, λ2_2 = _unsafe_mult_2by2_point(tinv_2, temp)

        out[n] = check_λ(λ1_1, λ2_1) * check_λ(λ1_2, λ2_2)
    end

    return p[out]
end

function _akltoussaint3{T<:Real}(p::Vector{Point{2, T}}, upts::Vector{Point{2, T}}, npts::Int)

    # Get the points from upts
    pt1, pt2, pt3 = upts

    # Allocate space for output of triangle
    out = Array(Bool, npts)

    # Create a "T inverse" matrix for each triangle
    tinv = create_t_inv(pt1, pt2, pt3)

    @inbounds for n=1:npts
        temp = (p[n] - pt3)

        λ1, λ2 = _unsafe_mult_2by2_point(tinv, temp)

        out[n] = check_λ(λ1, λ2)
    end

    return p[out]
end

#
# Graham Scan Algorithm
#
function _grahamscan{T<:Real}(p::Vector{Point{2, T}})

    # Get yminind
    xminind = indmin(p)
    xmin = p[xminind]

    # Create function to sort by angles
    lt(a, b) = xmin == a ? true : xmin == b ? false : ccw(xmin, a, b)
    psort = sort!(p, lt=lt)

    # Add the starting point at end so we go full circle
    push!(psort, xmin)

    # Call function that works around points
    ep = wrappoints(psort)

    return ep[1:end-1]
end

#
# Monotone Chain
#
function _monotonechain{T<:Real}(p::Vector{Point{2, T}})

    # Sort points
    psort = sort!(p)

    # Split into upper and lower
    xl, xu = split_xL_xU(psort)

    lh = wrappoints(xl)
    uh = wrappoints(reverse!(xu))

    return [lh[1:end-1]; uh[1:end-1]]
end

# function _quickhull{T<:Real}(p::Vector{Point{2, T}})

#     # Note: I suspect this algorithm can be implemented recursively
#     #       by calling it on two points, separating them into two
#     #       sets and then calling it on the endpoints of those
#     #       sets. Think carefully about how to implement this the
#     #       right way.
#     # Note: In a loop, at same time, check whether a point lies within the
#     #       the set of points that should be removed and whether it is the
#     #       furthest point from line.
# end


text:
module ParallelTreeSampling

using Distributions
using ImportanceWeightedRiskMetrics
using POMDPModelTools
using POMDPPolicies
using POMDPs
using ProgressMeter
using Random
using StatsBase

include("common_utils.jl")

include("tree.jl")

export PISSolver, PISPlanner
include("solver.jl")

include("rollout.jl")
include("search.jl")

export GenericDiscreteNonParametric
include("distributions.jl")

include("visualization.jl")

export plot_estimates, plot_histogram
include("plots.jl")

end


text:
# These are the tests for deprecated features, they should be deleted along with them

@testset "Trie: path iterator" begin
    t = Trie{Int}()
    t["rob"] = 27
    t["roger"] = 52
    t["kevin"] = Int8(11)
    t0 = t
    t1 = t0.children['r']
    t2 = t1.children['o']
    t3 = t2.children['b']
    @test collect(path(t, "b")) == [t0]
    @test collect(path(t, "rob")) == [t0, t1, t2, t3]
    @test collect(path(t, "robb")) == [t0, t1, t2, t3]
    @test collect(path(t, "ro")) == [t0, t1, t2]
    @test collect(path(t, "roa")) == [t0, t1, t2]
end

@testset "top" begin
    hh = BinaryMinHeap{Float64}([1,2,3])
    @test top(hh) == 1
end

function test_reverse_iter(it::T) where T
    arr = [i for i in it]
    index = length(arr)
    for i in reverse_iter(it)
        @test arr[index] == i
        index -= 1
    end

    @test reverse(arr) == [i for i in reverse_iter(it)]
end
@testset "reverse_iter" begin
    @testset "Queue" begin
        q = Queue{Int}(); push!(q, 1); push!(q, 2)
        test_reverse_iter(q)
    end
    @testset "Stack" begin
        s = Stack{Int}(); push!(s, 1); push!(s, 2)
        test_reverse_iter(s)
    end
end

@testset "enqueue! dequeue!" begin
    s = Queue{Int}(5)
    n = 100

    @test length(s) == 0
    @test eltype(s) == Int
    @test eltype(typeof(s)) == Int
    @test isempty(s)
    @test_throws ArgumentError first(s)
    @test_throws ArgumentError last(s)
    @test_throws ArgumentError dequeue!(s)

    for i = 1 : n
        enqueue!(s, i)
        @test first(s) == 1
        @test last(s) == i
        @test !isempty(s)
        @test length(s) == i
    end

    for i = 1 : n
        x = dequeue!(s)
        @test x == i
        if i < n
            @test first(s) == i + 1
            @test last(s) == n
        else
            @test_throws ArgumentError first(s)
            @test_throws ArgumentError last(s)
        end
        @test isempty(s) == (i == n)
        @test length(s) == n - i
    end
end


text:
module TestModule
    using GeometricAlgebra

    @basis "++-" prefix=g export_pseudoscalar=false modname=GA
end

@testset "Basis macro" begin
    @test v123 == Blade(1, 8)

    @test isdefined(TestModule, :g)
    @test isdefined(TestModule, :g1)
    @test isdefined(TestModule, :g123)
    @test isdefined(TestModule, :g13)
    @test !isdefined(TestModule, :I)
    @test isdefined(TestModule.GA, :I)
    @test :I ∉ names(TestModule.GA)
    @test TestModule.g123 == v123
    @test @isdefined(I)

    @test_throws ArgumentError("Macro options must be `(name=value)` pairs") @basis "+++" unknown_option
    @test_throws ArgumentError("Unknown option 'unknown_option'. Available options are: prefix, export_symbols, export_metadata, export_pseudoscalar, modname") @basis "+++" unknown_option=nothing
end


text:
"""
Farquhar–von Caemmerer–Berry (FvCB) model for C3 photosynthesis (Farquhar et al., 1980;
von Caemmerer and Farquhar, 1981).

Iterative implementation, i.e. the assimilation is computed iteratively over Cᵢ.

For more details on arguments, see [`Fvcb`](@ref).
This structure has several more parameters:

- `iter_A_max::Int`: maximum number of iterations allowed for the iteration on the assimilation.
- `ΔT_A::T = 1`: threshold bellow which the assimilation is considered constant. Given in
percent of change, *i.e.* 1% means that two successive assimilations with less than 1%
difference in value are considered the same value.
"""
struct FvcbIter{T} <: AbstractAModel
    Tᵣ::T
    VcMaxRef::T
    JMaxRef::T
    RdRef::T
    Eₐᵣ::T
    O₂::T
    Eₐⱼ::T
    Hdⱼ::T
    Δₛⱼ::T
    Eₐᵥ::T
    Hdᵥ::T
    Δₛᵥ::T
    α::T
    θ::T
    iter_A_max::Int
    ΔT_A::T
end

function FvcbIter(;Tᵣ = 25.0, VcMaxRef = 200.0, JMaxRef = 250.0, RdRef = 0.6, Eₐᵣ = 46390.0,
    O₂= 210.0, Eₐⱼ = 29680.0, Hdⱼ = 200000.0, Δₛⱼ = 631.88, Eₐᵥ = 58550.0, Hdᵥ = 200000.0,
    Δₛᵥ = 629.26, α = 0.425, θ = 0.90, iter_A_max = 20, ΔT_A = 1.0)

    # Add type promotion in case we want to use e.g. Measurements for one parameter only and
    # we don't want to set each parameter to ± 0 by hand.
    param_float = promote(Tᵣ, VcMaxRef, JMaxRef, RdRef, Eₐᵣ, O₂, Eₐⱼ, Hdⱼ, Δₛⱼ, Eₐᵥ, Hdᵥ, Δₛᵥ, α, θ, ΔT_A)
    FvcbIter(
        param_float[1:end-1]...,
        iter_A_max,
        param_float[end]
    )
end

function inputs(::FvcbIter)
    (:PPFD,:Tₗ,:Gbc)
end

function outputs(::FvcbIter)
    (:A,:Gₛ,:Cᵢ,:Cₛ)
end

Base.eltype(x::FvcbIter) = typeof(x).parameters[1]

"""
    assimilation!(leaf::LeafModels{I,E,<:FvcbIter,<:AbstractGsModel,S}, meteo, constants = Constants())

Photosynthesis using the Farquhar–von Caemmerer–Berry (FvCB) model for C3 photosynthesis
 (Farquhar et al., 1980; von Caemmerer and Farquhar, 1981).
Computation is made following Farquhar & Wong (1984), Leuning et al. (1995), and the
Archimed model.

Iterative implementation, i.e. the assimilation is computed iteratively over Cᵢ. For the
analytical resolution, see [`Fvcb`](@ref).

# Returns

Modify the first argument in place for A, Gₛ and Cᵢ:

- A: carbon assimilation (μmol[CO₂] m-2 s-1)
- Gₛ: stomatal conductance for CO₂ (mol[CO₂] m-2 s-1)
- Cᵢ: intercellular CO₂ concentration (ppm)

# Arguments

- `leaf::LeafModels{.,.,<:FvcbIter,<:AbstractGsModel,.}`: A [`LeafModels`](@ref) struct holding the parameters for
the model with initialisations for:
    - `Tₗ` (°C): leaf temperature
    - `PPFD` (μmol m-2 s-1): absorbed Photosynthetic Photon Flux Density
    - `Gbc` (mol m-2 s-1): boundary conductance for CO₂
    - `Dₗ` (mol m-2 s-1): vapour pressure difference between the surface and the air saturation
    vapour pressure in case you're using the stomatal conductance model of [`Medlyn`](@ref).
- `meteo`: meteorology structure, see [`Atmosphere`](@ref)
- `constants = Constants()`: physical constants. See [`Constants`](@ref) for more details

# Note

`Tₗ`, `PPFD`, `Gbc` (and `Dₗ` if you use [`Medlyn`](@ref)) must be initialised by providing
them as keyword arguments (see examples). If in doubt, it is simpler to compute the energy
balance of the leaf with the photosynthesis to get those variables. See
[`energy_balance`](@ref) for more details.

# Examples

```julia
meteo = Atmosphere(T = 20.0, Wind = 1.0, P = 101.3, Rh = 0.65)

leaf = LeafModels(photosynthesis = FvcbIter(),
            stomatal_conductance = Medlyn(0.03, 12.0),
            Tₗ = 25.0, PPFD = 1000.0, Gbc = 0.67, Dₗ = meteo.VPD)
# NB: we need  to initalise Tₗ, PPFD and Gbc.

assimilation!(leaf,meteo,Constants())
leaf.status.A
leaf.status.Cᵢ
```

# References

Baldocchi, Dennis. 1994. « An analytical solution for coupled leaf photosynthesis and
stomatal conductance models ». Tree Physiology 14 (7-8‑9): 1069‑79.
https://doi.org/10.1093/treephys/14.7-8-9.1069.

Farquhar, G. D., S. von von Caemmerer, et J. A. Berry. 1980. « A biochemical model of
photosynthetic CO2 assimilation in leaves of C3 species ». Planta 149 (1): 78‑90.

Leuning, R., F. M. Kelliher, DGG de Pury, et E.D. Schulze. 1995. Leaf nitrogen,
photosynthesis, conductance and transpiration: scaling from leaves to canopies ». Plant,
Cell & Environment 18 (10): 1183‑1200.
"""
function assimilation!(leaf::LeafModels{I,E,<:FvcbIter,<:AbstractGsModel,S}, meteo,
    constants = Constants()) where {I,E,S}

    # Start with a probable value for Cₛ and Cᵢ:
    leaf.status.Cₛ = meteo.Cₐ
    leaf.status.Cᵢ = leaf.status.Cₛ * 0.75

    # Tranform Celsius temperatures in Kelvin:
    Tₖ = leaf.status.Tₗ - constants.K₀
    Tᵣₖ = leaf.photosynthesis.Tᵣ - constants.K₀

    # Temperature dependence of the parameters:
    Γˢ = Γ_star(Tₖ,Tᵣₖ,constants.R) # Gamma star (CO2 compensation point) in μmol mol-1
    Km = get_km(Tₖ,Tᵣₖ,leaf.photosynthesis.O₂,constants.R) # effective Michaelis–Menten coefficient for CO2

    # Maximum electron transport rate at the given leaf temperature (μmol m-2 s-1):
    JMax = arrhenius(leaf.photosynthesis.JMaxRef,leaf.photosynthesis.Eₐⱼ,Tₖ,Tᵣₖ,
                        leaf.photosynthesis.Hdⱼ,leaf.photosynthesis.Δₛⱼ,constants.R)
    # Maximum rate of Rubisco activity at the given leaf temperature (μmol m-2 s-1):
    VcMax = arrhenius(leaf.photosynthesis.VcMaxRef,leaf.photosynthesis.Eₐᵥ,Tₖ,Tᵣₖ,
                        leaf.photosynthesis.Hdᵥ,leaf.photosynthesis.Δₛᵥ,constants.R)
    # Rate of mitochondrial respiration at the given leaf temperature (μmol m-2 s-1):
    Rd = arrhenius(leaf.photosynthesis.RdRef,leaf.photosynthesis.Eₐᵣ,Tₖ,Tᵣₖ,constants.R)
    # Rd is also described as the CO2 release in the light by processes other than the PCO
    # cycle, and termed "day" respiration, or "light respiration" (Harley et al., 1986).

    # Actual electron transport rate (considering intercepted PAR and leaf temperature):
    J = get_J(leaf.status.PPFD, JMax, leaf.photosynthesis.α, leaf.photosynthesis.θ) # in μmol m-2 s-1
    # RuBP regeneration
    Vⱼ = J / 4

    # First iteration to initialise the values for A and Gₛ:
    # Net assimilation (μmol m-2 s-1)
    leaf.status.A = Fvcb_net_assimiliation(leaf.status.Cᵢ,Vⱼ,Γˢ,VcMax,Km,Rd)

    iter = true
    iter_inc = 1

    while iter
        # Stomatal conductance (mol[CO₂] m-2 s-1)
        leaf.status.Gₛ = gs(leaf,meteo)
        # Surface CO₂ concentration (ppm):
        leaf.status.Cₛ = min(meteo.Cₐ, meteo.Cₐ - leaf.status.A / leaf.status.Gbc)
        # Intercellular CO₂ concentration (ppm):
        leaf.status.Cᵢ = min(leaf.status.Cₛ, leaf.status.Cₛ - leaf.status.A / leaf.status.Gₛ)

        if leaf.status.Cᵢ <= zero(leaf.status.Cᵢ)
            leaf.status.Cᵢ = 1e-9
            A_new = -Rd
        else
            # Net assimilation (μmol m-2 s-1):
            A_new = Fvcb_net_assimiliation(leaf.status.Cᵢ,Vⱼ,Γˢ,VcMax,Km,Rd)
        end

        if abs(A_new - leaf.status.A) / leaf.status.A <= leaf.photosynthesis.ΔT_A ||
            iter_inc == leaf.photosynthesis.iter_A_max

            iter = false
        end

        leaf.status.A = A_new

        iter_inc += 1
    end
end

"""
    Fvcb_net_assimiliation(Cᵢ,Vⱼ,Γˢ,VcMax,Km,Rd)

Net assimilation following the Farquhar–von Caemmerer–Berry (FvCB) model for C3 photosynthesis
(Farquhar et al., 1980; von Caemmerer and Farquhar, 1981)
"""
function Fvcb_net_assimiliation(Cᵢ,Vⱼ,Γˢ,VcMax,Km,Rd)
    # Electron-transport-limited rate of CO₂ assimilation (RuBP regeneration-limited):
    Wⱼ = Vⱼ * (Cᵢ - Γˢ) / (Cᵢ + 2.0 * Γˢ)
    # See Von Caemmerer, Susanna. 2000. Biochemical models of leaf photosynthesis.
    # Csiro publishing, eq. 2.23.
    # NB: here the equation is modified because we use Vⱼ instead of J, but it is the same.

    # Rubisco-carboxylation-limited rate of CO₂ assimilation (RuBP activity-limited):
    Wᵥ = VcMax * (Cᵢ - Γˢ) / (Cᵢ + Km)

    # Net assimilation (μmol m-2 s-1):
    A = min(Wᵥ,Wⱼ) - Rd
    return A
end


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

mutable struct ItemSimilarityRecommender_ConnectedItem <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function ItemSimilarityRecommender_ConnectedItem(; kwargs...)
        obj = new(meta(ItemSimilarityRecommender_ConnectedItem), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct ItemSimilarityRecommender_ConnectedItem
const __meta_ItemSimilarityRecommender_ConnectedItem = Ref{ProtoMeta}()
function meta(::Type{ItemSimilarityRecommender_ConnectedItem})
    ProtoBuf.metalock() do
        if !isassigned(__meta_ItemSimilarityRecommender_ConnectedItem)
            __meta_ItemSimilarityRecommender_ConnectedItem[] = target = ProtoMeta(ItemSimilarityRecommender_ConnectedItem)
            allflds = Pair{Symbol,Union{Type,String}}[:itemId => UInt64, :similarityScore => Float64]
            meta(target, ItemSimilarityRecommender_ConnectedItem, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_ItemSimilarityRecommender_ConnectedItem[]
    end
end
function Base.getproperty(obj::ItemSimilarityRecommender_ConnectedItem, name::Symbol)
    if name === :itemId
        return (obj.__protobuf_jl_internal_values[name])::UInt64
    elseif name === :similarityScore
        return (obj.__protobuf_jl_internal_values[name])::Float64
    else
        getfield(obj, name)
    end
end

mutable struct ItemSimilarityRecommender_SimilarItems <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function ItemSimilarityRecommender_SimilarItems(; kwargs...)
        obj = new(meta(ItemSimilarityRecommender_SimilarItems), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct ItemSimilarityRecommender_SimilarItems
const __meta_ItemSimilarityRecommender_SimilarItems = Ref{ProtoMeta}()
function meta(::Type{ItemSimilarityRecommender_SimilarItems})
    ProtoBuf.metalock() do
        if !isassigned(__meta_ItemSimilarityRecommender_SimilarItems)
            __meta_ItemSimilarityRecommender_SimilarItems[] = target = ProtoMeta(ItemSimilarityRecommender_SimilarItems)
            allflds = Pair{Symbol,Union{Type,String}}[:itemId => UInt64, :similarItemList => Base.Vector{ItemSimilarityRecommender_ConnectedItem}, :itemScoreAdjustment => Float64]
            meta(target, ItemSimilarityRecommender_SimilarItems, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_ItemSimilarityRecommender_SimilarItems[]
    end
end
function Base.getproperty(obj::ItemSimilarityRecommender_SimilarItems, name::Symbol)
    if name === :itemId
        return (obj.__protobuf_jl_internal_values[name])::UInt64
    elseif name === :similarItemList
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{ItemSimilarityRecommender_ConnectedItem}
    elseif name === :itemScoreAdjustment
        return (obj.__protobuf_jl_internal_values[name])::Float64
    else
        getfield(obj, name)
    end
end

mutable struct ItemSimilarityRecommender <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function ItemSimilarityRecommender(; kwargs...)
        obj = new(meta(ItemSimilarityRecommender), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct ItemSimilarityRecommender
const __meta_ItemSimilarityRecommender = Ref{ProtoMeta}()
function meta(::Type{ItemSimilarityRecommender})
    ProtoBuf.metalock() do
        if !isassigned(__meta_ItemSimilarityRecommender)
            __meta_ItemSimilarityRecommender[] = target = ProtoMeta(ItemSimilarityRecommender)
            fnum = Int[1,2,3,10,11,12,13,20,21]
            allflds = Pair{Symbol,Union{Type,String}}[:itemItemSimilarities => Base.Vector{ItemSimilarityRecommender_SimilarItems}, :itemStringIds => StringVector, :itemInt64Ids => Int64Vector, :itemInputFeatureName => AbstractString, :numRecommendationsInputFeatureName => AbstractString, :itemRestrictionInputFeatureName => AbstractString, :itemExclusionInputFeatureName => AbstractString, :recommendedItemListOutputFeatureName => AbstractString, :recommendedItemScoreOutputFeatureName => AbstractString]
            meta(target, ItemSimilarityRecommender, allflds, ProtoBuf.DEF_REQ, fnum, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_ItemSimilarityRecommender[]
    end
end
function Base.getproperty(obj::ItemSimilarityRecommender, name::Symbol)
    if name === :itemItemSimilarities
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{ItemSimilarityRecommender_SimilarItems}
    elseif name === :itemStringIds
        return (obj.__protobuf_jl_internal_values[name])::StringVector
    elseif name === :itemInt64Ids
        return (obj.__protobuf_jl_internal_values[name])::Int64Vector
    elseif name === :itemInputFeatureName
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :numRecommendationsInputFeatureName
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :itemRestrictionInputFeatureName
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :itemExclusionInputFeatureName
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :recommendedItemListOutputFeatureName
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :recommendedItemScoreOutputFeatureName
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    else
        getfield(obj, name)
    end
end

export ItemSimilarityRecommender_ConnectedItem, ItemSimilarityRecommender_SimilarItems, ItemSimilarityRecommender


text:
# DTU02435 Decision Making under Uncertainty
# Assignment 3, Task 2
# Solve by heuristics and Benders Algorithm / Gurobi
# Edward J. Xu
# April 27th, 2019
push!(LOAD_PATH, "$(homedir())/Desktop/optim")
cd("$(homedir())/Desktop/optim")
using JuMP
using GLPKMathProgInterface
using Gurobi
using CSV
using DataFrames
using BlackBoxOptim

function processData()
    datf_numMax = CSV.read("MaxPlanes.csv", header = ["p", "num"], delim = ',')
    datf_cost = CSV.read("OperationCosts2.csv", header = true, delim = ',')
    datf_cap = CSV.read("PlaneCapacity.csv", header = ["p", "cap"], delim = ',')
    datf_route = CSV.read("DemandsRoute.csv", header = true, delim = ',')
    datf_costAdmin = CSV.read("CostAdquisition.csv", header = ["p", "cost"], delim = ',')
    datf_prob = CSV.read("Probability.csv", header = ["prob"], delim = ',')
    #
    datf_demand = CSV.read("Scenario.csv", header = true, delim = ',')
    mat4_demand = zeros(4, 5, 12, 10)
    for i = 1: 4
        for j = 1: 5
            for w = 1: 10
                mat4_demand[i, j, :, w] = datf_demand[w] * datf_route[i, j+1]
            end
        end
    end
    return (datf_numMax.num, datf_cost, datf_cap.cap, datf_costAdmin.cost, mat4_demand, datf_prob.prob)
end

function solve_sub(vec_yBar, datf_cost, mat_demand, vec_cap)
    model_sub = Model(solver = GurobiSolver(Presolve = 0, OutputFlag = 0, gurobi_env))
    # model = Model(solver = GLPKSolverMIP())
    @variable(model_sub, mat3_z[1: 4, 1: 5, 1: 10] >= 0, Int)
    @variable(model_sub, mat_m[1: 4, 1: 5] >= 0)
    @objective(model_sub, Min, sum([datf_cost.cost[.&((datf_cost[:i] .== i), (datf_cost[:j] .== j),
        (datf_cost[:p] .== p))][1] * mat3_z[i, j, p] + mat_m[i, j] * 1000] for p = 1: 10, i = 1: 4,
        j = 1: 5)[1])
    @constraint(model_sub, [i = 1: 4, j = 1: 5], mat_demand[i, j] -
        sum([mat3_z[i, j, p] * vec_cap[p]] for p = 1: 10)[1] <= mat_m[i, j])
    @constraint(model_sub, [p = 1: 10], sum([mat3_z[i, j, p]] for i = 1: 4, j = 1: 5)[1] <=
        vec_yBar[p])
    solve(model_sub)
    obj = getobjectivevalue(model_sub)
    vec_result_z = getvalue(mat3_z)
    vec_result_m = getvalue(mat_m)
    return (obj, vec_result_z, vec_result_m)
end

function solve_sub_all(vec_yBarRaw; datf_cost = datf_cost, vec_cap = vec_cap, vec_costAdmin = vec_costAdmin,
        mat4_demand = mat4_demand, vec_prob = vec_prob)
# function solve_sub_all(vec_yBarRaw)
    vec_yBar = round.(vec_yBarRaw)
    println("--------------------------------------------------------------------------------\n",
            "------------------------------ Solve Sub-Problem -------------------------------\n",
            "Input: vec_yBarRaw = $(vec_yBarRaw)\n",
            "Input: vec_yBar = $(vec_yBar)\n",
            "--------------------------------------------------------------------------------\n")
    mat_obj = zeros(10, 12)
    for w = 1: 10
        for t = 1: 12
            println("Input: w = $(w), t = $(t), step = $(12 * (w - 1) + t) / 120 ")
            (mat_obj[w, t], mat_result_z, mat_result_m) = solve_sub(
                vec_yBar, datf_cost, mat4_demand[:, :, t, w], vec_cap)
        end
    end
    obj = sum([vec_yBar[p] * vec_costAdmin[p] + vec_prob[w] * mat_obj[w, t]] for p = 1: 10, w = 1: 5, t = 1: 12)[1]
    println("--------------------------------------------------------------------------------\n",
            "Result: obj = $(obj)\n",
            "--------------------------------------------------------------------------------\n")
    return obj
end

function main()
    # Begin Optimization
    println("################################################################################\n",
            "############################ 1/2. Begin Optimization ###########################\n",
            "################################################################################\n")
    mod_master = bboptimize(solve_sub_all; SearchRange = [(0.0, 12.0), (0.0, 12.0), (0.0, 12.0), (0.0, 12.0),
        (0.0, 12.0), (0.0, 12.0), (0.0, 12.0), (0.0, 8.0), (0.0, 8.0), (0.0, 8.0)], NumDimensions = 10,
        MaxTime = false)
    println("################################################################################\n",
            "############################## 2/2. Nominal Ending #############################\n",
            "################################################################################\n")
    # Print the Result
    # println("Result: obj = $(obj).")
    # println("Result: vec_y = $(vec_result_y).")
    # println("Result: mat5_z = $(vec_result_z).")
    # println("Result: mat4_m = $(vec_result_m).")
    # return (vec_result_z, vec_result_m)
end

const gurobi_env = Gurobi.Env()
const (vec_numMax, datf_cost, vec_cap, vec_costAdmin, mat4_demand, vec_prob) = processData()
main()


text:
"""
Generic template for evaluating residuals.
This assumes that projections have been handled externally.
The following flags control how residuals are evaluated:

+ `need_main`: If `true`, evaluates regression residuals.
+ `need_dist`: If `true`, evaluates difference between parameters and their projection.

**Note**: The values for each flag should be known at compile-time!
"""
function __evaluate_residuals__(problem, ϵ, extras, need_main::Bool, need_dist::Bool, need_z::Bool)
    @unpack Y, coeff, proj, res = problem
    @unpack Z = extras
    T = floattype(problem)
    X = get_design_matrix(problem)

    if need_main
        # main residuals, √(1/n) * (Y - X*B)
        a = 1 / sqrt(size(Y, 1))
        mul!(res.main.all, X, coeff.all)
        need_z && copyto!(Z, res.main.all)
        axpby!(a, Y, -a, res.main.all)

        # weighted residuals, W^1/2 * (Y - X*B)
        @inbounds for i in axes(Y, 1)
            yᵢ = view(Y, i, :)
            rᵢ = res.main.sample[i]
            wrᵢ = res.weighted.sample[i]
            normrᵢ = norm(rᵢ) / a
            wᵢ = ifelse(normrᵢ ≤ ϵ, zero(T), (normrᵢ-ϵ)/normrᵢ)
            @. wrᵢ = wᵢ * rᵢ
            if need_z
                # zᵢ = X*B if norm(rᵢ) ≤ ϵ
                # zᵢ = wᵢ*yᵢ + (1-wᵢ)*X*B otherwise
                zᵢ = view(Z, i, :)
                axpby!(wᵢ, yᵢ, 1-wᵢ, zᵢ)
            end
        end
    end

    if need_dist
        # res_dist = P(B) - B
        copyto!(res.dist.all, proj.all)
        axpby!(-one(T), coeff.all, one(T), res.dist.all)
    end

    return nothing
end

"""
Evaluate the gradiant of the regression problem. Assumes residuals have been evaluated.
"""
function __evaluate_gradient__(problem, ρ, extras)
    @unpack res, grad = problem
    n, _, _ = probdims(problem)
    X = get_design_matrix(problem)

    # ∇g_ρ(B ∣ Bₘ)ⱼ = -[aXᵀ bⱼI] * Rₘ
    a, b = 1/sqrt(n), ρ
    mul!(grad.all, X', res.weighted.all)
    axpby!(-b, res.dist.all, -a, grad.all)

    return nothing
end

function __evaluate_reg_gradient__(problem, λ, extras)
    @unpack res, grad = problem
    n, _, _ = probdims(problem)
    X = get_design_matrix(problem)

    # ∇q_λ(B ∣ Bₘ) = -aXᵀRₘ + λB
    a, b = 1/sqrt(n), λ
    mul!(grad.all, X', res.weighted.all)
    axpby!(b, problem.coeff.all, -a, grad.all)

    return nothing
end

"""
Evaluate the penalized least squares criterion. Also updates the gradient.
This assumes that projections have been handled externally.
"""
function __evaluate_objective__(problem, ϵ, ρ, extras)
    @unpack res, grad = problem

    __evaluate_residuals__(problem, ϵ, extras, true, true, false)
    __evaluate_gradient__(problem, ρ, extras)

    R, D, G = res.weighted.all, res.dist.all, grad.all
    loss = dot(R, R) # 1/n * ∑ᵢ (Zᵢ - Bᵀxᵢ)²
    dist = dot(D, D) # ∑ⱼ (P(B)ⱼ - Bⱼ)²
    obj = 1//2 * (loss + ρ * dist)
    gradsq = dot(G, G)

    return IterationResult(loss, obj, dist, gradsq)
end

function __evaluate_reg_objective__(problem, ϵ, λ, extras)
    @unpack coeff, res, grad = problem

    __evaluate_residuals__(problem, ϵ, extras, true, false, false)
    __evaluate_reg_gradient__(problem, λ, extras)

    R, B, G = res.weighted.all, coeff.all, grad.all
    loss = dot(R, R) # 1/n * ∑ᵢ (Zᵢ - Bᵀxᵢ)²
    objective = 1//2 * (loss + λ * dot(B, B))
    gradsq = dot(G, G)

    return IterationResult(loss, objective, 0.0, gradsq)
end

"""
Apply acceleration to the current iterate `x` based on the previous iterate `y`
according to Nesterov's method with parameter `r=3` (default).
"""
function __apply_nesterov__!(x, y, iter::Integer, needs_reset::Bool, r::Int=3)
    if needs_reset # Reset acceleration scheme
        copyto!(y, x)
        iter = 1
    else # Nesterov acceleration
        γ = (iter - 1) / (iter + r - 1)
        @inbounds for i in eachindex(x)
            xi, yi = x[i], y[i]
            zi = xi + γ * (xi - yi)
            x[i], y[i] = zi, xi
        end
        iter += 1
    end

    return iter
end

"""
Map a sparsity level `s` to an integer `k`, assuming `n` elements.
"""
sparsity_to_k(problem::MVDAProblem, s) = __sparsity_to_k__(problem.kernel, problem, s)
__sparsity_to_k__(::Nothing, problem::MVDAProblem, s) = round(Int, (1-s) * problem.p)
__sparsity_to_k__(::Kernel, problem::MVDAProblem, s) = round(Int, (1-s)*problem.n)

get_projection_indices(problem::MVDAProblem) = __get_projection_indices__(problem.kernel, problem)
__get_projection_indices__(::Nothing, problem::MVDAProblem) = 1:problem.p
__get_projection_indices__(::Kernel, problem::MVDAProblem) = 1:problem.n

"""
Apply a projection to model coefficients.
"""
function apply_projection(projection, problem, k)
    idx = get_projection_indices(problem)
    @unpack coeff, proj = problem
    @unpack scores = projection
    copyto!(proj.all, coeff.all)

    # projection step, might not be unique
    if problem.intercept
        projection(view(proj.all, idx, :), k)
    else
        projection(proj.all, k)
    end

    return proj.all
end

"""
Define a geometric progression recursively by the rule
```
    rho_new = min(rho_max, rho * multiplier)
```
The result is guaranteed to have type `typeof(rho)`.
"""
function geometric_progression(rho, iter, rho_max, multiplier::Real=1.2)
    return convert(typeof(rho), min(rho_max, rho * multiplier))
end

"""
Default error message for missing methods. For internal use only.
"""
not_implemented(alg, msg) = error(string(msg, " not implemented for ", alg, "."))

"""
Placeholder for callbacks in main functions.
"""
__do_nothing_callback__(iter, problem, lambda, rho, k, history) = nothing
__do_nothing_callback__(fold, problem, train_problem, data, lambda, sparsity, model_size, result) = nothing

__svd_wrapper__(A::StridedMatrix) = svd(A)
__svd_wrapper__(A::AbstractMatrix) = svd!(copy(A))

function prediction_error(model::MVDAProblem, train_set, validation_set, test_set)
    # Extract number of features to make predictions consistent.
    @unpack p = model

    # Extract data for each set.
    Tr_Y, Tr_X = train_set
    V_Y, V_X = validation_set
    T_Y, T_X = test_set

    # Helper function to make predictions on each subset and evaluate errors.
    _error = function(model, Y, X)
        # Sanity check: Y and X have the same number of rows.
        size(Y, 1) != size(X, 1) && error("Vertex assignments Y $(size(Y)) not compatible with data X $(size(X)).")
        # Predict response in vertex space; may use @batch.
        Yhat = predict(model, X)

        # Sweep through predictions and tally the mistakes.
        nincorrect = 0
        @views for i in axes(Yhat, 1)
            @inbounds yᵢ, yhatᵢ = Y[i,:], Yhat[i,:]
            true_label = model.vertex2label[yᵢ]
            predicted_label = __classify__(model, yhatᵢ)
            nincorrect += predicted_label != true_label
        end

        return 100 * nincorrect / size(Y, 1)
    end

    Tr = _error(model, Tr_Y, view(Tr_X, :, 1:p))
    V = _error(model, V_Y, view(V_X, :, 1:p))
    T = _error(model, T_Y, view(T_X, :, 1:p))

    return (Tr, V, T)
end

struct IterationResult
    loss::Float64
    objective::Float64
    distance::Float64
    gradient::Float64
end

# destructuring
Base.iterate(r::IterationResult) = (r.loss, Val(:objective))
Base.iterate(r::IterationResult, ::Val{:objective}) = (r.objective, Val(:distance))
Base.iterate(r::IterationResult, ::Val{:distance}) = (r.distance, Val(:gradient))
Base.iterate(r::IterationResult, ::Val{:gradient}) = (r.gradient, Val(:done))
Base.iterate(r::IterationResult, ::Val{:done}) = nothing

struct SubproblemResult
    iters::Int
    loss::Float64
    objective::Float64
    distance::Float64
    gradient::Float64
end

function SubproblemResult(iters, r::IterationResult)
    return SubproblemResult(iters, r.loss, r.objective, r.distance, r.gradient)
end

# destructuring
Base.iterate(r::SubproblemResult) = (r.iters, Val(:loss))
Base.iterate(r::SubproblemResult, ::Val{:loss}) = (r.loss, Val(:objective))
Base.iterate(r::SubproblemResult, ::Val{:objective}) = (r.objective, Val(:distance))
Base.iterate(r::SubproblemResult, ::Val{:distance}) = (r.distance, Val(:gradient))
Base.iterate(r::SubproblemResult, ::Val{:gradient}) = (r.gradient, Val(:done))
Base.iterate(r::SubproblemResult, ::Val{:done}) = nothing

function save_cv_results(filename, problem_info, replicates, grids)
    ϵ_grid, s_grid = grids
    example, subsets, nfeatures, nclasses, sparse2dense, kernel = problem_info
    ntrain, nvalidate, ntest = subsets
    delim = ','
    header = ["example" "ntrain" "nvalidate" "ntest" "nfeatures" "nclasses" "sparse2dense" "kernel" "replicate" "fold" "epsilon" "sparsity" "time" "train" "validation" "test"]
    open(filename, "w+") do io
        writedlm(io, header, delim)
        for (r, replicate) in enumerate(replicates)
            for k in eachindex(replicate[1]), j in eachindex(ϵ_grid), i in eachindex(s_grid)
                ϵ, s = ϵ_grid[j], 100*s_grid[i]
                Tr, V, T, t = [arr[k][i,j] for arr in replicate]
                writedlm(io, Any[example ntrain nvalidate ntest nfeatures nclasses sparse2dense kernel r k ϵ s t Tr V T], delim)
            end
        end
    end
    return nothing
end

function cv_error(df::DataFrame)
    # Group replicates based on hyperparameter pairs (ϵ, s).
    gdf = groupby(df, [:replicate, :epsilon, :sparsity])

    # Aggregate over folds.
    f(a,b,c,d) = (time=sum(a), train=mean(b), validation=mean(c), test=mean(d))
    out = combine(gdf, [:time, :train, :validation, :test] => f => AsTable)

    return out
end

function plot_cv_paths(df::DataFrame, col::Symbol)
    # Group by replicate and initialize figure based on selected column.
    gdf = groupby(df, :replicate)
    n = length(gdf)
    fig = plot(
        xlabel="Sparsity (%)",
        ylabel=col == :time ? "Time (s)" : "Error (%)",
        ylim=col == :time ? nothing : (-1,101),
        xlim=(-1,101),
        xticks=0:10:100,
        yticks=col == :time ? nothing : 0:10:100,
    )

    # Plot paths from each replicate.
    foreach(path -> plot!(path[!,:sparsity], path[!,col], lw=3, color=:black, alpha=1/n, label=nothing), gdf)

    # Highlight the mean and median paths.
    f(a) = (; median=median(a), mean=mean(a))
    tmp = combine(groupby(df, [:epsilon, :sparsity]), [col] => f => AsTable)
    plot!(tmp.sparsity, tmp.mean, lw=3, color=:red, ls=:dash, label=nothing)   # mean
    plot!(tmp.sparsity, tmp.median, lw=3, color=:blue, ls=:dot, label=nothing) # median

    return fig
end

"""
Returns the row index `j` corresponding to the optimal model.

The input `df` must contain cross-validation errors (see [`MVDA.cv_error`](@ref)).
Optimality is determined by the following:

- Robustness: maximal deadzone radius, `epsilon`.
- Parsimony: maximal `sparsity`.
- Predictive: minimal classification error in `validation`.
"""
function optimal_model(df::AbstractDataFrame)
    itr = zip(df.validation, df.sparsity, df.epsilon)
    adjusted_score = [(error, 100-s, 1/ϵ) for (error, s, ϵ) in itr]
    j = argmin(adjusted_score)
end

function credible_intervals(df::DataFrame, credibility=19/20)
    # Identify the optimal point in each replicate.
    gdf = groupby(df, :replicate)
    s_opt = zeros(length(gdf))
    for (r, replicate) in enumerate(gdf)
        j = optimal_model(replicate)
        s_opt[r] = replicate.sparsity[j]
    end

    # Compute parameter for equal-tailed interval and define functions to aggregate along path.
    α = (1 - credibility) / 2
    estimate_interval(data, _α) = median(data), quantile(data, _α), quantile(data, 1-_α)
    f = function (a, b, c, d)
        time_md, time_lo, time_hi = estimate_interval(a, α)
        train_md, train_lo, train_hi = estimate_interval(b, α)
        validation_md, validation_lo, validation_hi = estimate_interval(c, α)
        test_md, test_lo, test_hi = estimate_interval(d, α)
        return (;
            time_md=time_md, time_lo=time_lo, time_hi=time_hi,
            train_md=train_md, train_lo=train_lo, train_hi=train_hi,
            validation_md=validation_md, validation_lo=validation_lo, validation_hi=validation_hi,
            test_md=test_md, test_lo=test_lo, test_hi=test_hi,
        )
    end

    # Group by hyperparameter pairs (ϵ, s) and aggregate over replicates.
    out = combine(groupby(df, [:epsilon, :sparsity]), [:time, :train, :validation, :test] => f => AsTable)

    # Add optimal point credible interval to DataFrame.
    model_md, model_lo, model_hi = estimate_interval(s_opt, α)
    out[!, :model_md] .= model_md
    out[!, :model_lo] .= model_lo
    out[!, :model_hi] .= model_hi

    return out
end

function plot_credible_intervals(df::DataFrame, col::Symbol)
    # Plot the credible interval for the selected metric.
    ys, lo, hi = df[!,Symbol(col, :_md)], df[!,Symbol(col, :_lo)], df[!,Symbol(col, :_hi)]
    xs = df.sparsity
    lower, upper = ys - lo, hi - ys

    fig = plot(
        xlabel="Sparsity (%)",
        ylabel=col == :time ? "Time (s)" : "Error (%)",
        ylim=col == :time ? nothing : (-1,101),
        xlim=(-1,101),
        xticks=0:10:100,
    )

    # Add a point highlighting the optimal point + its credible interval.
    s_opt, s_lo, s_hi = df[1, [:model_md, :model_lo, :model_hi]]
    j = findlast(≤(s_opt), xs)
    error_bars = [(s_opt - s_lo, s_hi - s_opt)]
    scatter!((s_opt, ys[j]), xerr=error_bars, color=:black, markersize=6, markerstrokewidth=3, label="optimal model")

    annotate!([ (0.0, 90.0, ("Sparsity: $(round(s_opt, digits=4))%", 10, :left)) ])
    annotate!([ (0.0, 80.0, ("Error: $(round(ys[j], digits=4))%", 10, :left)) ])

    plot!(xs, ys, lw=3, ribbon=(lower, upper), label="95% credible interval", ls=:dash)

    return fig
end


text:
@inline fadj(g::MetaGraph, x...) = fadj(g.graph, x...)
@inline badj(g::MetaGraph, x...) = badj(g.graph, x...)

eltype(g::MetaGraph) = eltype(g.graph)
edgetype(g::MetaGraph) = edgetype(g.graph)
nv(g::MetaGraph) = nv(g.graph)
vertices(g::MetaGraph) = vertices(g.graph)

ne(g::MetaGraph) = ne(g.graph)
edges(g::MetaGraph) = edges(g.graph)

has_vertex(g::MetaGraph, x...) = has_vertex(g.graph, x...)
@inline has_edge(g::MetaGraph, x...) = has_edge(g.graph, x...)

inneighbors(g::MetaGraph, v::Integer) = inneighbors(g.graph, v)
outneighbors(g::MetaGraph, v::Integer) = fadj(g.graph, v)  # TODO: why not outneighbors?

issubset(g::T, h::T) where {T<:MetaGraph} = issubset(g.graph, h.graph)

"""
    add_edge!(g, u, v, val)

Add an edge `(u, v)` to MetaGraph `g` having value `val`.

Return true if the edge has been added, false otherwise.
"""
function add_edge!(g::MetaGraph, u::Integer, v::Integer, val)
    added = add_edge!(g.graph, u, v)
    if added
        metaindex = g.metaindex
        g.eprops[arrange(g, metaindex[u], metaindex[v], u, v)] = val
    end
    added
end

@inline function rem_edge!(g::MetaGraph, u::Integer, v::Integer)
    metaindex = g.metaindex
    removed = rem_edge!(g.graph, u, v)
    if removed
        delete!(g.eprops, arrange(g, metaindex[u], metaindex[v], u, v))
    end
    removed
end

"""
    add_vertex!(g, label, val)

Add a vertex to MetaGraph `g` with label `label` having value `val`.

Return true if the vertex has been added, false otherwise.
"""
function add_vertex!(g::MetaGraph, label, val)
    added = add_vertex!(g.graph)
    if added
        v = nv(g)
        g.vprops[label] = (v, val)
        g.metaindex[v] = label
    end
    added
end

function _rem_vertex!(g, label, v)
    vprops = g.vprops
    eprops = g.eprops
    metaindex = g.metaindex
    lastv = nv(g)
    for n in outneighbors(g, v)
        delete!(eprops, arrange(g, label, metaindex[n], v, n))
    end
    for n in inneighbors(g, v)
        delete!(eprops, arrange(g, metaindex[n], label, n, v))
    end
    removed = rem_vertex!(g.graph, v)
    if removed
        if v != lastv # ignore if we're removing the last vertex.
            lastl = metaindex[lastv]
            _, lastvprops = vprops[lastl]
            vprops[lastl] = v, lastvprops
            metaindex[v] = lastl
        end
        delete!(vprops, label)
        delete!(metaindex, lastv)
    end
    removed
end

function rem_vertex!(g::MetaGraph, v::Integer)
    exists = has_vertex(g, v)
    if exists
        _rem_vertex!(g, g.metaindex[v], v)
    else
        false
    end
end

function induced_subgraph(
    g::T,
    v::AbstractVector{U},
) where {T<:MetaGraph} where {U<:Integer}
    inducedgraph, vmap = induced_subgraph(g.graph, v)
    newg = MetaGraph(
        inducedgraph,
        empty(g.vprops),
        empty(g.eprops),
        g.gprops,
        g.weightfunction,
        g.defaultweight,
        empty(g.metaindex),
    )
    _copy_props!(g, newg, vmap)
    return newg, vmap
end


text:
module DiagTree
using PyCall, Printf
# export PropagatorKT, Weight, addChild
export Diagrams, Momentum, Propagator, addMomentum!, addPropagator!, Node, addChild

mutable struct Momentum
    basis::Vector{Int} #loop basis of the momentum
    curr::Vector{Float64}
    new::Vector{Float64}

    version::Int128
    excited::Bool #if set to excited, then the current weight needs to be replaced with the new weight

    function Momentum(loopbasis)
        _current = loopbasis .* 0
        _new = _current
        return new(loopbasis, _current, _new, 1, false)
    end
end

mutable struct Propagator{W}
    type::Int #1: Green's function, 2: interaction
    order::Int #the propagator may have an internal order (say, a Green's function diagram with multiple self-energy sub-diagrams)
    Kidx::Int #loop basis of the momentum
    Tidx::Tuple{Int,Int}
    factor::Float64 #additional factor, for example, if the propagator is an exchange interaction, then factor = -1.0

    version::Int128
    excited::Bool #if set to excited, then the current weight needs to be replaced with the new weight
    curr::W
    new::W

    function Propagator{W}(_type::Int, _order::Int, _Kidx, _Tidx, _factor = 1.0) where {W}
        return new{W}(_type, _order, _Kidx, Tuple(_Tidx), _factor, 0, false, W(0), W(0))
    end
end

mutable struct Node{W}
    id::Int
    operation::Int #1: multiply, 2: add
    factor::Float64 #symmetry factor, Fermi factor, spin factor
    version::Int128
    excited::Bool #if set to excited, then the current weight needs to be replaced with the new weight
    extK::Vector{Int}
    extT::Vector{Int}
    curr::W
    new::W

    #### link to the other nodes ##########
    parent::Int
    propagators::Vector{Int}
    nodes::Vector{Int} #if the Node is a leaf, then child stores the index of propagator, otherwise, it stores the indices of the child nodes

    function Node{W}(_parent = 0) where {W}
        new{W}(1, 0, 1.0, 0, false, [], [], W(0), W(0), _parent, [], [])
    end
    function Node{W}(id, operation::Int, factor, propagators = [], nodes = []; extK = [], extT = [], parent = 0) where {W}
        new{W}(id, operation, factor, 0, false, collect(extK), collect(extT), W(0), W(0), parent, propagators, nodes)
    end
end

isleaf(node) = (length(node.nodes) == 0) #if a node doesn't have child nodes, it must be a leaf (it must be ponting to propagators)

mutable struct Diagrams{W}
    momenta::Vector{Momentum}
    propagators::Vector{Propagator}
    tree::Vector{Node{W}}
    root::Vector{Int} #index of the root of the diagram tree
    # loopNum::Int
    # tauNum::Int
    function Diagrams{W}() where {W}
        momenta = Vector{Momentum}(undef, 0)
        propagators = Vector{Propagator}(undef, 0)
        tree = Vector{Node{W}}(undef, 0)
        return new(momenta, propagators, tree, [])
    end
end

Base.show(io::IO, w::Node{W}) where {W} = print(io, "Type: $(w.type), Prop: $(w.propagatorIdx), curr: $(w.current)")

function addChild(tree::Vector{Node{W}}, _parent) where {W}
    push!(tree, Node{W}(_parent))
    idx = length(tree) #the index of the last element, which is the new weight node
    push!(tree[_parent].child, idx)
    return idx
end

function addNode!(diagrams::Diagrams{W}, operation::Int, factor, propagators = [], nodes = []; isRoot = false, extT = [], extK = []) where {W}
    if isempty(propagators) && isempty(nodes)
        return -1
    end
    tree = diagrams.tree
    idx = length(tree) + 1
    # make sure that the propagators already exist
    for p in propagators
        @assert (1 <= p <= length(diagrams.propagators)) || p == -1
    end
    # make sure that the nodes already exist
    for n in nodes
        @assert (1 <= n <= length(diagrams.tree)) || n == -1
    end
    node = Node{W}(idx, operation, factor, propagators, nodes; extT = extT, extK = extK)
    push!(tree, node)
    if isRoot
        push!(diagrams.root, idx)
    end
    return idx #index of the new node
end

compareTidx(Tidx1, Tidx2, hasTimeReversal) = hasTimeReversal ? ((Tidx1 == Tidx2) || (Tidx1 == (Tidx2[2], Tidx2[1]))) : Tidx1 == Tidx2
compareKidx(Kidx1, Kidx2, hasMirrorSymmetry) = hasMirrorSymmetry ? ((Kidx1 ≈ Kidx2) || (Kidx1 ≈ -Kidx2)) : Kidx1 ≈ Kidx2

function addMomentum!(diagrams::Diagrams, _Kbasis, _symmetry)
    momenta = diagrams.momenta
    for (i, K) in enumerate(momenta)
        if compareKidx(K.basis, _Kbasis, :mirror in _symmetry)
            return i, false #existing momentum
        end
    end
    push!(momenta, Momentum(_Kbasis))
    return length(momenta), true #new momentum
end

# add new propagators to the propagator list
function addPropagator!(diagrams::Diagrams{W}, type::Int, order::Int, _Kbasis, _Tidx, symmetry = [], factor = 1.0) where {W}
    propagators = diagrams.propagators
    _Kidx, isNewK = addMomentum!(diagrams, _Kbasis, symmetry)

    for (i, p) in enumerate(propagators)
        if p.type == type && p.order == order
            Tflag = compareTidx(p.Tidx, _Tidx, (:particlehole in symmetry || :timereversal in symmetry))
            if Tflag && p.Kidx == _Kidx && p.factor ≈ factor
                return i, false #existing propagator
            end
        end
    end
    push!(propagators, Propagator{W}(type, order, _Kidx, _Tidx, factor))
    return length(propagators), true #new propagator
end

# function addPropagator!(diagrams::Diagrams{W}, type::Int, order::Int, _Kbasis, _Tidx::Vector{Tuple{Int,Int}}, _symmetry = [], _factor = 1.0) where {W}
#     @assert length(_Kbasis) == length(_Tidx)
#     pidx = []
#     for (ti, tidx) in enumerate(_Tidx)
#         kbasis = _Kbasis[ti]
#     end
# end

"""
    showTree(ver4, para::Para; verbose=0, depth=999)

    Visualize the diagram tree using ete3 python package

#Arguments
- `ver4`: the 4-vertex diagram tree to visualize
- `para`: parameters
- `verbose=0`: the amount of information to show
- `depth=999`: deepest level of the diagram tree to show
"""
function showTree(diag::Diagrams, _root = diag.root[end]; verbose = 0, depth = 999)
    tree = diag.tree
    K = diag.momenta
    root = tree[_root]

    # pushfirst!(PyVector(pyimport("sys")."path"), @__DIR__) #comment this line if no need to load local python module
    ete = PyCall.pyimport("ete3")

    function info(node)
        s = "N$(node.id):"
        if isempty(node.extT) == false
            s *= "T$(Tuple(node.extT)), "
        end

        if node.operation == 1
            if (node.factor ≈ 1.0) == false
                s *= @sprintf("%3.1f", node.factor)
            end
            s *= "× "
        elseif node.operation == 2
            @assert node.factor ≈ 1.0
            s *= "+ "
        else
            error("not implemented!")
        end
        return s
    end

    function treeview(node, level, t = nothing)
        if isnothing(t)
            t = ete.Tree(name = " ")
        end

        nt = t.add_child(name = info(node))
        name_face = ete.TextFace(nt.name, fgcolor = "black", fsize = 10)
        nt.add_face(name_face, column = 0, position = "branch-top")

        for child in node.nodes
            if child != -1
                treeview(tree[child], level + 1, nt)
            else
                nnt = nt.add_child(name = "0")
            end
        end
        for pidx in node.propagators
            if pidx != -1
                p = diag.propagators[pidx]
                factor = @sprintf("%3.1f", p.factor)
                nnt = nt.add_child(name = "P$pidx: typ $(p.type), K$(K[p.Kidx].basis), T$(p.Tidx), $factor")
            else
                nnt = nt.add_child(name = "0")
            end

            # name_face = ete.TextFace(nnt.name, fgcolor = "black", fsize = 10)
            # nnt.add_face(name_face, column = 0, position = "branch-top")
        end

        return t
    end

    t = treeview(root, 1)
    # style = ete.NodeStyle()
    # style["bgcolor"] = "Khaki"
    # t.set_style(style)

    ts = ete.TreeStyle()
    ts.show_leaf_name = true
    # ts.show_leaf_name = True
    # ts.layout_fn = my_layout
    ####### show tree vertically ############
    # ts.rotation = 90 #show tree vertically

    ####### show tree in an arc  #############
    # ts.mode = "c"
    # ts.arc_start = -180
    # ts.arc_span = 180
    # t.write(outfile="/home/kun/test.txt", format=8)
    t.show(tree_style = ts)
end

include("treeEval.jl")
end


text:
#=
- button http://jsxgraph.org/docs/symbols/Button.html
- checkbox http://jsxgraph.org/docs/symbols/Checkbox.html
=#

mutable struct Slider <: Object
    name::String
    vals::Vector{Vector{<:Real}}
    opts::Option{LittleDict{Symbol,Any}}
end

"""
    slider(name, vals; opts...)

Create a new slider, add it to the current board.

* name - name of the slider ex: `"slider1"`
* vals - position and values of the slider ex: `[[0,0],[3,0],[0,1.5,3]]`
          first element is `[xa, ya]` (position of minimum value point)
          second element is `[xb, yb]` (position of maximum value point)
          third element is `[va, v0, vb]` (min value, default value, max value)
"""
function slider(name, vals; kw...)
    @assert length(vals) == 3 "`slider`::expected three subarrays in `vals`."
    @assert length(vals[1]) == 2 && length(vals[2]) == 2 &&
            length(vals[3]) == 3 "`slider`::subvector dims should be 2,2,3."
    s = Slider(name, vals, dict(;kw...))
    return s
end
slider(name, a, b, v; kw...) = slider(name, [a, b, v]; kw...)

mutable struct Button{X<:FR,Y<:FR,F<:JSFun} <: Object
    name::String
    label::String
    x::X
    y::Y
    f::F
    opts::Option{LittleDict{Symbol,Any}}
end

"""
    button(name, label, x, y, f; opts...)
    button(label, x, y, f; opts)
"""
button(n, l, x, y, f; kw...) = Button(n, l, x, y, f, dict(;kw...))
button(l, x, y, f; kw...) = button("button_"*randstring(3), l, x, y, f; kw...)

# ---------------------------------------------------------------------------

val(s::Slider) = s.vals[3][2] # midvalue of last vector [va, v0, vb]

# ---------------------------------------------------------------------------

function str(s::Slider, b::Board)
    opts = get_opts(s)
    jss = js".create('slider', $(s.vals), $opts);"
    return "var $(s.name)=" * b.name * jss.s
end

function str(bu::Button, b::Board)
    xs, xss, xrp = strf(bu.x, "FX", b)
    ys, yss, yrp = strf(bu.y, "FY", b)
    fs, fss, frp = strf(bu.f, "FF", b)
    opts = get_opts(bu)
    jss = js".create('button', [$xss, $yss, $(bu.label), $fss], $opts);"
    return xs * ys * fs * b.name * replacefn(jss.s, xrp, yrp, frp)
end


text:
# The functions in this module are a direct translation of the original
# Matlab functions by J. J. Moré and S. M. Wild.
# See http://www.mcs.anl.gov/~more/dfo.
#
# D. Orban, 2016.

"""Compute data ratios used to produce a data profile.

There is normally no need to call this function directly.
See the documentation of `data_profile()` for more information.
"""
function data_ratios(H :: Array{Float64,3}, N :: Vector{Float64}; τ :: Float64=1.0e-3)

  (nf, np, ns) = size(H)
  H[isinf(H)] = NaN;
  H[H .< 0] = NaN;
  for j = 1 : ns
    for i = 2 : nf
      H[i, :, j] = min(H[i, :, j], H[i-1, :, j])
    end
  end

  prob_min = minimum(minimum(H, 1), 3)  # min for each problem
  prob_max = H[1, :, 1]      # starting value for each problem

  # For each problem and solver, determine the number of costly operations
  # (e.g, objective evaluations) required to reach the cutoff value.
  T = zeros(np, ns)
  for p = 1 : np
    cutoff = prob_min[p] + τ * (prob_max[p] - prob_min[p])
    for s = 1 : ns
      nfevs = findfirst(H[:, p, s] .≤ cutoff)
      T[p, s] = nfevs > 0 ? nfevs/N[p] : NaN
    end
  end

  # Replace all NaNs with twice the max ratio and sort.
  max_data = maximum(T)
  T[isnan(T)] = 2 * max_data
  T = sort(T, 1)
  return (T, max_data)
end

data_ratios{TH <: Number, TN <: Number, Tt <: Number}(H :: Array{TH,3}, N :: Vector{TN}; τ :: Tt=convert(Tt, 1.0e-3)) = data_ratios(convert(Array{Float64,3}, H), convert(Vector{Float64}, N), τ=convert(Float64, τ))

"""Produce a data profile.

The 3-dimensional array `H` gives the performance data for each solver
and each problem (smaller is better). `H[k,p,s]` is the `k`-th costly
operation (e.g., function evaluation) for problem `p` and solver `s`.
Failures on a given problem are represented by a negative value, an
infinite value, or `NaN`.

The vector `N` gives a scaling associated to each problem. If the number
of simplex gradients is being measured, `N[p]` should be `n(p) + 1`
where `n(p)` is the number of variables of problem `p`.

The threshold `τ` determines the tolerance in the convergence criterion

  f(x) ≤ fL + τ (f0 - fL),

where for each problem f(x) is a measure recorded (e.g., the objective
value), f0 is the measure at the initial point, and fL is the best measure
obtained by any solver.
"""
function data_profile(H :: Array{Float64,3}, N :: Vector{Float64},
                      labels :: Vector{AbstractString}=AbstractString[];
                      τ :: Float64=1.0e-3,
                      operations :: AbstractString="function evaluations",
                      title :: AbstractString="",
                      kwargs...)

  (T, max_data) = data_ratios(H, N, τ=τ)
  (np, ns) = size(T)
  xs = [1:np;] / np
  length(labels) == 0 && (labels = [@sprintf("column %d", col) for col = 1 : ns])
  profile = Plots.plot(; kwargs...)  # initial empty plot
  for s = 1 : ns
    Plots.plot!(T[:, s], xs, t=:steppost, label=labels[s])
  end
  Plots.xlims!(0.0, 1.1 * max_data)
  Plots.ylims!(0, 1.1)
  Plots.xlabel!("Number of " * strip(operations) * @sprintf(" (thresh %7.1e)", τ))
  Plots.ylabel!("Proportion of problems")
  Plots.title!(title)
  return profile
end

data_profile{TH <: Number, TN <: Number, S <: AbstractString}(H :: Array{TH,3}, N :: Vector{TN}, labels :: Vector{S}; kwargs...) =
  data_profile(convert(Array{Float64,3}, H), convert(Vector{Float64}, N), convert(Vector{AbstractString}, labels); kwargs...)


text:
@node Gamma Stochastic [ out, (α, aliases = [ shape ]), (θ, aliases = [ scale ]) ]

@average_energy Gamma (q_out::Any, q_α::PointMass, q_θ::Any) = begin
    loggammamean(q_α) + mean(q_α) * logmean(q_θ) - (mean(q_α) - 1.0) * logmean(q_out) + mean(q_out)/mean(q_θ)
end

@average_energy Gamma (q_out::Any, q_α::GammaDistributionsFamily, q_θ::Any) = begin
    loggammamean(q_α) + mean(q_α) * logmean(q_θ) - (mean(q_α) - 1.0) * logmean(q_out) + mean(q_out)/mean(q_θ)
end