text:
### A Pluto.jl notebook ###
# v0.15.1

using Markdown
using InteractiveUtils

# ╔═╡ 1e037376-8252-41fb-ae95-93bc9011b7fa
# we need to install a special package in order to use the gamma function
using SpecialFunctions

# ╔═╡ 76d6c2d4-e370-4e83-9d04-ca6fc546bdd9
# we write a function that returns the exact volume of a hypersphere in d dimensions
function exactVol(d)
    vol = pi^(d/2) / gamma(1 + d/2)
    return vol
end

# ╔═╡ b58cf6d8-291c-49f4-a61a-51f1a7454651
n = 10000;

# ╔═╡ 2f4a8b5e-f7c9-4abc-b04d-3286207a97b3
dim = 20;

# ╔═╡ 4da795ef-b8d3-4c94-bb89-7eeaca9602fb
md"Exact volume is:"

# ╔═╡ 8d8e10d0-2597-4198-92b1-9ab2ebe4b3c8
exV = exactVol(dim)

# ╔═╡ 6a8c6c37-2a48-47c7-8d40-ca276c8a09d2
# we write a function that returns a random vector from [0,1]^d with exponential density
function randExp(d, lam)

	# we generate a uniformly random vector
	y = rand(d)

	# we transform (elementwise) to random numbers in [exp(-lambda),1]
	y = exp(-lam) .+ (1-exp(-lam))*y

	# we now get numbers with exponential density
	x = -log.(y)/lam
	return x
end

# ╔═╡ 8de1cad7-ac2c-42ec-9581-1d7ff06d576d
# we write an improoved function that uses Monte Carlo with exponential importance sampling
# to estimate the volume of the unit hypersphere in d dimensions.
# We choose our points x with independent Cartesian components and prob density exp(-lam.xi) / [ lam*((1-exp(-lam)) ].
# We achieve this by taking xi = (1/lam) log (1/y) where y is chosen uniformly in (exp(-lam),1)
function sphereVolMCImpExp(d, n, lam)
	count = 0.0
	for i = 1:n

		# we generate x according to our distribution
		x = randExp(d, lam)

		# we calculate the squared norm of x and check if x is inside the hypersphere
		norm2 = x' * x
		if norm2 < 1

			# we compute the density at x
			rho = 1
			for i = 1:d

				# note that y(i) = exp(-lam * x(i))
				rho *= exp(-lam*x[i]) * lam / (1-exp(-lam))
			end
			count += 1/rho
		end
	end
	vol =  2^d * count / n
	return vol
end

# ╔═╡ c8381d29-3b9f-4399-9a60-c0b6b620e99c
sphereVolMCImpExp(20, n, 2)

# ╔═╡ 0e1580e1-cdae-498e-ab14-261f8a96a802
md" We do the calculation $m$ times and check mean/variance:"

# ╔═╡ 770e2019-420d-415b-b498-96493ac676bf
m = 1000;

# ╔═╡ dc6e6729-7ce5-45a8-8ecc-b320d283340a
data = [sphereVolMCImpExp(20, n, 2) for i in 1:m];

# ╔═╡ fb6c8ff3-d200-4e48-b9ad-aa34d188e675
mean = sum(data)/m

# ╔═╡ 48b51297-b3fb-4bce-8d41-7af6120b553c
error = mean - exV

# ╔═╡ bb18d766-cdad-45b7-aebb-c59763896bd4
variance = sum(broadcast(x -> (x-mean)^2, data))/m

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"

[compat]
SpecialFunctions = "~1.6.1"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[ChainRulesCore]]
deps = ["Compat", "LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "bdc0937269321858ab2a4f288486cb258b9a0af7"
uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
version = "1.3.0"

[[Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "727e463cfebd0c7b999bbf3e9e7e16f254b94193"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.34.0"

[[CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[DocStringExtensions]]
deps = ["LibGit2"]
git-tree-sha1 = "a32185f5428d3986f47c2ab78b1f216d5e6cc96f"
uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
version = "0.8.5"

[[Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[IrrationalConstants]]
git-tree-sha1 = "f76424439413893a832026ca355fe273e93bce94"
uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
version = "0.1.0"

[[JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "642a199af8b68253517b80bd3bfd17eb4e84df6e"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.3.0"

[[LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[LinearAlgebra]]
deps = ["Libdl"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[LogExpFunctions]]
deps = ["DocStringExtensions", "IrrationalConstants", "LinearAlgebra"]
git-tree-sha1 = "3d682c07e6dd250ed082f883dc88aee7996bf2cc"
uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
version = "0.3.0"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[OpenSpecFun_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
version = "0.5.5+0"

[[Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[Preferences]]
deps = ["TOML"]
git-tree-sha1 = "00cfd92944ca9c760982747e9a1d0d5d86ab1e5a"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.2.2"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[SpecialFunctions]]
deps = ["ChainRulesCore", "LogExpFunctions", "OpenSpecFun_jll"]
git-tree-sha1 = "a322a9493e49c5f3a10b50df3aedaf1cdb3244b7"
uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
version = "1.6.1"

[[Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╠═1e037376-8252-41fb-ae95-93bc9011b7fa
# ╠═76d6c2d4-e370-4e83-9d04-ca6fc546bdd9
# ╠═b58cf6d8-291c-49f4-a61a-51f1a7454651
# ╠═2f4a8b5e-f7c9-4abc-b04d-3286207a97b3
# ╟─4da795ef-b8d3-4c94-bb89-7eeaca9602fb
# ╠═8d8e10d0-2597-4198-92b1-9ab2ebe4b3c8
# ╠═6a8c6c37-2a48-47c7-8d40-ca276c8a09d2
# ╠═8de1cad7-ac2c-42ec-9581-1d7ff06d576d
# ╠═c8381d29-3b9f-4399-9a60-c0b6b620e99c
# ╟─0e1580e1-cdae-498e-ab14-261f8a96a802
# ╠═770e2019-420d-415b-b498-96493ac676bf
# ╠═dc6e6729-7ce5-45a8-8ecc-b320d283340a
# ╠═fb6c8ff3-d200-4e48-b9ad-aa34d188e675
# ╠═48b51297-b3fb-4bce-8d41-7af6120b553c
# ╠═bb18d766-cdad-45b7-aebb-c59763896bd4
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

max_stars_repo_path: 1a/Monte Carlo integration/mcTest.jl

max_issues_repo_path: 1a/Monte Carlo integration/mcTest.jl

max_forks_repo_path: 1a/Monte Carlo integration/mcTest.jl


text:
################################################################################
#                             BASIC REFERENCE METHODS
################################################################################
"""
    JuMP.owner_model(cref::InfOptConstraintRef)::InfiniteModel

Extend [`JuMP.owner_model`](@ref JuMP.owner_model(::JuMP.ConstraintRef)) to
return the infinite model associated with `cref`.

**Example**
```julia-repl
julia> model = owner_model(cref)
An InfiniteOpt Model
Minimization problem with:
Finite Parameters: 0
Infinite Parameters: 3
Variables: 3
Derivatives: 0
Measures: 0
Objective function type: GeneralVariableRef
`GenericAffExpr{Float64,GeneralVariableRef}`-in-`MathOptInterface.EqualTo{Float64}`: 1 constraint
Names registered in the model: g, t, h, x
Optimizer model backend information:
Model mode: AUTOMATIC
CachingOptimizer state: NO_OPTIMIZER
Solver name: No optimizer attached.
```
"""
JuMP.owner_model(cref::InfOptConstraintRef)::InfiniteModel = cref.model

"""
    JuMP.index(cref::InfOptConstraintRef)::ConstraintIndex

Extend [`JuMP.index`](@ref JuMP.index(::JuMP.ConstraintRef)) to return the index
of an `InfiniteOpt` constraint `cref`.

**Example**
```julia-repl
julia> index(cref)
ConstrainIndex(2)
```
"""
JuMP.index(cref::InfOptConstraintRef)::ConstraintIndex = cref.index

# Extend Base and JuMP functions
function Base.:(==)(v::InfOptConstraintRef, w::InfOptConstraintRef)::Bool
    return v.model === w.model && v.index == w.index
end
Base.broadcastable(cref::InfOptConstraintRef) = Ref(cref)
JuMP.constraint_type(m::InfiniteModel) = InfOptConstraintRef

################################################################################
#                             CORE OBJECT METHODS
################################################################################
# Extend _add_data_object
function _add_data_object(model::InfiniteModel,
                          object::ConstraintData)::ConstraintIndex
    return MOIUC.add_item(model.constraints, object)
end

# Extend _data_dictionary
function _data_dictionary(cref::InfOptConstraintRef
    )::MOIUC.CleverDict{ConstraintIndex, ConstraintData}
    return JuMP.owner_model(cref).constraints
end

# Extend _data_object
function _data_object(cref::InfOptConstraintRef)::ConstraintData
  object = get(_data_dictionary(cref), JuMP.index(cref), nothing)
  object === nothing && error("Invalid constraint reference, cannot find " *
                        "corresponding constraint in the model. This is likely " *
                        "caused by using the reference of a deleted constraint.")
  return object
end

# Return the core constraint object
function _core_constraint_object(cref::InfOptConstraintRef)::JuMP.AbstractConstraint
    return _data_object(cref).constraint
end

# Update the core constraint object
function _set_core_constraint_object(cref::InfOptConstraintRef,
                                     constr::JuMP.AbstractConstraint)::Nothing
    _data_object(cref).constraint = constr
    return
end

# Extend _object_numbers
function _object_numbers(cref::InfOptConstraintRef)::Vector{Int}
    return _data_object(cref).object_nums
end

# Extend _measure_dependencies
function _measure_dependencies(cref::InfOptConstraintRef)::Vector{MeasureIndex}
    return _data_object(cref).measure_indices
end

# Return if this constraint is an info constraint
function _is_info_constraint(cref::InfOptConstraintRef)::Bool
    return _data_object(cref).is_info_constraint
end

# Extend _delete_data_object
function _delete_data_object(cref::InfOptConstraintRef)::Nothing
    delete!(_data_dictionary(cref), JuMP.index(cref))
    return
end

################################################################################
#                             DEFINITION METHODS
################################################################################
# This might not be necessary...
function JuMP.build_constraint(_error::Function,
    v::GeneralVariableRef,
    set::MOI.AbstractScalarSet;
    parameter_bounds::ParameterBounds{GeneralVariableRef} = ParameterBounds()
    )::JuMP.AbstractConstraint
    # make the constraint
    if length(parameter_bounds) != 0
        _check_bounds(parameter_bounds)
        return BoundedScalarConstraint(v, set, parameter_bounds,
                                       copy(parameter_bounds))
    else
        return JuMP.ScalarConstraint(v, set)
    end
end

"""
    JuMP.build_constraint(_error::Function, expr::JuMP.AbstractJuMPScalar,
                          set::MOI.AbstractScalarSet;
                          [parameter_bounds::ParameterBounds = ParameterBounds()]
                          )::JuMP.AbstractConstraint

Extend `JuMP.build_constraint` to accept the ```parameter_bounds``` argument
and return a [`BoundedScalarConstraint`](@ref) if the ```parameter_bounds``` keyword
argument is specifed or return a [`JuMP.ScalarConstraint`](@ref) otherwise. This is
primarily intended to work as an internal function for constraint macros.

**Example**
```jldoctest; setup = :(using JuMP, InfiniteOpt; model = InfiniteModel())
julia> @infinite_parameter(model, t in [0, 10]);

julia> @infinite_variable(model, g(t));

julia> @hold_variable(model, x);

julia> constr = build_constraint(error, g + x, MOI.EqualTo(42.0),
              parameter_bounds = ParameterBounds((t => IntervalSet(0, 1),)));

julia> isa(constr, BoundedScalarConstraint)
true
```
"""
function JuMP.build_constraint(_error::Function,
    expr::Union{JuMP.GenericAffExpr{C, V}, JuMP.GenericQuadExpr{C, V}},
    set::MOI.AbstractScalarSet;
    parameter_bounds::ParameterBounds{GeneralVariableRef} = ParameterBounds()
    )::JuMP.AbstractConstraint where {C, V <: GeneralVariableRef}
    # make the constraint
    offset = JuMP.constant(expr)
    JuMP.add_to_expression!(expr, -offset)
    if !isempty(parameter_bounds)
        _check_bounds(parameter_bounds)
        return BoundedScalarConstraint(expr, MOIU.shift_constant(set, -offset),
                                       parameter_bounds, copy(parameter_bounds))
    else
        return JuMP.ScalarConstraint(expr, MOIU.shift_constant(set, -offset))
    end
end

## Perfrom bound checks and update them if needed, then return the updated constraint
# BoundedScalarConstraint
function _check_and_update_bounds(model::InfiniteModel,
                                  c::BoundedScalarConstraint,
                                  vrefs::Vector{GeneralVariableRef}
                                  )::JuMP.AbstractConstraint
    # look for bounded hold variables and update bounds
    for vref in vrefs
        _update_var_bounds(vref, c.bounds)
    end
    # now validate and return
    _validate_bounds(model, c.bounds)
    # TODO should we check that bounds don't violate point variables?
    return c
end

# ScalarConstraint
function _check_and_update_bounds(model::InfiniteModel,
                                  c::JuMP.ScalarConstraint,
                                  vrefs::Vector{GeneralVariableRef}
                                  )::JuMP.AbstractConstraint
    bounds = ParameterBounds()
    # check for bounded hold variables and build the intersection of the bounds
    for vref in vrefs
        _update_var_bounds(vref, bounds)
    end
    # if we added bounds, change to a bounded constraint and validate
    if !isempty(bounds)
        new_c = BoundedScalarConstraint(c.func, c.set, bounds, ParameterBounds())
        _validate_bounds(model, new_c.bounds)
        return new_c
    end
    # TODO should we check that bounds don't violate point variables?
    return c
end

# Used to update the model.var_to_constrs field
function _update_var_constr_mapping(vrefs::Vector{GeneralVariableRef},
                                    cref::InfOptConstraintRef)::Nothing
    for vref in vrefs
        dvref = dispatch_variable_ref(vref)
        push!(_constraint_dependencies(dvref), JuMP.index(cref))
        if dvref isa MeasureRef
            push!(_measure_dependencies(cref), JuMP.index(dvref))
        end
        # TODO maybe make mapping for derivatives...
    end
    return
end

# Extend functions for bounded constraints
JuMP.shape(c::BoundedScalarConstraint) = JuMP.shape(JuMP.ScalarConstraint(c.func, c.set))
JuMP.jump_function(c::BoundedScalarConstraint) = c.func
JuMP.moi_set(c::BoundedScalarConstraint) = c.set
parameter_bounds(c::BoundedScalarConstraint) = c.bounds
original_parameter_bounds(c::BoundedScalarConstraint) = c.orig_bounds

"""
    JuMP.add_constraint(model::InfiniteModel, c::JuMP.AbstractConstraint,
                        [name::String = ""])::InfOptConstraintRef

Extend [`JuMP.add_constraint`](@ref JuMP.add_constraint(::JuMP.Model, ::JuMP.AbstractConstraint, ::String))
to add a constraint `c` to an infinite model
`model` with name `name`. Returns an appropriate constraint reference whose type
depends on what variables are used to define the constraint. Errors if a vector
constraint is used, the constraint only constains parameters, or if any
variables do not belong to `model`. This is primarily used as an internal
method for the cosntraint macros.

**Example**
```jldoctest; setup = :(using JuMP, InfiniteOpt; model = InfiniteModel())
julia> @infinite_parameter(model, t in [0, 10]);

julia> @infinite_variable(model, g(t));

julia> @hold_variable(model, x);

julia> constr = build_constraint(error, g + x, MOI.EqualTo(42));

julia> cref = add_constraint(model, constr, "name")
name : g(t) + x = 42.0, ∀ t ∈ [0, 10]
```
"""
function JuMP.add_constraint(model::InfiniteModel,
                             c::JuMP.AbstractConstraint,
                             name::String = "";
                             is_info_constr::Bool = false
                             )::InfOptConstraintRef
    # gather the unique list of variable references for testing and mapping
    vrefs = _all_function_variables(JuMP.jump_function(c))
    # test in the model
    for vref in vrefs
        JuMP.check_belongs_to_model(vref, model)
    end
    # test parameter bounds and update with hold variable bounds
    if model.has_hold_bounds
        c = _check_and_update_bounds(model, c, vrefs)
    elseif c isa BoundedScalarConstraint
        _validate_bounds(model, parameter_bounds(c))
    end
    # get the parameter object numbers
    object_nums = sort!(_object_numbers(vrefs))
    # add the constaint to the model
    constr_object = ConstraintData(c, object_nums, name, MeasureIndex[],
                                   is_info_constr)
    cindex = _add_data_object(model, constr_object)
    cref = InfOptConstraintRef(model, cindex)
    # update the variable mappings and model status
    _update_var_constr_mapping(vrefs, cref)
    set_optimizer_model_ready(model, false)
    return cref
end

# Fallback for vector constraints
function JuMP.add_constraint(model::InfiniteModel, c::JuMP.VectorConstraint,
                             name::String = "";
                             is_info_constr::Bool = false)
    error("Vector constraints not supported in InfiniteOpt.")
end

################################################################################
#                                 JuMP METHODS
################################################################################
"""
    JuMP.is_valid(model::InfiniteModel, cref::InfOptConstraintRef)::Bool

Extend [`JuMP.is_valid`](@ref JuMP.is_valid(::JuMP.Model, ::JuMP.ConstraintRef{JuMP.Model}))
to return `Bool` whether an `InfiniteOpt` constraint reference is valid.

**Example**
```jldoctest; setup = :(using JuMP, InfiniteOpt; model = InfiniteModel(); cref = @constraint(model, 2 == 2))
julia> is_valid(model, cref)
true
```
"""
function JuMP.is_valid(model::InfiniteModel, cref::InfOptConstraintRef)::Bool
    return (model === JuMP.owner_model(cref) &&
            JuMP.index(cref) in keys(_data_dictionary(cref)))
end

"""
    JuMP.constraint_object(cref::InfOptConstraintRef)::JuMP.AbstractConstraint

Extend [`JuMP.constraint_object`](@ref JuMP.constraint_object)
to return the constraint object associated with `cref`.

**Example**
```jldoctest; setup = :(using JuMP, InfiniteOpt; model = InfiniteModel())
julia> @infinite_parameter(model, t in [0, 10]);

julia> @hold_variable(model, x <= 1);

julia> cref = UpperBoundRef(x);

julia> obj = constraint_object(cref)
ScalarConstraint{GeneralVariableRef,MathOptInterface.LessThan{Float64}}(x,
MathOptInterface.LessThan{Float64}(1.0))
```
"""
function JuMP.constraint_object(cref::InfOptConstraintRef)::JuMP.AbstractConstraint
    return _core_constraint_object(cref)
end

"""
    JuMP.name(cref::InfOptConstraintRef)::String

Extend [`JuMP.name`](@ref JuMP.name(::JuMP.ConstraintRef{JuMP.Model,<:JuMP._MOICON})
to return the name of an `InfiniteOpt` constraint.

**Example**
```jldoctest; setup = :(using JuMP, InfiniteOpt; model = InfiniteModel(); cref = @constraint(model, constr_name, 2 == 2))
julia> name(cref)
"constr_name"
```
"""
function JuMP.name(cref::InfOptConstraintRef)::String
    object = get(_data_dictionary(cref), JuMP.index(cref), nothing)
    return object === nothing ? "" : object.name
end

"""
    JuMP.set_name(cref::InfOptConstraintRef, name::String)::Nothing

Extend [`JuMP.set_name`](@ref JuMP.set_name(::JuMP.ConstraintRef{JuMP.Model,<:JuMP._MOICON}, ::String))
to specify the name of a constraint `cref`.

**Example**
```jldoctest; setup = :(using JuMP, InfiniteOpt; model = InfiniteModel(); cref = @constraint(model, 2 == 2))
julia> set_name(cref, "new_name")

julia> name(cref)
"new_name"
```
"""
function JuMP.set_name(cref::InfOptConstraintRef, name::String)::Nothing
    _data_object(cref).name = name
    JuMP.owner_model(cref).name_to_constr = nothing
    return
end

# Return a constraint set with an updated value
function _set_set_value(set::S, value::Real) where {T, S <: Union{MOI.LessThan{T},
                                            MOI.GreaterThan{T}, MOI.EqualTo{T}}}
    return S(convert(T, value))
end

"""
    JuMP.set_normalized_rhs(cref::InfOptConstraintRef, value::Real)::Nothing

Set the right-hand side term of `constraint` to `value`.
Note that prior to this step, JuMP will aggregate all constant terms onto the
right-hand side of the constraint. For example, given a constraint `2x + 1 <=
2`, `set_normalized_rhs(con, 4)` will create the constraint `2x <= 4`, not `2x +
1 <= 4`.

```julia-repl
julia> @constraint(model, con, 2x + 1 <= 2)
con : 2 x ≤ 1.0

julia> set_normalized_rhs(con, 4)

julia> con
con : 2 x ≤ 4.0
```
"""
function JuMP.set_normalized_rhs(cref::InfOptConstraintRef, value::Real)::Nothing
    old_constr = JuMP.constraint_object(cref)
    new_set = _set_set_value(JuMP.moi_set(old_constr), value)
    if old_constr isa BoundedScalarConstraint
        new_constr = BoundedScalarConstraint(JuMP.jump_function(old_constr),
                        new_set, parameter_bounds(old_constr),
                        original_parameter_bounds(old_constr))
    else
        new_constr = JuMP.ScalarConstraint(JuMP.jump_function(old_constr),
                                           new_set)
    end
    _set_core_constraint_object(cref, new_constr)
    return
end

"""
    JuMP.normalized_rhs(cref::InfOptConstraintRef)::Float64

Return the right-hand side term of `cref` after JuMP has converted the
constraint into its normalized form.
"""
function JuMP.normalized_rhs(cref::InfOptConstraintRef)::Float64
    constr = JuMP.constraint_object(cref)
    return MOI.constant(JuMP.moi_set(constr))
end

"""
    JuMP.add_to_function_constant(cref::InfOptConstraintRef, value::Real)::Nothing

Add `value` to the function constant term.
Note that for scalar constraints, JuMP will aggregate all constant terms onto the
right-hand side of the constraint so instead of modifying the function, the set
will be translated by `-value`. For example, given a constraint `2x <=
3`, `add_to_function_constant(c, 4)` will modify it to `2x <= -1`.
```
"""
function JuMP.add_to_function_constant(cref::InfOptConstraintRef,
                                       value::Real)::Nothing
    current_value = JuMP.normalized_rhs(cref)
    JuMP.set_normalized_rhs(cref, current_value - value)
    return
end

"""
    JuMP.set_normalized_coefficient(cref::InfOptConstraintRef,
                                    variable::GeneralVariableRef,
                                    value::Real)::Nothing

Set the coefficient of `variable` in the constraint `constraint` to `value`.
Note that prior to this step, JuMP will aggregate multiple terms containing the
same variable. For example, given a constraint `2x + 3x <= 2`,
`set_normalized_coefficient(con, x, 4)` will create the constraint `4x <= 2`.

```julia-repl
julia> con
con : 5 x ≤ 2.0

julia> set_normalized_coefficient(con, x, 4)

julia> con
con : 4 x ≤ 2.0
```
"""
function JuMP.set_normalized_coefficient(cref::InfOptConstraintRef,
                                         variable::GeneralVariableRef,
                                         value::Real)::Nothing
    # update the constraint expression and update the constraint
    old_constr = JuMP.constraint_object(cref)
    new_expr = _set_variable_coefficient!(JuMP.jump_function(old_constr),
                                          variable, value)
    if old_constr isa BoundedScalarConstraint
        new_constr = BoundedScalarConstraint(new_expr, JuMP.moi_set(old_constr),
                        parameter_bounds(old_constr),
                        original_parameter_bounds(old_constr))
    else
        new_constr = JuMP.ScalarConstraint(new_expr, JuMP.moi_set(old_constr))
    end
    _set_core_constraint_object(cref, new_constr)
    return
end

"""
    JuMP.normalized_coefficient(cref::InfOptConstraintRef,
                                variable::GeneralVariableRef)::Float64

Return the coefficient associated with `variable` in `constraint` after JuMP has
normalized the constraint into its standard form.
"""
function JuMP.normalized_coefficient(cref::InfOptConstraintRef,
                                     variable::GeneralVariableRef)::Float64
    constr = JuMP.constraint_object(cref)
    expr = JuMP.jump_function(constr)
    if expr isa GeneralVariableRef && expr == variable
        return 1.0
    elseif expr isa GeneralVariableRef
        return 0.0
    else
        return JuMP._affine_coefficient(expr, variable)
    end
end

# Return the appropriate constraint reference given the index and model
function _make_constraint_ref(model::InfiniteModel,
                              index::ConstraintIndex)::InfOptConstraintRef
    constr = model.constraints[index].constraint
    return InfOptConstraintRef(model, index)
end

"""
    JuMP.constraint_by_name(model::InfiniteModel,
                            name::String)::Union{InfOptConstraintRef, Nothing}

Extend [`JuMP.constraint_by_name`](@ref JuMP.constraint_by_name)
to return the constraint reference
associated with `name` if one exists or returns nothing. Errors if more than
one constraint uses the same name.

**Example**
```julia-repl
julia> constraint_by_name(model, "constr_name")
constr_name : x + pt = 3.0
```
"""
function JuMP.constraint_by_name(model::InfiniteModel,
    name::String)::Union{InfOptConstraintRef, Nothing}
    if model.name_to_constr === nothing
        # Inspired from MOI/src/Utilities/model.jl
        model.name_to_constr = Dict{String, Int}()
        for (index, data_object) in model.constraints
            constr_name = data_object.name
            if haskey(model.name_to_constr, constr_name)
                model.name_to_constr[constr_name] = ConstraintIndex(-1)
            else
                model.name_to_constr[constr_name] = index
            end
        end
    end
    index = get(model.name_to_constr, name, nothing)
    if index === nothing
        return nothing
    elseif index == ConstraintIndex(-1)
        error("Multiple constraints have the name $name.")
    else
        return _make_constraint_ref(model, index)
    end
end

"""
    JuMP.num_constraints(model::InfiniteModel,
                         function_type::Type{<:JuMP.AbstractJuMPScalar},
                         set_type::Type{<:MOI.AbstractSet})::Int

Extend [`JuMP.num_constraints`](@ref JuMP.num_constraints(::JuMP.Model, ::Type{<:Union{JuMP.AbstractJuMPScalar, Vector{<:JuMP.AbstractJuMPScalar}}}, ::Type{<:MOI.AbstractSet}))
to return the number of constraints
with a partiuclar function type and set type.

**Example**
```julia-repl
julia> num_constraints(model, HoldVariableRef, MOI.LessThan)
1
```
"""
function JuMP.num_constraints(model::InfiniteModel,
                              function_type::Type{<:JuMP.AbstractJuMPScalar},
                              set_type::Type{<:MOI.AbstractSet})::Int
    counter = 0
    for (index, data_object) in model.constraints
        if isa(JuMP.jump_function(data_object.constraint), function_type) &&
           isa(JuMP.moi_set(data_object.constraint), set_type)
            counter += 1
        end
    end
    return counter
end

"""
    JuMP.num_constraints(model::InfiniteModel,
                         function_type::Type{<:JuMP.AbstractJuMPScalar})::Int

Extend [`JuMP.num_constraints`](@ref JuMP.num_constraints(::JuMP.Model, ::Type{<:Union{JuMP.AbstractJuMPScalar, Vector{<:JuMP.AbstractJuMPScalar}}}, ::Type{<:MOI.AbstractSet}))
to search by function types for all MOI
sets and return the total number of constraints with a particular function type.

```julia-repl
julia> num_constraints(model, HoldVariableRef)
3
```
"""
function JuMP.num_constraints(model::InfiniteModel,
                            function_type::Type{<:JuMP.AbstractJuMPScalar})::Int
    return JuMP.num_constraints(model, function_type, MOI.AbstractSet)
end

"""
    JuMP.num_constraints(model::InfiniteModel,
                         function_type::Type{<:MOI.AbstractSet})::Int

Extend [`JuMP.num_constraints`](@ref JuMP.num_constraints(::JuMP.Model, ::Type{<:Union{JuMP.AbstractJuMPScalar, Vector{<:JuMP.AbstractJuMPScalar}}}, ::Type{<:MOI.AbstractSet}))
to search by MOI set type for all function
types and return the total number of constraints that use a particular MOI set
type.

```julia-repl
julia> num_constraints(model, MOI.LessThan)
2
```
"""
function JuMP.num_constraints(model::InfiniteModel,
                              set_type::Type{<:MOI.AbstractSet})::Int
    return JuMP.num_constraints(model, JuMP.AbstractJuMPScalar, set_type)
end

"""
    JuMP.num_constraints(model::InfiniteModel)::Int

Extend [`JuMP.num_constraints`](@ref JuMP.num_constraints(::JuMP.Model, ::Type{<:Union{JuMP.AbstractJuMPScalar, Vector{<:JuMP.AbstractJuMPScalar}}}, ::Type{<:MOI.AbstractSet}))
to return the total number of constraints
in an infinite model `model`.

```julia-repl
julia> num_constraints(model)
4
```
"""
function JuMP.num_constraints(model::InfiniteModel)::Int
    return length(model.constraints)
end

"""
    JuMP.all_constraints(model::InfiniteModel,
                         function_type::Type{<:JuMP.AbstractJuMPScalar},
                         set_type::Type{<:MOI.AbstractSet}
                         )::Vector{InfOptConstraintRef}

Extend [`JuMP.all_constraints`](@ref JuMP.all_constraints(::JuMP.Model, ::Type{<:Union{JuMP.AbstractJuMPScalar, Vector{<:JuMP.AbstractJuMPScalar}}}, ::Type{<:MOI.AbstractSet}))
to return a list of all the constraints with a particular function type and set type.

```julia-repl
julia> all_constraints(model, GeneralVariableRef, MOI.LessThan)
1-element Array{InfOptConstraintRef,1}:
 x ≤ 1.0
```
"""
function JuMP.all_constraints(model::InfiniteModel,
                              function_type::Type{<:JuMP.AbstractJuMPScalar},
                              set_type::Type{<:MOI.AbstractSet}
                              )::Vector{InfOptConstraintRef}
    constr_list = Vector{InfOptConstraintRef}(undef,
                           JuMP.num_constraints(model, function_type, set_type))
    counter = 1
    for (index, object) in model.constraints
        if isa(JuMP.jump_function(object.constraint), function_type) &&
           isa(JuMP.moi_set(object.constraint), set_type)
            constr_list[counter] = _make_constraint_ref(model, index)
            counter += 1
        end
    end
    return constr_list
end

"""
    JuMP.all_constraints(model::InfiniteModel,
                         function_type::Type{<:JuMP.AbstractJuMPScalar}
                         )::Vector{InfOptConstraintRef}

Extend [`JuMP.all_constraints`](@ref JuMP.all_constraints(::JuMP.Model, ::Type{<:Union{JuMP.AbstractJuMPScalar, Vector{<:JuMP.AbstractJuMPScalar}}}, ::Type{<:MOI.AbstractSet}))
to search by function types for all MOI
sets and return a list of all constraints use a particular function type.

```julia-repl
julia> all_constraints(model, GeneralVariableRef)
3-element Array{InfOptConstraintRef,1}:
 x ≥ 0.0
 x ≤ 3.0
 x integer
```
"""
function JuMP.all_constraints(model::InfiniteModel,
                              function_type::Type{<:JuMP.AbstractJuMPScalar}
                              )::Vector{InfOptConstraintRef}
    return JuMP.all_constraints(model, function_type, MOI.AbstractSet)
end

"""
    JuMP.all_constraints(model::InfiniteModel,
                         set_type::Type{<:MOI.AbstractSet}
                         )::Vector{InfOptConstraintRef}

Extend [`JuMP.all_constraints`](@ref JuMP.all_constraints(::JuMP.Model, ::Type{<:Union{JuMP.AbstractJuMPScalar, Vector{<:JuMP.AbstractJuMPScalar}}}, ::Type{<:MOI.AbstractSet}))
to search by MOI set type for all function
types and return a list of all constraints that use a particular set type.

```julia-repl
julia> all_constraints(model, MOI.GreaterThan)
3-element Array{InfOptConstraintRef,1}:
 x ≥ 0.0
 g(t) ≥ 0.0, ∀ t ∈ [0, 6]
 g(0.5) ≥ 0.0
```
"""
function JuMP.all_constraints(model::InfiniteModel,
                              set_type::Type{<:MOI.AbstractSet}
                              )::Vector{InfOptConstraintRef}
    return JuMP.all_constraints(model, JuMP.AbstractJuMPScalar, set_type)
end

"""
    JuMP.all_constraints(model::InfiniteModel)::Vector{InfOptConstraintRef}

Extend [`JuMP.all_constraints`](@ref JuMP.all_constraints(::JuMP.Model, ::Type{<:Union{JuMP.AbstractJuMPScalar, Vector{<:JuMP.AbstractJuMPScalar}}}, ::Type{<:MOI.AbstractSet}))
to return all a list of all the constraints
in an infinite model `model`.

```julia-repl
julia> all_constraints(model)
5-element Array{InfOptConstraintRef,1}:
 x ≥ 0.0
 x ≤ 3.0
 x integer
 g(t) ≥ 0.0, ∀ t ∈ [0, 6]
 g(0.5) ≥ 0.0
```
"""
function JuMP.all_constraints(model::InfiniteModel)::Vector{InfOptConstraintRef}
    return [_make_constraint_ref(model, index)
            for (index, object) in model.constraints]
end

"""
    JuMP.list_of_constraint_types(model::InfiniteModel)::Vector{Tuple{DataType, DataType}}

Extend [`JuMP.list_of_constraint_types`](@ref JuMP.list_of_constraint_types(::JuMP.Model))
to return a list of tuples that
contain all the used combinations of function types and set types in the model.

```julia-repl
julia> all_constraints(model)
3-element Array{Tuple{DataType,DataType},1}:
 (GeneralVariableRef, MathOptInterface.LessThan{Float64})
 (GeneralVariableRef, MathOptInterface.GreaterThan{Float64})
 (GeneralVariableRef, MathOptInterface.Integer)
```
"""
function JuMP.list_of_constraint_types(model::InfiniteModel
    )::Vector{Tuple{DataType, DataType}}
    type_set = Set{Tuple{DataType, DataType}}()
    for (index, object) in model.constraints
        push!(type_set, (typeof(JuMP.jump_function(object.constraint)),
                          typeof(JuMP.moi_set(object.constraint))))
    end
    return collect(type_set)
end

################################################################################
#                         PARAMETER REFERENCE METHODS
################################################################################
"""
    parameter_refs(cref::InfOptConstraintRef)::Tuple

Return the tuple of infinite parameter references that determine the infinite
dependencies of `cref`.

**Example**
```julia-repl
julia> parameter_refs(cref)
(t,)
```
"""
function parameter_refs(cref::InfOptConstraintRef)::Tuple
    model = JuMP.owner_model(cref)
    obj_indices = _param_object_indices(model)[_object_numbers(cref)]
    return Tuple(_make_param_tuple_element(model, idx) for idx in obj_indices)
end

################################################################################
#                           PARAMETER BOUND METHODS
################################################################################
"""
    has_parameter_bounds(cref::InfOptConstraintRef)::Bool

Return a `Bool` indicating if `cref` is limited to a sub-domain as defined
by a [`ParameterBounds`](@ref) object.

**Example**
```julia-repl
julia> has_parameter_bounds(cref)
true
```
"""
function has_parameter_bounds(cref::InfOptConstraintRef)::Bool
    if JuMP.constraint_object(cref) isa BoundedScalarConstraint
        return !isempty(JuMP.constraint_object(cref).bounds)
    else
        return false
    end
end

"""
    parameter_bounds(cref::InfOptConstraintRef)::ParameterBounds{GeneralVariableRef}

Return the [`ParameterBounds`](@ref) object associated with the constraint
`cref`. Errors if `cref` does not have parameter bounds.

**Example**
```julia-repl
julia> parameter_bounds(cref)
Subdomain bounds (1): t ∈ [0, 2]
```
"""
function parameter_bounds(cref::InfOptConstraintRef)::ParameterBounds{GeneralVariableRef}
    !has_parameter_bounds(cref) && error("$cref does not have parameter bounds.")
    return JuMP.constraint_object(cref).bounds
end

# Internal function used to change the parameter bounds of a constraint
function _update_constr_param_bounds(cref::InfOptConstraintRef,
                                     bounds::ParameterBounds{GeneralVariableRef},
                                     orig_bounds::ParameterBounds{GeneralVariableRef}
                                     )::Nothing
    old_constr = JuMP.constraint_object(cref)
    func = JuMP.jump_function(old_constr)
    set = JuMP.moi_set(old_constr)
    if !isempty(bounds)
        new_constr = BoundedScalarConstraint(func, set, bounds, orig_bounds)
    else
        new_constr = JuMP.ScalarConstraint(func, set)
    end
    _set_core_constraint_object(cref, new_constr)
    return
end

"""
    set_parameter_bounds(cref::InfOptConstraintRef,
                         bounds:ParameterBounds{GeneralVariableRef};
                         [force::Bool = false])::Nothing

Specify a new [`ParameterBounds`](@ref) object `bounds` for the constraint `cref`.
This is meant to be primarily used by [`@set_parameter_bounds`](@ref) which
provides a more intuitive syntax.

**Example**
```julia-repl
julia> set_parameter_bounds(cref, ParameterBounds((t => IntervalSet(0, 2),)))

julia> parameter_bounds(cref)
Subdomain bounds (1): t ∈ [0, 2]
```
"""
function set_parameter_bounds(cref::InfOptConstraintRef,
                              bounds::ParameterBounds{GeneralVariableRef};
                              force::Bool = false, _error::Function = error
                              )::Nothing
    if has_parameter_bounds(cref) && !force
        _error("$cref already has parameter bounds. Consider adding more using " *
               "`add_parameter_bounds` or overwriting them by setting " *
               "the keyword argument `force = true`")
    else
        # check that bounds are valid and add support(s) if necessary
        _check_bounds(bounds, _error = _error)
        _validate_bounds(JuMP.owner_model(cref), bounds, _error = _error)
        # consider hold variables
        orig_bounds = copy(bounds)
        vrefs = _all_function_variables(JuMP.jump_function(JuMP.constraint_object(cref)))
        for vref in vrefs
            _update_var_bounds(vref, bounds)
        end
        # set the new bounds
        _update_constr_param_bounds(cref, bounds, orig_bounds)
        # update status
        set_optimizer_model_ready(JuMP.owner_model(cref), false)
    end
    return
end

"""
    add_parameter_bounds(cref::InfOptConstraintRef,
                         new_bounds::ParameterBounds{GeneralVariableRef}
                         )::Nothing

Add an additional parameter bound to `cref` such that it is defined over the
sub-domain based on `pref` from `lower` to `upper`. This is primarily meant to be
used by [`@add_parameter_bounds`](@ref).

```julia-repl
julia> add_parameter_bounds(cref, ParameterBounds((t => IntervalSet(0, 2),))

julia> parameter_bounds(cref)
Subdomain bounds (1): t ∈ [0, 2]
```
"""
function add_parameter_bounds(cref::InfOptConstraintRef,
                              new_bounds::ParameterBounds{GeneralVariableRef};
                              _error::Function = error)::Nothing
    # check the new bounds
    _check_bounds(new_bounds, _error = _error)
    _validate_bounds(JuMP.owner_model(cref), new_bounds, _error = _error)
    # add the bounds
    if JuMP.constraint_object(cref) isa BoundedScalarConstraint
        _update_bounds(parameter_bounds(cref), new_bounds, _error = _error)
        _update_bounds(original_parameter_bounds(JuMP.constraint_object(cref)),
                       new_bounds, _error = _error)
    else
        _update_constr_param_bounds(cref, new_bounds, copy(new_bounds))
    end
    # update the optimizer model status
    set_optimizer_model_ready(JuMP.owner_model(cref), false)
    return
end

"""
    delete_parameter_bounds(cref::InfOptConstraintRef)::Nothing

Delete all the parameter bounds of the constraint `cref`. Note any bounds that
are needed for hold variables inside in `cref` will be unaffected.

**Example**
```julia-repl
julia> @BDconstraint(model, c1(x == 0), y <= 42)
c1 : y(x) ≤ 42, ∀ x[1] = 0, x[2] = 0

julia> delete_parameter_bounds(c1)

julia> c1
c1 : y(x) ≤ 42, ∀ x[1] ∈ [-1, 1], x[2] ∈ [-1, 1]
```
"""
function delete_parameter_bounds(cref::InfOptConstraintRef)::Nothing
    # check if has bounds on pref from the constraint
    constr = JuMP.constraint_object(cref)
    if has_parameter_bounds(cref) && !isempty(original_parameter_bounds(constr))
        # consider hold variables
        new_bounds = ParameterBounds()
        vrefs = _all_function_variables(JuMP.jump_function(constr))
        for vref in vrefs
            _update_var_bounds(vref, new_bounds)
        end
        # set the new bounds
        _update_constr_param_bounds(cref, new_bounds, ParameterBounds())
    end
    return
end

################################################################################
#                                 DELETION
################################################################################
"""
    JuMP.delete(model::InfiniteModel, cref::InfOptConstraintRef)::Nothing

Extend [`JuMP.delete`](@ref JuMP.delete(::JuMP.Model, ::JuMP.ConstraintRef{JuMP.Model}))
to delete an `InfiniteOpt` constraint and all associated information. Errors
if `cref` is invalid.

**Example**
```julia-repl
julia> print(model)
Min measure(g(t)*t) + z
Subject to
 z ≥ 0.0
 g(t) + z ≥ 42.0, ∀ t ∈ [0, 6]

julia> delete(model, cref)

julia> print(model)
Min measure(g(t)*t) + z
Subject to
 z ≥ 0.0
```
"""
function JuMP.delete(model::InfiniteModel, cref::InfOptConstraintRef)::Nothing
    # check valid reference
    @assert JuMP.is_valid(model, cref) "Invalid constraint reference."
    # update variable dependencies
    constr = _core_constraint_object(cref)
    all_vrefs = _all_function_variables(JuMP.jump_function(constr))
    for vref in all_vrefs
        filter!(e -> e != JuMP.index(cref), _constraint_dependencies(vref))
    end
    # delete constraint information
    _delete_data_object(cref)
    # reset optimizer model status
    set_optimizer_model_ready(model, false)
    return
end

ext: jl


text:
@testset "Unconstrained" begin
    include("double_msd.jl")
    include("AJPR14e54.jl")
    include("AP12e21.jl")
    include("JCG14e61.jl")
    include("JCG14e63.jl")
    include("PJ08e28.jl")
    include("PJ08e54.jl")
end
@testset "Constrained" begin
    include("PEDJ16s4.jl")
end


text:
module VLM

struct CPdata
    chord::Array{Float64, 1}
    twist::Array{Float64, 1}
    tc::Array{Float64, 1}
    x::Array{Float64, 1}
    y::Array{Float64, 1}
    z::Array{Float64, 1}
    dihedral::Array{Float64, 1}
    sweep::Array{Float64, 1}
    ds::Array{Float64, 1}
end

struct position
    x::Array{Float64, 1}
    y::Array{Float64, 1}
    z::Array{Float64, 1}
end

struct wingsection
    span
    chord
    twist
    tc
    sweep
    dihedral
    N
end

struct fs_def
    alpha
    beta
    rho
    Vinf
    # CL
    # method
end
struct ref_def
    S
    c
    CLmax
end

struct velocity
    x
    y
    z
end
# struct pdrag_def
#     polar
#     alt
#     xt
#     method
# end
# struct mvr_def
#     qN
#     n
#     kbar
# end

## -------- geometry ---------------

function geometry(wing::wingsection)

    # ----------------- rename for convenience --------------------
    b = wing.span
    Lambda = wing.sweep
    chord = wing.chord
    twist = wing.twist
    tc = wing.tc
    phi = wing.dihedral
    N = sum(wing.N)
    P = wing.N

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

    # P = round.(Int64, b/sum(b)*N) # divide up panels

    # -------------  Quarter Chord Locations --------------------------
    M = 1 + sum(P)

    c = zeros(M)
    t = zeros(M)
    thickness = zeros(M)
    QCx = zeros(M)
    QCy = zeros(M)
    QCz = zeros(M)
    LEx = zeros(M)
    LEy = zeros(M)
    LEz = zeros(M)
    TEx = zeros(M)
    TEy = zeros(M)
    TEz = zeros(M)
    CPx = zeros(N)
    CPy = zeros(N)
    CPz = zeros(N)

    last = 1

    for i = 1:length(b)
        first = last
        last = first + P[i]
        eta = linspace(0, b[i], P[i]+1)

        QCx[first:last] = QCx[first] + eta*tan(Lambda[i])
        QCy[first:last] = QCy[first] + eta*cos(phi[i])
        QCz[first:last] = QCz[first] + eta*sin(phi[i])
        c[first:last] = chord[i] + eta*(chord[i+1]-chord[i])/b[i] # chord
        t[first:last] = twist[i] + eta*(twist[i+1]-twist[i])/b[i] # twist
        thickness[first:last] = tc[i]*chord[i] + eta*(tc[i+1]*chord[i+1]-tc[i]*chord[i])/b[i] # thickness
    end
    QC = position(QCx, QCy, QCz)
    # --------------------------------------------------------------

    # ------------ Trailing Edge Locations --------------------------
    TEx = QCx + 3.0/4*c
    TEy = QCy
    TEz = QCz
    TE = position(TEx, TEy, TEz)
    # ----------------------------------------------------------------

    # -------------- Leading Edge locations ------------------------
    LEx = QCx - c/4
    LEy = QCy
    LEz = QCz
    LE = position(LEx, LEy, LEz)
    # -----------------------------------------------------------------

    # ------------- Control Point Locations --------------------------
    chord = 0.5*(c[1:N] + c[2:N+1])
    twist = 0.5*(t[1:N] + t[2:N+1])
    tc = 0.5*(thickness[1:N] + thickness[2:N+1])./chord
    CPx = 0.5*(QC.x[1:N] + QC.x[2:N+1]) + 0.5*chord
    CPy = 0.5*(QC.y[1:N] + QC.y[2:N+1])
    CPz = 0.5*(QC.z[1:N] + QC.z[2:N+1])

    dihedral = zeros(N)
    sweep = zeros(N)

    last = 0
    for i = 1:length(b)
      first = last + 1
      last = first + P[i] - 1
      dihedral[first:last] = phi[i]
      sweep[first:last] = Lambda[i]
    end
    ds = sqrt.((QC.y[2:N+1]-QC.y[1:N]).^2 + (QC.z[2:N+1]-QC.z[1:N]).^2)

    CP = CPdata(chord, twist, tc, CPx, CPy, CPz, dihedral, sweep, ds)

    return QC, TE, CP, LE

end

## -------- influence coefficients ---------------

#=
--------------------
lift influence coefficients

Author: S. Andrew Ning
---------------------
=#
function getLIC(CP::CPdata, rho::Float64, Vinf::Array{Float64, 1})
    LIC = 2*rho*Vinf.*cos.(CP.dihedral).*CP.ds
end

#=
-----------------------------------
Pitching moment coefficients about (x) c.g. location

Author: S. Andrew Ning
Updates: 1/30/09 - removed dependence on QC since this is known
-------------------------------------
=#
function getMIC(CP::CPdata, rho::Float64, Vinf::Float64, xcg::Float64)

    xQC = CP.x - 0.5*CP.chord # quarter chord locations

    MIC = -(xQC - xcg)*2*rho*Vinf.*cos.(CP.dihedral).*CP.ds

end

#=
------------------------------------------------
Induced drag influence coefficients with a Trefftz plane calculation
Di = gamma'*DIC*gamma

Author: S. Andrew Ning

Updated feb. 13 2008
1 is control points, 2 is vorticies.  If you leave out last two
arguments the control points and vorticies are on one wing.  Include them
if you need effect of wing on tail for example.
-------------------------------------------------
=#

function getDIC(yFF::Array{Float64, 1}, zFF::Array{Float64, 1}, rho::Float64)
    getDIC(yFF, zFF, rho, yFF, zFF)
end

function getDIC(yFF::Array{Float64, 1}, zFF::Array{Float64, 1}, rho::Float64, yFF2::Array{Float64, 1}, zFF2::Array{Float64, 1})

    n = length(yFF)-1
    n2 = length(yFF2)-1

     # ------------ define useful variables ------------------
    yc = 0.5*(yFF[2:n+1] + yFF[1:n])   # center of panels
    zc = 0.5*(zFF[2:n+1] + zFF[1:n])
    ny = zFF[2:n+1] - zFF[1:n]
    nz = yFF[2:n+1] - yFF[1:n]
     # ----------------------------------------------------

    DIC = zeros(n,n2)
     #------- normal wash calculation ----------------
    for i = 1:n
        for j = 1:n2
            ry = yFF2[j] - yc[i]
            rz = zFF2[j] - zc[i]
            r2 = ry^2 + rz^2
            DIC[i, j] = -rho/(2*pi*r2)*(rz*ny[i]+ry*nz[i])

             # add other side of panel
            ry = yFF2[j+1] - yc[i]
            rz = zFF2[j+1] - zc[i]
            r2 = ry^2 + rz^2
            DIC[i, j] += rho/(2*pi*r2)*(rz*ny[i]+ry*nz[i])
        end
    end

    yFF2 = -yFF2 # add left side

    for i = 1:n
        for j = 1:n2
            ry = yFF2[j] - yc[i]
            rz = zFF2[j] - zc[i]
            r2 = ry^2 + rz^2
            DIC[i,j] += rho/(2*pi*r2)*(rz*ny[i]+ry*nz[i])

            # add other side of panel
            ry = yFF2[j+1] - yc[i]
            rz = zFF2[j+1] - zc[i]
            r2 = ry^2 + rz^2
            DIC[i,j] -= rho/(2*pi*r2)*(rz*ny[i]+ry*nz[i])
        end
    end

    return DIC
end

#=
---------------------
assumes a parabolic variation in viscous drag with section lift coefficient
Dp = D0 + D1*gamma + D2*gamma.^2

Author: S. Andrew Ning
-----------------------
=#

function getViscousDrag(cd1::Float64, cd2::Float64, CP::CPdata, rho::Float64, Vinf::Float64)

  D1 = cd1*rho*Vinf*CP.ds
  D2 = cd2*2*rho./CP.chord.*CP.ds

  return D1, D2
end

#=
----------------------------------------
 Calculates velocity induced by vortex (with unit vorticity)
 starting at rA ending at rB at control point rC. (see Bertin/Smith)

 Author: S. Andrew Ning
------------------------------------------
=#

function vortex(xA::Array{Float64, 1}, yA::Array{Float64, 1}, zA::Array{Float64, 1},
        xB::Array{Float64, 1}, yB::Array{Float64, 1}, zB::Array{Float64, 1},
        xC::Array{Float64, 1}, yC::Array{Float64, 1}, zC::Array{Float64, 1})

    m = length(yC)
    n = length(yA)
    v = zeros(m, n)
    w = zeros(m, n)

    for i = 1:m
        for j = 1:n
            # define repeatedly used values
            x1 = xC[i] - xA[j]
            x2 = xC[i] - xB[j]
            x21 = xB[j] - xA[j]
            y1 = yC[i] - yA[j]
            y2 = yC[i] - yB[j]
            y21 = yB[j] - yA[j]
            z1 = zC[i] - zA[j]
            z2 = zC[i] - zB[j]
            z21 = zB[j] - zA[j]

            denom1 = (y1*z2-y2*z1)^2 + (x1*z2-x2*z1)^2 + (x1*y2-x2*y1)^2
            # i1 = y1*z2 - y2*z1
            j1 = -x1*z2 + x2*z1
            k1 = x1*y2 - x2*y1
            frac1_1 = (x21*x1 + y21*y1 + z21*z1)/sqrt(x1^2 + y1^2 + z1^2)
            frac1_2 = (x21*x2 + y21*y2 + z21*z2)/sqrt(x2^2 + y2^2 + z2^2)
            frac1 = frac1_1 - frac1_2
            # u[i, j] = 1.0/(4*pi)*i1/denom1*frac1
            v[i, j] = 1.0/(4*pi)*j1/denom1*frac1
            w[i, j] = 1.0/(4*pi)*k1/denom1*frac1
        end
    end

    return v, w
end

#=
------------------------------------------------
This function calculates the induced velocies at control points
xCP,yCP,zCP from the vorticies which are bound at the quarter chord,
follow the wing to the trailing edge, and leave the wing parallel to the
freestream (drag-free wake).  All values are computed for unit vorticity

Author: Andrew Ning
------------------------------------------------
=#
function velocity3d(QC::position, TE::position, CP::CPdata)

    # first calculate induced velocities from vorticies on right side of wing.
    m = length(CP.y)
    n = length(QC.y)-1

    # u = zeros(m, n)
    v = zeros(m, n)
    w = zeros(m, n)

    # Induced velocities due to bound vortex (BC)
    dv, dw = vortex(QC.x[1:n], QC.y[1:n], QC.z[1:n], QC.x[2:n+1], QC.y[2:n+1], QC.z[2:n+1],
        CP.x, CP.y, CP.z)
    # u += du
    v += dv
    w += dw

    # Induced velocities due to AB - left vortex bound to wing
    dv, dw = vortex(TE.x[1:n], TE.y[1:n], TE.z[1:n], QC.x[1:n], QC.y[1:n], QC.z[1:n],
        CP.x, CP.y, CP.z)
    # u += du
    v += dv
    w += dw

    # Induced velocities due to CD - right vortex bound to wing
    dv, dw = vortex(QC.x[2:n+1], QC.y[2:n+1], QC.z[2:n+1], TE.x[2:n+1], TE.y[2:n+1], TE.z[2:n+1],
        CP.x, CP.y, CP.z)
    # u += du
    v += dv
    w += dw

    # Induced velocities due to trailing vortex from A
    l = 1e9  # some large number to represent infinity
    # xINF = TE.x + l*cos(alpha)*ones(size(TE.x))
    # yINF = TE.y
    # zINF = TE.z + l*sin(alpha)*ones(size(TE.z))
    xINF = TE.x + l
    yINF = TE.y
    zINF = TE.z

    dv, dw = vortex(xINF[1:n], yINF[1:n], zINF[1:n], TE.x[1:n], TE.y[1:n], TE.z[1:n],
        CP.x, CP.y, CP.z)
    # u += du
    v += dv
    w += dw

    # Induced velocities due to trailing vortex from D
    dv, dw = vortex(TE.x[2:n+1], TE.y[2:n+1], TE.z[2:n+1], xINF[2:n+1], yINF[2:n+1], zINF[2:n+1],
        CP.x, CP.y, CP.z)
    # u += du
    v += dv
    w += dw

    # ---------- Now repeat with contribution from left side of the wing -----
    # y switched signs (except control points of course)
    # QCy = -QC.y
    # TEy = -TE.y
    # yINF = -yINF
    # Also all induced velocities have opposite sign because vorticies travel
    # in opposite direction.

    # Induced velocities due to bound vortex (BC)
    dv, dw = vortex(QC.x[1:n], -QC.y[1:n], QC.z[1:n], QC.x[2:n+1], -QC.y[2:n+1], QC.z[2:n+1],
        CP.x, CP.y, CP.z)
    # u -= du
    v -= dv
    w -= dw

    # Induced velocities due to AB - left vortex bound to wing
    dv, dw = vortex(TE.x[1:n], -TE.y[1:n], TE.z[1:n], QC.x[1:n], -QC.y[1:n], QC.z[1:n],
        CP.x, CP.y, CP.z)
    # u -= du
    v -= dv
    w -= dw

    # Induced velocities due to CD - right vortex bound to wing
    dv, dw = vortex(QC.x[2:n+1], -QC.y[2:n+1], QC.z[2:n+1], TE.x[2:n+1], -TE.y[2:n+1], TE.z[2:n+1],
        CP.x, CP.y, CP.z)
    # u -= du
    v -= dv
    w -= dw

    # Induced velocities due to trailing vortex from A

    dv, dw = vortex(xINF[1:n], -yINF[1:n], zINF[1:n], TE.x[1:n], -TE.y[1:n], TE.z[1:n],
        CP.x, CP.y, CP.z)
    # u -= du
    v -= dv
    w -= dw

    # Induced velocities due to trailing vortex from D
    dv, dw = vortex(TE.x[2:n+1], -TE.y[2:n+1], TE.z[2:n+1], xINF[2:n+1], -yINF[2:n+1], zINF[2:n+1],
        CP.x, CP.y, CP.z)
    # u -= du
    v -= dv
    w -= dw

    return v, w
end

#=
------------------------------------
Aerodynamic influence coefficients
Vn = AIC*gamma

Author: S. Andrew Ning
Updates: summer 08 - simplified so that twist/camber are only applied in boundary conditions
                     consistent with assumptions of VLM.
--------------------------------------
=#

function getAIC(QC::position, TE::position, CP::CPdata)

    v, w =  velocity3d(QC, TE, CP)

    m = length(CP.y)
    n = length(QC.y)-1

    AIC = zeros(m, n)

    for j = 1:n
        AIC[:, j] = -v[:, j].*sin.(CP.dihedral) + w[:, j].*cos.(CP.dihedral)
    end

    return AIC
end

#=
-------------------------------------------------
this function computes integrated bending moment over thickness and the
bending moment distribution
W = WIC*gamma;
bending moment distribution = BMM*gamma

Author: S. Andrew Ning
Updates: 1/23/09 - corrected integration along spar and vectorized for speed
         1/30/09 - allows for thickness distribution
---------------------------------------------------
=#

function getWIC(CP::CPdata, rho::Float64, Vinf::Float64)

    x = CP.x - 0.5*CP.chord
    y = CP.y
    z = CP.z

    # rename for convenience
    N = length(y)
    phi = CP.dihedral
    Lambda = CP.sweep
    chord = CP.chord
    tc = CP.tc
    ds = CP.ds

    R = zeros(N, N)

    for i = 1:N
        for j = i+1:N
            cij = cos(CP.dihedral[j])*cos(CP.dihedral[i]) + sin(CP.dihedral[j])*sin(CP.dihedral[i])
            R[i, j] = (x[j] - x[i])*cij*sin(CP.sweep[i]) +
                     (y[j] - y[i])*cos(CP.dihedral[j])*cos(CP.sweep[i]) +
                     (z[j] - z[i])*sin(CP.dihedral[j])*cos(CP.sweep[i])

            # if j <= i, R(i,j) = 0
        end
    end

    # integrate along structural span
    ds_str = CP.ds./cos.(CP.sweep)

    BMM = rho*Vinf*R.*repmat(ds.', N, 1)  #(ones(N, 1)*ds)
    WIC = (2*(ds_str./(tc.*chord)).'*BMM).'

    return WIC, BMM
end

# TODO: haven't tested which is faster in Julia

# # compute R matrix
# PhiM = ((cos(phi)*cos(phi)') + (sin(phi)*sin(phi)'))
# Rx = (sin(Lambda)*x' - (x.*sin(Lambda))*ones(1,N)).*PhiM
# Ry = cos(Lambda)*(y.*cos(phi))' - (y.*cos(Lambda))*cos(phi)'
# Rz = cos(Lambda)*(z.*sin(phi))' - (z.*cos(Lambda))*sin(phi)'
# R = Rx + Ry + Rz
# R = R - tril(R) # R(i,j) = 0 for j <= i

# % below version is exactly the same (slower since its in loops, but left
# % here since it's easier to read then vectorized version above)
#
# % R = zeros(N);
#
# % for i = 1:N
# %     for j = 1:N
# %         cij = cos(CP.dihedral(j))*cos(CP.dihedral(i)) + sin(CP.dihedral(j))*sin(CP.dihedral(i));
# %         R(i,j) = (x(j) - x(i))*cij*sin(CP.sweep(i)) + ...
# %                  (y(j) - y(i))*cos(CP.dihedral(j))*cos(CP.sweep(i)) + ...
# %                  (z(j) - z(i))*sin(CP.dihedral(j))*cos(CP.sweep(i));
# %         if j <= i
# %             R(i,j) = 0;
# %         end
# %     end
# % end

# % ds_str = CP.ds./cos(CP.sweep);
#
# % WIC = zeros(1,N);
# %
# % for j = 1:N
# %     WIC(j) = 2*rho*U/tc*CP.ds(j).*ds_str./CP.chord*R(:,j);
# % end

#=
-----------------------------------------------------
This function computes angle of attack necessary to match a specified
lift

Author: S. Andrew Ning
Updates: 1/28/09 - created
         9/2/09 - changed the way alpha was computed - seems more robust
                 to larger angles of attack by using arcsin rather than arccos
---------------------------------------------------
=#

# function getAlpha(Ltarget::Float64, CP::CPdata, LIC::Array{Float64, 1}, AIC::Array{Float64, 2}, Vinf::Float64)
#
#     # solve quadratic equation for alpha
#     qq = LIC/AIC
#     a = -U*qq*sin.(CP.twist)'
#     b = -U*qq*(cos.(CP.twist).*cos.(CP.dihedral))'
#     c = Ltarget
#     # alpha = acos((a*c + b*sqrt(a^2+b^2-c^2))/(a^2+b^2));
#     alpha = asin(c/b - a/b*(a*c + b*sqrt(a^2+b^2-c^2))/(a^2 + b^2))
#
#     return alpha
# end

## -------- vortex lattice method ---------------

#=
----------------------------------------------------------------------
This version assumes (one) symmetric wing.

Author: S. Andrew Ning
Updates: repackaged with features from several different versions - 1/30/09
        2/2/09 - fixed computation of t/c_avg to be properly area weighted
                (used in compressibility and Re dependent parasite drag)
               - clmax is now a function of t/c (polynomial fit from 241 data)
        2/3/09 - added viscous dependent induced drag if PASS method is used to compute
                 parasite drag
        2/17/09 - added area dependent weight, correct area in quadratic parasite
                  drag calculation, parameter to choose pdrag method
----------------------------------------------------------------------
=#

function run(wing, fs, Vext)

    # freestream properties
    # if fs.method == "CL" # specify CL
    #     CLref = fs.CL
    #     alpha = 0.0
    # else
    #     alpha = fs.alpha
    #     CLref = 0.0
    # end
    alpha = fs.alpha
    beta = fs.beta
    rho = fs.rho
    Vinf = fs.Vinf

    # # mach = fs.mach
    # rho = 1.0
    # Vinf = 1.0
    # q = 0.5*rho*Vinf^2

    # reference quantities
    # Sref = ref.S
    # cref = ref.c
    # CLmax = ref.CLmax

    # # viscous drag (2 possible methods)
    # if pdrag.method == "pass" # Reynolds number dependent method
    #     alt = pdrag.alt
    #     xt = pdrag.xt
    #
    #     # not used
    #     cd0 = 0.0
    #     cd1 = 0.0
    #     cd2 = 0.0
    #
    # else # quadratic variation with section lift coefficient
    #     cd0 = pdrag.polar[1]
    #     cd1 = pdrag.polar[2]
    #     cd2 = pdrag.polar[3]
    #
    #     # not used
    #     alt = 0.0
    #     xt = 0.0
    # end

    # # structures
    # qmvrN = mvr.qN # ratio of maneuver dynamic pressure to cruise dynamic pressure
    # n = mvr.n # load factor
    # kbar = mvr.kbar # coefficient used in computing area dependent weight (3.57e4 or 0)

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

    # ------------- geometry ------------------
    QC, TE, CP, LE = geometry(wing)

    # compute wing area
    # S = 2*sum(CP.chord.*CP.ds)
    # ------------------------------------------

    # ---------- force influence coefficients --------------------
    # aerodynamic influence coefficients
    AIC = getAIC(QC, TE, CP)

    # lift
    # LIC = getLIC(CP, rho, Vinfeff)

    # induced drag
    DIC = getDIC(TE.y, TE.z, rho)

    # viscous drag
    # D1, D2 = getViscousDrag(cd1, cd2, CP, rho, Vinf)

    # weight
    # WIC, BMM = getWIC(CP, rho, Vinf)
    # -----------------------------------------------

    # # ---------- compute angle of attack necessary to match CL -----
    # if fs.method == "CL"
    #     alpha = getAlpha(CLref*q*Sref, CP, LIC, AIC, Vinf)
    # end
    # # ------------------------------------------------------------------

    # -------- compute circulation -----------------
    Vinfn = (Vinf*cos(alpha)*cos(beta)*sin.(CP.twist) +
        Vinf*sin(beta)*cos.(CP.twist).*sin.(CP.dihedral) +
        Vinf*sin(alpha)*cos(beta)*cos.(CP.twist).*cos.(CP.dihedral))
    Vextn = (Vext.x.*cos(beta).*sin.(CP.twist) +
        -Vext.y.*cos.(CP.twist).*sin.(CP.dihedral) +
        Vext.z.*cos.(CP.twist).*cos.(CP.dihedral))
    Vn = -(Vinfn + Vextn)

    gamma = AIC\Vn
    # ----------------------------------------------

    # total velocity in direction of Vinf
    Vinfeff = Vinf + Vext.x*cos(alpha)*cos(beta) + Vext.y*sin(beta) + Vext.z*sin(alpha)*cos(beta)

    # effective angle of attack
    alphaeff = atan.(-Vn/Vinf)

    # --------- aerodynamic forces ------------------
    # L = dot(LIC, gamma)
    # Kutta Joukowski
    L = sum(2*rho*Vinfeff.*gamma.*cos.(CP.dihedral).*CP.ds)

    Di = gamma'*DIC*gamma

    # add near field drag from external normalwash
    Di += sum(-rho*Vextn.*gamma.*CP.ds)

    # Dp = cd0*q*S + D1'*gamma + D2'*gamma.^2

    # CL = L/q/Sref
    # CDi = Di/q/Sref
    # CDp = Dp/q/Sref
    # ------------------------------------------------

    # # --------- weight (integrated bending moment over thickness --------
    # circulation at maneuver load
    # bbb = AIC\cos.(CP.dihedral)
    # LL = dot(LIC, bbb)
    # gamma_mvr = gamma + ((n/qmvrN-1)*(LIC'*gamma)/LL*bbb')'
    #
    # # compute weight
    # W = qmvrN*WIC'*gamma_mvr
    #
    # # add area dependent weight
    # W = W + kbar*S
    #
    # # weight coefficient
    # CW = W/q/Sref/cref
    # # -----------------------------------------------------------------

    # # ----- estimate compressibility and parasite drag (strip theory + PASS method) -----
    # supercrit = 1
    # P = round(wing.span/sum(wing.span)*wing.N) # number of panels in each section
    # cdc = zeros(1,length(wing.span))
    # cdp = zeros(1,length(wing.span))
    # area = zeros(1,length(wing.span))
    # start = 1
    #
    # for i = 1:length(wing.span)
    #   finish = start + P[i] - 1
    #
    #   # rename for convenience
    #   cr = wing.chord[i]
    #   ct = wing.chord[i+1]
    #   cbar = 1/2*(cr + ct)
    #   tcbar = (wing.tc[i]*wing.chord[i] + wing.tc[i+1]*wing.chord[i+1])/(wing.chord[i]+wing.chord[i+1])
    #   area[i] = cbar*wing.span[i]
    #   mac = 2/3*(cr + ct - cr*ct/(cr+ct))
    #   CL_local = 0.1*sum(gamma[start:finish]'.*CP.ds[start:finish])*2/U/area[i]
    #
    #   # compressibility drag
    #   cdc[i] = Cdrag(CL_local,wing.sweep[i],tcbar,mach,supercrit)
    #
    #   # parasite drag
    #   cdp[i] = Pdrag(alt,mach,xt,mac,wing.sweep[i],tcbar)
    #
    #   start = finish + 1
    # end
    #
    # # compressibility drag - area weighted average
    # CDc = 2*sum(cdc.*area)/Sref[1]
    #
    # if pass_mthd
    #   # parasite drag - area weighted average
    #   CDp = 2*sum(cdp.*area)/Sref[1]
    #
    #   # add viscous dependent induced drag
    #   Lambda_bar = area*wing.sweep'/sum(area)
    #   CDi = CDi + 0.38*CDp*CL^2/cos(Lambda_bar)^2
    # end
    # # ----------------------------------------------------------------------

    # ----------- cl distribution at CLmax ----------------
    cl = 2/Vinf*gamma./CP.chord
    cllocal = cl*Vinf./Vinfeff
    # clmax_dist = cl + rho*Vinf*(CLmax*Sref-L/q)/LL*bbb./CP.chord
    #
    # # clmax as a function of thickness - polynomial fit
    # tc = CP.tc*100
    # clmax = -1.748 + 0.8013*tc - 0.06567*tc.^2 + 0.0022307*tc.^3 - 2.7634e-5*tc.^4
    #
    # cl_margin = clmax - clmax_dist
    # ----------------------------------------

    # # --------- pitching moment about a.c. --------------------
    # # pitching moment about quarter chord
    # MIC = getMIC(CP, rho, Vinf, QC.x[1])
    #
    # # find aerodynamic center
    # dRHS = -sin(alpha)*sin.(CP.twist) + cos(alpha)*cos.(CP.twist).*cos.(CP.dihedral)
    # dbc = AIC\dRHS
    # dMda = MIC'*dbc
    # dLda = LIC'*dbc
    # xac = -dMda/dLda + QC.x[1]
    #
    # # find Moment about a.c.
    # MICac = getMIC(CP, rho, Vinf, xac[1])
    # Mac = MICac'*gamma
    #
    # Cmac = Mac/q/Sref/cref
    # # ----------------------------------------------------------

    # # --------------- structures --------------------------
    # # bending moment distribution
    # Mb = qmvrN*BMM*gamma
    #
    # # distance along structural span
    # ds_str = CP.ds./cos.(CP.sweep)
    # eta_str = [0; cumsum(ds_str, 2)]
    # eta_str = 0.5*(eta_str[1:end-1] + eta_str[2:end])
    # # --------------------------------------------------------------

    # if (plots)
    #   # ------------- plots --------------------
    #   # plot wings
    # #       plot_wing(LE,QC,TE,CP)
    #   N = length(QC.x)
    #   #       axis equal
    #   PyPlot.figure()
    # #       if TE[3] < 0:
    # #         TE.y = -TE.y
    # #       end
    #   for i = 1:N-1
    #     PyPlot.plot([LE.y[i], LE.y[i+1]], -[LE.x[i], LE.x[i+1]], "b")
    #     PyPlot.plot([TE.y[i], TE.y[i+1]],-[TE.x[i], TE.x[i+1]], "b")
    #   end
    #   PyPlot.plot([LE.y[1], TE.y[1]],-[LE.x[1], TE.x[1]], "b")
    # #       PyPlot.plot([LE.y[end], TE.y[end]],-[LE.x[end], TE.x[end]], "b")
    #   for i = 1:N-1
    #       PyPlot.plot([-LE.y[i], -LE.y[i+1]],-[LE.x[i], LE.x[i+1]],"b")
    #       PyPlot.plot([-TE.y[i], -TE.y[i+1]],-[TE.x[i], TE.x[i+1]],"b")
    #   end
    #   PyPlot.plot([-LE.y[1], -TE.y[1]],-[LE.x[1], TE.x[1]], "b")
    # #       PyPlot.plot([-LE.y[end], -TE.y[end]],-[LE.x[end], TE.x[end]],"b")
    #   PyPlot.xlabel("x")
    #   PyPlot.ylabel("y")
    #   PyPlot.title("Plot of wing")
    #
    #   # plot lift
    #   #       figure(50) hold on
    #   PyPlot.figure()
    #   l = 2*gamma/U/(ref.c)
    #   l_mvr = 2*gamma_mvr/U/(ref.c)
    #   eta = linspace(0,0.5,length(l))
    #   PyPlot.plot(eta,l,"b")
    #   PyPlot.plot(eta,l_mvr,"r")
    #   PyPlot.xlabel("xi / b")
    #   PyPlot.ylabel("c_l c / c_ref")
    # #   PyPlot.title("")
    #
    #   # plot cl
    #   PyPlot.figure()
    #   PyPlot.plot(eta,cl,"b")
    #   PyPlot.plot(eta,clmax_dist,"r")
    #   PyPlot.plot(eta,clmax,"r--")
    #   PyPlot.xlabel("xi / b")
    #   PyPlot.ylabel("c_l")
    #   PyPlot.title("Plot of c_l")
    #
    #   # plot bending over thickness
    #   PyPlot.figure()
    #   PyPlot.plot(eta_str/eta_str[end]*0.5,Mb./(CP.tc.*CP.chord))
    #   PyPlot.xlabel("xi_str/b_str")
    #   PyPlot.ylabel("M_b/t")
    #   PyPlot.title("Plot of bending over thickness")
    #   PyPlot.show()
    #   # -------------------------------------------------
    # end

    return L, Di, cl, cllocal, Vinfeff, alphaeff, Vextn
end

end


text:
# Automatically generated using Clang.jl

const CIMGUIPLOT_INCLUDED = nothing

# Skipping MacroDefinition: API __declspec ( dllexport )
# Skipping MacroDefinition: EXTERN extern
# Skipping MacroDefinition: CIMGUI_API EXTERN API
# Skipping MacroDefinition: CONST const

const ImPlotMarker = Cint

struct ImPlotNextItemData
    Colors::NTuple{5, ImVec4}
    LineWeight::Cfloat
    Marker::ImPlotMarker
    MarkerSize::Cfloat
    MarkerWeight::Cfloat
    FillAlpha::Cfloat
    ErrorBarSize::Cfloat
    ErrorBarWeight::Cfloat
    DigitalBitHeight::Cfloat
    DigitalBitGap::Cfloat
    RenderLine::Bool
    RenderFill::Bool
    RenderMarkerLine::Bool
    RenderMarkerFill::Bool
    HasHidden::Bool
    Hidden::Bool
    HiddenCond::ImGuiCond
end

struct ImPlotTick
    PlotPos::Cdouble
    PixelPos::Cfloat
    LabelSize::ImVec2
    TextOffset::Cint
    Major::Bool
    ShowLabel::Bool
    Level::Cint
end

struct ImVector_ImPlotTick
    Size::Cint
    Capacity::Cint
    Data::Ptr{ImPlotTick}
end

struct ImPlotTickCollection
    Ticks::ImVector_ImPlotTick
    TextBuffer::ImGuiTextBuffer
    TotalWidth::Cfloat
    TotalHeight::Cfloat
    MaxWidth::Cfloat
    MaxHeight::Cfloat
    Size::Cint
end

struct ImPlotAnnotation
    Pos::ImVec2
    Offset::ImVec2
    ColorBg::ImU32
    ColorFg::ImU32
    TextOffset::Cint
    Clamp::Bool
end

struct ImVector_ImPlotAnnotation
    Size::Cint
    Capacity::Cint
    Data::Ptr{ImPlotAnnotation}
end

struct ImPlotAnnotationCollection
    Annotations::ImVector_ImPlotAnnotation
    TextBuffer::ImGuiTextBuffer
    Size::Cint
end

struct ImPlotPointError
    X::Cdouble
    Y::Cdouble
    Neg::Cdouble
    Pos::Cdouble
end

struct ImPlotColormapMod
    Colormap::Ptr{ImVec4}
    ColormapSize::Cint
end

struct ImPlotTime
    S::Cint
    Us::Cint
end

const ImPlotDateFmt = Cint
const ImPlotTimeFmt = Cint

struct ImPlotDateTimeFmt
    Date::ImPlotDateFmt
    Time::ImPlotTimeFmt
    UseISO8601::Bool
    Use24HourClock::Bool
end

struct ImBufferWriter
    Buffer::Cstring
    Size::Cint
    Pos::Cint
end

struct ImPlotRange
    Min::Cdouble
    Max::Cdouble
end

struct ImPlotNextPlotData
    XRangeCond::ImGuiCond
    YRangeCond::NTuple{3, ImGuiCond}
    X::ImPlotRange
    Y::NTuple{3, ImPlotRange}
    HasXRange::Bool
    HasYRange::NTuple{3, Bool}
    ShowDefaultTicksX::Bool
    ShowDefaultTicksY::NTuple{3, Bool}
    FitX::Bool
    FitY::NTuple{3, Bool}
    LinkedXmin::Ptr{Cdouble}
    LinkedXmax::Ptr{Cdouble}
    LinkedYmin::NTuple{3, Ptr{Cdouble}}
    LinkedYmax::NTuple{3, Ptr{Cdouble}}
end

const ImPlotFlags = Cint
const ImPlotAxisFlags = Cint
const ImPlotOrientation = Cint

struct ImPlotAxis
    Flags::ImPlotAxisFlags
    PreviousFlags::ImPlotAxisFlags
    Range::ImPlotRange
    Pixels::Cfloat
    Orientation::ImPlotOrientation
    Dragging::Bool
    ExtHovered::Bool
    AllHovered::Bool
    Present::Bool
    HasRange::Bool
    LinkedMin::Ptr{Cdouble}
    LinkedMax::Ptr{Cdouble}
    PickerTimeMin::ImPlotTime
    PickerTimeMax::ImPlotTime
    PickerLevel::Cint
    ColorMaj::ImU32
    ColorMin::ImU32
    ColorTxt::ImU32
    RangeCond::ImGuiCond
    HoverRect::ImRect
end

struct ImVector_int
    Size::Cint
    Capacity::Cint
    Data::Ptr{Cint}
end

struct ImPlotLegendData
    Indices::ImVector_int
    Labels::ImGuiTextBuffer
end

struct ImPlotItem
    ID::ImGuiID
    Color::ImVec4
    NameOffset::Cint
    Show::Bool
    LegendHovered::Bool
    SeenThisFrame::Bool
end

struct ImVector_ImPlotItem
    Size::Cint
    Capacity::Cint
    Data::Ptr{ImPlotItem}
end

struct ImPool_ImPlotItem
    Buf::ImVector_ImPlotItem
    Map::ImGuiStorage
    FreeIdx::ImPoolIdx
end

const ImPlotLocation = Cint

struct ImPlotPlot
    ID::ImGuiID
    Flags::ImPlotFlags
    PreviousFlags::ImPlotFlags
    XAxis::ImPlotAxis
    YAxis::NTuple{3, ImPlotAxis}
    LegendData::ImPlotLegendData
    Items::ImPool_ImPlotItem
    SelectStart::ImVec2
    QueryStart::ImVec2
    QueryRect::ImRect
    Selecting::Bool
    Querying::Bool
    Queried::Bool
    DraggingQuery::Bool
    LegendHovered::Bool
    LegendOutside::Bool
    LegendFlipSideNextFrame::Bool
    FrameHovered::Bool
    PlotHovered::Bool
    ColormapIdx::Cint
    CurrentYAxis::Cint
    MousePosLocation::ImPlotLocation
    LegendLocation::ImPlotLocation
    LegendOrientation::ImPlotOrientation
    FrameRect::ImRect
    CanvasRect::ImRect
    PlotRect::ImRect
    AxesRect::ImRect
end

const ImPlotAxisColor = Cvoid
const ImPlotAxisState = Cvoid

struct ImPlotInputMap
    PanButton::ImGuiMouseButton
    PanMod::ImGuiKeyModFlags
    FitButton::ImGuiMouseButton
    ContextMenuButton::ImGuiMouseButton
    BoxSelectButton::ImGuiMouseButton
    BoxSelectMod::ImGuiKeyModFlags
    BoxSelectCancelButton::ImGuiMouseButton
    QueryButton::ImGuiMouseButton
    QueryMod::ImGuiKeyModFlags
    QueryToggleMod::ImGuiKeyModFlags
    HorizontalMod::ImGuiKeyModFlags
    VerticalMod::ImGuiKeyModFlags
end

struct ImPlotStyle
    LineWeight::Cfloat
    Marker::Cint
    MarkerSize::Cfloat
    MarkerWeight::Cfloat
    FillAlpha::Cfloat
    ErrorBarSize::Cfloat
    ErrorBarWeight::Cfloat
    DigitalBitHeight::Cfloat
    DigitalBitGap::Cfloat
    PlotBorderSize::Cfloat
    MinorAlpha::Cfloat
    MajorTickLen::ImVec2
    MinorTickLen::ImVec2
    MajorTickSize::ImVec2
    MinorTickSize::ImVec2
    MajorGridSize::ImVec2
    MinorGridSize::ImVec2
    PlotPadding::ImVec2
    LabelPadding::ImVec2
    LegendPadding::ImVec2
    LegendInnerPadding::ImVec2
    LegendSpacing::ImVec2
    MousePosPadding::ImVec2
    AnnotationPadding::ImVec2
    FitPadding::ImVec2
    PlotDefaultSize::ImVec2
    PlotMinSize::ImVec2
    Colors::NTuple{24, ImVec4}
    AntiAliasedLines::Bool
    UseLocalTime::Bool
    UseISO8601::Bool
    Use24HourClock::Bool
end

struct ImPlotLimits
    X::ImPlotRange
    Y::ImPlotRange
end

struct ImPlotPoint
    x::Cdouble
    y::Cdouble
end

struct ImVector_ImPlotPlot
    Size::Cint
    Capacity::Cint
    Data::Ptr{ImPlotPlot}
end

struct ImPool_ImPlotPlot
    Buf::ImVector_ImPlotPlot
    Map::ImGuiStorage
    FreeIdx::ImPoolIdx
end

const ImPlotScale = Cint

struct ImVector_ImPlotColormapMod
    Size::Cint
    Capacity::Cint
    Data::Ptr{ImPlotColormapMod}
end

struct ImPlotContext
    Plots::ImPool_ImPlotPlot
    CurrentPlot::Ptr{ImPlotPlot}
    CurrentItem::Ptr{ImPlotItem}
    PreviousItem::Ptr{ImPlotItem}
    XTicks::ImPlotTickCollection
    YTicks::NTuple{3, ImPlotTickCollection}
    YAxisReference::NTuple{3, Cfloat}
    Annotations::ImPlotAnnotationCollection
    Scales::NTuple{3, ImPlotScale}
    PixelRange::NTuple{3, ImRect}
    Mx::Cdouble
    My::NTuple{3, Cdouble}
    LogDenX::Cdouble
    LogDenY::NTuple{3, Cdouble}
    ExtentsX::ImPlotRange
    ExtentsY::NTuple{3, ImPlotRange}
    FitThisFrame::Bool
    FitX::Bool
    FitY::NTuple{3, Bool}
    RenderX::Bool
    RenderY::NTuple{3, Bool}
    ChildWindowMade::Bool
    Style::ImPlotStyle
    ColorModifiers::ImVector_ImGuiColorMod
    StyleModifiers::ImVector_ImGuiStyleMod
    Colormap::Ptr{ImVec4}
    ColormapSize::Cint
    ColormapModifiers::ImVector_ImPlotColormapMod
    Tm::Cint
    VisibleItemCount::Cint
    DigitalPlotItemCnt::Cint
    DigitalPlotOffset::Cint
    NextPlotData::ImPlotNextPlotData
    NextItemData::ImPlotNextItemData
    InputMap::ImPlotInputMap
    MousePos::NTuple{3, ImPlotPoint}
end

const ImPlotCol = Cint
const ImPlotStyleVar = Cint
const ImPlotColormap = Cint
const ImPlotYAxis = Cint
const ImPlotTimeUnit = Cint

struct ImVector_ImS16
    Size::Cint
    Capacity::Cint
    Data::Ptr{ImS16}
end

struct ImVector_ImS32
    Size::Cint
    Capacity::Cint
    Data::Ptr{ImS32}
end

struct ImVector_ImS64
    Size::Cint
    Capacity::Cint
    Data::Ptr{ImS64}
end

struct ImVector_ImS8
    Size::Cint
    Capacity::Cint
    Data::Ptr{ImS8}
end

struct ImVector_ImU16
    Size::Cint
    Capacity::Cint
    Data::Ptr{ImU16}
end

struct ImVector_ImU64
    Size::Cint
    Capacity::Cint
    Data::Ptr{ImU64}
end

struct ImVector_ImU8
    Size::Cint
    Capacity::Cint
    Data::Ptr{ImU8}
end

struct ImVector_double
    Size::Cint
    Capacity::Cint
    Data::Ptr{Cdouble}
end

@cenum ImPlotFlags_::UInt32 begin
    ImPlotFlags_None = 0
    ImPlotFlags_NoTitle = 1
    ImPlotFlags_NoLegend = 2
    ImPlotFlags_NoMenus = 4
    ImPlotFlags_NoBoxSelect = 8
    ImPlotFlags_NoMousePos = 16
    ImPlotFlags_NoHighlight = 32
    ImPlotFlags_NoChild = 64
    ImPlotFlags_Equal = 128
    ImPlotFlags_YAxis2 = 256
    ImPlotFlags_YAxis3 = 512
    ImPlotFlags_Query = 1024
    ImPlotFlags_Crosshairs = 2048
    ImPlotFlags_AntiAliased = 4096
    ImPlotFlags_CanvasOnly = 31
end

@cenum ImPlotAxisFlags_::UInt32 begin
    ImPlotAxisFlags_None = 0
    ImPlotAxisFlags_NoGridLines = 1
    ImPlotAxisFlags_NoTickMarks = 2
    ImPlotAxisFlags_NoTickLabels = 4
    ImPlotAxisFlags_LogScale = 8
    ImPlotAxisFlags_Time = 16
    ImPlotAxisFlags_Invert = 32
    ImPlotAxisFlags_LockMin = 64
    ImPlotAxisFlags_LockMax = 128
    ImPlotAxisFlags_Lock = 192
    ImPlotAxisFlags_NoDecorations = 7
end

@cenum ImPlotCol_::UInt32 begin
    ImPlotCol_Line = 0
    ImPlotCol_Fill = 1
    ImPlotCol_MarkerOutline = 2
    ImPlotCol_MarkerFill = 3
    ImPlotCol_ErrorBar = 4
    ImPlotCol_FrameBg = 5
    ImPlotCol_PlotBg = 6
    ImPlotCol_PlotBorder = 7
    ImPlotCol_LegendBg = 8
    ImPlotCol_LegendBorder = 9
    ImPlotCol_LegendText = 10
    ImPlotCol_TitleText = 11
    ImPlotCol_InlayText = 12
    ImPlotCol_XAxis = 13
    ImPlotCol_XAxisGrid = 14
    ImPlotCol_YAxis = 15
    ImPlotCol_YAxisGrid = 16
    ImPlotCol_YAxis2 = 17
    ImPlotCol_YAxisGrid2 = 18
    ImPlotCol_YAxis3 = 19
    ImPlotCol_YAxisGrid3 = 20
    ImPlotCol_Selection = 21
    ImPlotCol_Query = 22
    ImPlotCol_Crosshairs = 23
    ImPlotCol_COUNT = 24
end

@cenum ImPlotStyleVar_::UInt32 begin
    ImPlotStyleVar_LineWeight = 0
    ImPlotStyleVar_Marker = 1
    ImPlotStyleVar_MarkerSize = 2
    ImPlotStyleVar_MarkerWeight = 3
    ImPlotStyleVar_FillAlpha = 4
    ImPlotStyleVar_ErrorBarSize = 5
    ImPlotStyleVar_ErrorBarWeight = 6
    ImPlotStyleVar_DigitalBitHeight = 7
    ImPlotStyleVar_DigitalBitGap = 8
    ImPlotStyleVar_PlotBorderSize = 9
    ImPlotStyleVar_MinorAlpha = 10
    ImPlotStyleVar_MajorTickLen = 11
    ImPlotStyleVar_MinorTickLen = 12
    ImPlotStyleVar_MajorTickSize = 13
    ImPlotStyleVar_MinorTickSize = 14
    ImPlotStyleVar_MajorGridSize = 15
    ImPlotStyleVar_MinorGridSize = 16
    ImPlotStyleVar_PlotPadding = 17
    ImPlotStyleVar_LabelPadding = 18
    ImPlotStyleVar_LegendPadding = 19
    ImPlotStyleVar_LegendInnerPadding = 20
    ImPlotStyleVar_LegendSpacing = 21
    ImPlotStyleVar_MousePosPadding = 22
    ImPlotStyleVar_AnnotationPadding = 23
    ImPlotStyleVar_FitPadding = 24
    ImPlotStyleVar_PlotDefaultSize = 25
    ImPlotStyleVar_PlotMinSize = 26
    ImPlotStyleVar_COUNT = 27
end

@cenum ImPlotMarker_::Int32 begin
    ImPlotMarker_None = -1
    ImPlotMarker_Circle = 0
    ImPlotMarker_Square = 1
    ImPlotMarker_Diamond = 2
    ImPlotMarker_Up = 3
    ImPlotMarker_Down = 4
    ImPlotMarker_Left = 5
    ImPlotMarker_Right = 6
    ImPlotMarker_Cross = 7
    ImPlotMarker_Plus = 8
    ImPlotMarker_Asterisk = 9
    ImPlotMarker_COUNT = 10
end

@cenum ImPlotColormap_::UInt32 begin
    ImPlotColormap_Default = 0
    ImPlotColormap_Deep = 1
    ImPlotColormap_Dark = 2
    ImPlotColormap_Pastel = 3
    ImPlotColormap_Paired = 4
    ImPlotColormap_Viridis = 5
    ImPlotColormap_Plasma = 6
    ImPlotColormap_Hot = 7
    ImPlotColormap_Cool = 8
    ImPlotColormap_Pink = 9
    ImPlotColormap_Jet = 10
    ImPlotColormap_COUNT = 11
end

@cenum ImPlotLocation_::UInt32 begin
    ImPlotLocation_Center = 0
    ImPlotLocation_North = 1
    ImPlotLocation_South = 2
    ImPlotLocation_West = 4
    ImPlotLocation_East = 8
    ImPlotLocation_NorthWest = 5
    ImPlotLocation_NorthEast = 9
    ImPlotLocation_SouthWest = 6
    ImPlotLocation_SouthEast = 10
end

@cenum ImPlotOrientation_::UInt32 begin
    ImPlotOrientation_Horizontal = 0
    ImPlotOrientation_Vertical = 1
end

@cenum ImPlotYAxis_::UInt32 begin
    ImPlotYAxis_1 = 0
    ImPlotYAxis_2 = 1
    ImPlotYAxis_3 = 2
end

@cenum ImPlotScale_::UInt32 begin
    ImPlotScale_LinLin = 0
    ImPlotScale_LogLin = 1
    ImPlotScale_LinLog = 2
    ImPlotScale_LogLog = 3
end

@cenum ImPlotTimeUnit_::UInt32 begin
    ImPlotTimeUnit_Us = 0
    ImPlotTimeUnit_Ms = 1
    ImPlotTimeUnit_S = 2
    ImPlotTimeUnit_Min = 3
    ImPlotTimeUnit_Hr = 4
    ImPlotTimeUnit_Day = 5
    ImPlotTimeUnit_Mo = 6
    ImPlotTimeUnit_Yr = 7
    ImPlotTimeUnit_COUNT = 8
end

@cenum ImPlotDateFmt_::UInt32 begin
    ImPlotDateFmt_None = 0
    ImPlotDateFmt_DayMo = 1
    ImPlotDateFmt_DayMoYr = 2
    ImPlotDateFmt_MoYr = 3
    ImPlotDateFmt_Mo = 4
    ImPlotDateFmt_Yr = 5
end

@cenum ImPlotTimeFmt_::UInt32 begin
    ImPlotTimeFmt_None = 0
    ImPlotTimeFmt_Us = 1
    ImPlotTimeFmt_SUs = 2
    ImPlotTimeFmt_SMs = 3
    ImPlotTimeFmt_S = 4
    ImPlotTimeFmt_HrMinSMs = 5
    ImPlotTimeFmt_HrMinS = 6
    ImPlotTimeFmt_HrMin = 7
    ImPlotTimeFmt_Hr = 8
end

const ImPlotPoint_getter = Ptr{Cvoid}


text:
# * Orbital equation

"""
    AtomicOrbitalEquation(atom, equation, orbital, ϕ, hamiltonian)

Governs the evolution of an atomic `orbital` belonging to an
`atom`. `equation` is the symbolic expression, from which
`hamiltonian` is constructed. `ϕ` is the `QuasiVector` representing
the radial orbital.
"""
mutable struct AtomicOrbitalEquation{T, B<:Basis,
                                     O<:AbstractOrbital,
                                     A<:Atom{T,B,O},
                                     Equation,
                                     M# <:OrbitalHamiltonian{T,B,O}
                                     }
    atom::A
    equation::Equation
    orbital::O
    index::Int
    ϕ::RadialOrbital{T,B}
    hamiltonian::M
end

SCF.hamiltonian(hfeq::AtomicOrbitalEquation) = hfeq.hamiltonian

function AtomicOrbitalEquation(atom::A, equation::Equation, orbital::O,
                               terms::Vector{<:OrbitalHamiltonianTerm{<:O,<:O,T}},
                               symmetry_orbitals) where {T,
                                                         B<:Basis,
                                                         O<:AbstractOrbital,
                                                         A<:Atom{T,B,O},
                                                         Equation}
    OV = typeof(view(atom,orbital))
    projector = Projector(OV[view(atom, other)
                             for other in symmetry_orbitals],
                          symmetry_orbitals)

    hamiltonian = OrbitalHamiltonian(radial_basis(atom), terms,
                                     atom.mix_coeffs, projector, orbital)

    AtomicOrbitalEquation(atom, equation, orbital, orbital_index(atom, orbital),
                          view(atom, orbital), hamiltonian)
end

"""
    SCF.energy(hfeq::AtomicOrbitalEquation[, which=:total])

Compute the orbital energy for the orbital governed by
`hfeq`. Optionally select which contribution is computed (`:total`,
`:onebody`, `:direct`, or `:exchange`).
"""
SCF.energy(hfeq::AtomicOrbitalEquation, which::Symbol=:total) =
    materialize(applied(*, hfeq.ϕ', hfeq.hamiltonian[which], hfeq.ϕ))

function Base.show(io::IO, hfeq::AtomicOrbitalEquation)
    EHa = energy(hfeq)
    write(io, "⟨$(hfeq.orbital)| 𝓗 |$(hfeq.orbital)⟩ = $(EHa) Ha = $(27.211EHa) eV")
end

function Base.show(io::IO, ::MIME"text/plain", hfeq::AtomicOrbitalEquation)
    write(io, "Hartree–Fock equation: E|$(hfeq.orbital)⟩ = ")
    show(io, hfeq.equation)
    write(io, "\n")
    show(io, hfeq)
    write(io, "\n")
end

import SCF: energy
export energy


text:
using Test
import ArchGDAL;
const AG = ArchGDAL;
import ImageCore

@testset "test_types.jl" begin
    @testset "Testing convert macro" begin
        @testset "Basic conversions" begin
            @test Base.convert(AG.GDALDataType, UInt8) == AG.GDT_Byte
        end

        @testset "Error on conversion non implemented" begin
            @test_throws MethodError Base.convert(AG.GDALDataType, Int64)
            @test_throws KeyError Base.convert(DataType, AG.GDT_TypeCount)
            @test_throws KeyError Base.convert(ImageCore.Normed, AG.GDT_Float32)
            @test_throws KeyError Base.convert(DataType, AG.OFSTMaxSubType)
        end

        @testset "Default types for non bijective conversion mapset" begin
            # Default type for AG.OFTInteger is Int32
            @test Base.convert(AG.OGRFieldType, Bool) == AG.OFTInteger
            @test Base.convert(AG.OGRFieldType, Int16) == AG.OFTInteger
            @test Base.convert(AG.OGRFieldType, Int32) == AG.OFTInteger
            @test Base.convert(DataType, AG.OFTInteger) == Int32

            # Default type for AG.OFTReal is Float64
            @test Base.convert(AG.OGRFieldType, Float32) == AG.OFTReal
            @test Base.convert(AG.OGRFieldType, Float64) == AG.OFTReal
            @test Base.convert(DataType, AG.OFTReal) == Float64
        end
    end

    @testset "Testing GDAL Type Methods" begin
        @testset "GDAL Open Flags" begin
            @test AG.OF_READONLY | 0x04 == 0x04
            @test 0x06 | AG.OF_READONLY == 0x06
            @test AG.OF_READONLY | AG.OF_GNM == AG.OF_READONLY | AG.OF_GNM
        end

        @testset "GDAL Data Types" begin
            @test AG.typesize(AG.GDT_UInt16) == 16
            @test AG.typename(AG.GDT_UInt16) == "UInt16"
            @test AG.gettype("UInt16") == AG.GDT_UInt16

            @test AG.typeunion(AG.GDT_UInt16, AG.GDT_Byte) == AG.GDT_UInt16
            @test AG.iscomplex(AG.GDT_Float32) == false

            @test AG.iscomplex(AG.GDT_CFloat32) == true
            @test AG.iscomplex(AG.GDT_CFloat64) == true
        end

        @testset "GDAL Colors and Palettes" begin
            @test AG.getname(AG.GARIO_COMPLETE) == "COMPLETE"
            @test AG.asyncstatustype("COMPLETE") == AG.GARIO_COMPLETE
            @test AG.asyncstatustype("ERROR") == AG.GARIO_ERROR
            @test AG.asyncstatustype("PENDING") == AG.GARIO_PENDING
            @test AG.asyncstatustype("UPDATE") == AG.GARIO_UPDATE

            @test AG.colorinterp("Hue") == AG.GCI_HueBand
            @test AG.colorinterp("Red") == AG.GCI_RedBand
            @test AG.colorinterp("Blue") == AG.GCI_BlueBand

            @test AG.getname(AG.GPI_Gray) == "Gray"
            @test AG.getname(AG.GPI_RGB) == "RGB"
            @test AG.getname(AG.GPI_CMYK) == "CMYK"
            @test AG.getname(AG.GPI_HLS) == "HLS"
        end

        @testset "GDAL Field Types" begin
            @test AG.getname(AG.OFTString) == "String"
            @test AG.getname(AG.OFTIntegerList) == "IntegerList"
            @test AG.getname(AG.OFSTBoolean) == "Boolean"
            @test AG.getname(AG.OFSTFloat32) == "Float32"

            @test AG.arecompatible(AG.OFTReal, AG.OFSTNone) == true
            @test AG.arecompatible(AG.OFTReal, AG.OFSTBoolean) == false
            @test AG.arecompatible(AG.OFTReal, AG.OFSTInt16) == false
            @test AG.arecompatible(AG.OFTReal, AG.OFSTFloat32) == true
        end

        @testset "Base Geometry Types" begin
            @test AG.basetype(AG.wkbUnknown) == AG.wkbUnknown
            @test AG.basetype(AG.wkbPoint) == AG.wkbPoint
            @test AG.basetype(AG.wkbLineString) == AG.wkbLineString
            @test AG.basetype(AG.wkbPolygon) == AG.wkbPolygon
            @test AG.basetype(AG.wkbMultiPoint) == AG.wkbMultiPoint
            @test AG.basetype(AG.wkbMultiLineString) == AG.wkbMultiLineString
            @test AG.basetype(AG.wkbMultiPolygon) == AG.wkbMultiPolygon
            @test AG.basetype(AG.wkbGeometryCollection) ==
                  AG.wkbGeometryCollection
            @test AG.basetype(AG.wkbCircularString) == AG.wkbCircularString
            @test AG.basetype(AG.wkbCompoundCurve) == AG.wkbCompoundCurve
            @test AG.basetype(AG.wkbCurvePolygon) == AG.wkbCurvePolygon
            @test AG.basetype(AG.wkbMultiCurve) == AG.wkbMultiCurve
            @test AG.basetype(AG.wkbMultiSurface) == AG.wkbMultiSurface
            @test AG.basetype(AG.wkbCurve) == AG.wkbCurve
            @test AG.basetype(AG.wkbSurface) == AG.wkbSurface
            @test AG.basetype(AG.wkbPolyhedralSurface) ==
                  AG.wkbPolyhedralSurface
            @test AG.basetype(AG.wkbTIN) == AG.wkbTIN
            @test AG.basetype(AG.wkbTriangle) == AG.wkbTriangle
            @test AG.basetype(AG.wkbNone) == AG.wkbNone
            @test AG.basetype(AG.wkbLinearRing) == AG.wkbLinearRing
            @test AG.basetype(AG.wkbCircularStringZ) == AG.wkbCircularString
            @test AG.basetype(AG.wkbCompoundCurveZ) == AG.wkbCompoundCurve
            @test AG.basetype(AG.wkbCurvePolygonZ) == AG.wkbCurvePolygon
            @test AG.basetype(AG.wkbMultiCurveZ) == AG.wkbMultiCurve
            @test AG.basetype(AG.wkbMultiSurfaceZ) == AG.wkbMultiSurface
            @test AG.basetype(AG.wkbCurveZ) == AG.wkbCurve
            @test AG.basetype(AG.wkbSurfaceZ) == AG.wkbSurface
            @test AG.basetype(AG.wkbPolyhedralSurfaceZ) ==
                  AG.wkbPolyhedralSurface
            @test AG.basetype(AG.wkbTINZ) == AG.wkbTIN
            @test AG.basetype(AG.wkbTriangleZ) == AG.wkbTriangle
            @test AG.basetype(AG.wkbPointM) == AG.wkbPoint
            @test AG.basetype(AG.wkbLineStringM) == AG.wkbLineString
            @test AG.basetype(AG.wkbPolygonM) == AG.wkbPolygon
            @test AG.basetype(AG.wkbMultiPointM) == AG.wkbMultiPoint
            @test AG.basetype(AG.wkbMultiLineStringM) == AG.wkbMultiLineString
            @test AG.basetype(AG.wkbMultiPolygonM) == AG.wkbMultiPolygon
            @test AG.basetype(AG.wkbGeometryCollectionM) ==
                  AG.wkbGeometryCollection
            @test AG.basetype(AG.wkbCircularStringM) == AG.wkbCircularString
            @test AG.basetype(AG.wkbCompoundCurveM) == AG.wkbCompoundCurve
            @test AG.basetype(AG.wkbCurvePolygonM) == AG.wkbCurvePolygon
            @test AG.basetype(AG.wkbMultiCurveM) == AG.wkbMultiCurve
            @test AG.basetype(AG.wkbMultiSurfaceM) == AG.wkbMultiSurface
            @test AG.basetype(AG.wkbCurveM) == AG.wkbCurve
            @test AG.basetype(AG.wkbSurfaceM) == AG.wkbSurface
            @test AG.basetype(AG.wkbPolyhedralSurfaceM) ==
                  AG.wkbPolyhedralSurface
            @test AG.basetype(AG.wkbTINM) == AG.wkbTIN
            @test AG.basetype(AG.wkbTriangleM) == AG.wkbTriangle
            @test AG.basetype(AG.wkbPointZM) == AG.wkbPoint
            @test AG.basetype(AG.wkbLineStringZM) == AG.wkbLineString
            @test AG.basetype(AG.wkbPolygonZM) == AG.wkbPolygon
            @test AG.basetype(AG.wkbMultiPointZM) == AG.wkbMultiPoint
            @test AG.basetype(AG.wkbMultiLineStringZM) == AG.wkbMultiLineString
            @test AG.basetype(AG.wkbMultiPolygonZM) == AG.wkbMultiPolygon
            @test AG.basetype(AG.wkbGeometryCollectionZM) ==
                  AG.wkbGeometryCollection
            @test AG.basetype(AG.wkbCircularStringZM) == AG.wkbCircularString
            @test AG.basetype(AG.wkbCompoundCurveZM) == AG.wkbCompoundCurve
            @test AG.basetype(AG.wkbCurvePolygonZM) == AG.wkbCurvePolygon
            @test AG.basetype(AG.wkbMultiCurveZM) == AG.wkbMultiCurve
            @test AG.basetype(AG.wkbMultiSurfaceZM) == AG.wkbMultiSurface
            @test AG.basetype(AG.wkbCurveZM) == AG.wkbCurve
            @test AG.basetype(AG.wkbSurfaceZM) == AG.wkbSurface
            @test AG.basetype(AG.wkbPolyhedralSurfaceZM) ==
                  AG.wkbPolyhedralSurface
            @test AG.basetype(AG.wkbTINZM) == AG.wkbTIN
            @test AG.basetype(AG.wkbTriangleZM) == AG.wkbTriangle
            @test AG.basetype(AG.wkbPoint25D) == AG.wkbPoint
            @test AG.basetype(AG.wkbLineString25D) == AG.wkbLineString
            @test AG.basetype(AG.wkbPolygon25D) == AG.wkbPolygon
            @test AG.basetype(AG.wkbMultiPoint25D) == AG.wkbMultiPoint
            @test AG.basetype(AG.wkbMultiLineString25D) == AG.wkbMultiLineString
            @test AG.basetype(AG.wkbMultiPolygon25D) == AG.wkbMultiPolygon
            @test AG.basetype(AG.wkbGeometryCollection25D) ==
                  AG.wkbGeometryCollection
        end
    end
end