text:
################################################################################
#                       INFINITE PARAMETER FUNCTION METHODS
################################################################################
# Extend dispatch_variable_ref
function dispatch_variable_ref(model::InfiniteModel,
                               index::ParameterFunctionIndex
                               )::ParameterFunctionRef
    return ParameterFunctionRef(model, index)
end

# Extend _add_data_object
function _add_data_object(model::InfiniteModel,
                          object::ParameterFunctionData
                          )::ParameterFunctionIndex
    return MOIUC.add_item(model.param_functions, object)
end

# Extend _data_dictionary (reference based)
function _data_dictionary(fref::ParameterFunctionRef
    )::MOIUC.CleverDict{ParameterFunctionIndex, ParameterFunctionData{InfiniteParameterFunction{GeneralVariableRef}}}
    return JuMP.owner_model(fref).param_functions
end

# Extend _data_object
function _data_object(fref::ParameterFunctionRef
    )::ParameterFunctionData{InfiniteParameterFunction{GeneralVariableRef}}
  object = get(_data_dictionary(fref), JuMP.index(fref), nothing)
  object === nothing && error("Invalid infinite parameter function reference, cannot find " *
                        "corresponding object in the model. This is likely " *
                        "caused by using the reference of a deleted function.")
  return object
end

# Extend _core_variable_object
function _core_variable_object(fref::ParameterFunctionRef
    )::InfiniteParameterFunction{GeneralVariableRef}
    return _data_object(fref).func
end

# Extend _object_numbers
function _object_numbers(fref::ParameterFunctionRef)::Vector{Int}
    return _core_variable_object(fref).object_nums
end

# Extend _parameter_numbers
function _parameter_numbers(fref::ParameterFunctionRef)::Vector{Int}
    return _core_variable_object(fref).parameter_nums
end

"""
    build_parameter_function(_error::Function, func::Function,
        parameter_refs::Union{GeneralVariableRef, AbstractArray{<:GeneralVariableRef}, Tuple}
        )::InfiniteParameterFunction{GeneralVariableRef}

Build an [`InfiniteParameterFunction`](@ref) object that employs a infinite
parameter function `func` that takes instances of the infinite parameter(s) as
input. This can ultimately by incorporated into expressions to enable nonlinear
infinite parameter behavior and/or incorporate data over infinite domains. Errors
if the infinite parameter tuple is formatted incorrectly. The allowed format
follows that of infinite variables. Also errors if the function doesn't accept
a support realization of the `parameter_refs` as input.

**Example**
```julia-repl
julia> f = build_parameter_function(error, sin, t)
InfiniteParameterFunction{GeneralVariableRef}(sin, (t,), [1], [1])
```
"""
function build_parameter_function(
    _error::Function,
    func::Function,
    parameter_refs::Union{GeneralVariableRef, AbstractArray{<:GeneralVariableRef}, Tuple}
    )::InfiniteParameterFunction{GeneralVariableRef}
    # process the parameter reference inputs
    prefs = VectorTuple(parameter_refs)
    # check the arguments
    _check_parameter_tuple(_error, prefs)
    _check_valid_function(_error, func, prefs)
    # get the parameter object numbers
    object_nums = Int[]
    for pref in prefs
        union!(object_nums, _object_number(pref))
    end
    # get the parameter numbers
    param_nums = [_parameter_number(pref) for pref in prefs]
    # make the variable and return
    return InfiniteParameterFunction(func, prefs, object_nums, param_nums)
end

# Used to update parameter-parameter function mappings
function _update_param_var_mapping(fref::ParameterFunctionRef,
                                   prefs::VectorTuple)::Nothing
    for pref in prefs
        dependency_list = _parameter_function_dependencies(pref)
        if !(JuMP.index(fref) in dependency_list)
            push!(dependency_list, JuMP.index(fref))
        end
    end
    return
end

"""
    add_parameter_function(model::InfiniteModel, pfunc::InfiniteParameterFunction,
                           [name::String = ""])::GeneralVariableRef

Add an [`InfiniteParameterFunction`](@ref) `pfunc` to the `model` using `name` for
printing and return a `GeneralVariableRef` such that it can be embedded in
expressions. Errors if the infinite parameters `pfunc` points to do not belong to
`model`.

**Example**
```julia-repl
julia> f = build_parameter_function(error, sin, t);

julia> fref = add_parameter_function(model, f, "sin")
sin(t)
```
"""
function add_parameter_function(
    model::InfiniteModel,
    pfunc::InfiniteParameterFunction,
    name::String = ""
    )::GeneralVariableRef
    _check_parameters_valid(model, pfunc.parameter_refs)
    if isempty(name)
        name = String(first(methods(pfunc.func)).name)
    end
    data_object = ParameterFunctionData(pfunc, name)
    findex = _add_data_object(model, data_object)
    fref = ParameterFunctionRef(model, findex)
    _update_param_var_mapping(fref, pfunc.parameter_refs)
    return _make_variable_ref(model, findex)
end

"""
    JuMP.name(fref::ParameterFunctionRef)::String

Extend [`JuMP.name`](@ref JuMP.name(::JuMP.VariableRef)) to return the base name of
`fref`.

**Example**
```julia-repl
julia> name(fref)
"func_name"
```
"""
function JuMP.name(fref::ParameterFunctionRef)::String
    object = get(_data_dictionary(fref), JuMP.index(fref), nothing)
    return object === nothing ? "" : object.name
end

"""
    JuMP.set_name(fref::ParameterFunctionRef, name::String)::Nothing

Extend [`JuMP.set_name`](@ref JuMP.set_name(::JuMP.VariableRef, ::String)) to set
names of infinite parameter functions.

**Example**
```julia-repl
julia> set_name(fref, "func_name")

julia> name(fref)
"func_name"
```
"""
function JuMP.set_name(fref::ParameterFunctionRef, name::String)::Nothing
    _data_object(fref).name = name
    return
end

"""
    parameter_function(func::Function,
                       pref_inputs::Union{GeneralVariableRef, AbstractArray{GeneralVariableRef}, Tuple};
                       [name::String = ""])::GeneralVariableRef

Make an infinite parameter function and return a `GeneralVariableRef` that can be
embedded in InfiniteOpt expressions. This serves as a convenient wrapper for
[`build_parameter_function`](@ref) and [`add_parameter_function`](@ref). Here
`func` denotes the function that will take a support of infinite parameters as
input (formatted like `pref_inputs`) and will return a scalar value. Errors if
`func` will not take such as support as input or if `pref_inputs` follow an invalid
input format.

**Example**
```julia-repl
julia> p_func = parameter_function(sin, t)
sin(t)
```
"""
function parameter_function(func::Function, pref_inputs;
                            name::String = "")::GeneralVariableRef
    f = build_parameter_function(error, func, pref_inputs)
    model = JuMP.owner_model(first(f.parameter_refs))
    return add_parameter_function(model, f, name)
end

"""
    raw_parameter_refs(fref::ParameterFunctionRef)::VectorTuple{GeneralVariableRef}

Return the raw [`VectorTuple`](@ref) of the parameter references that `fref`
depends on. This is primarily an internal method where
[`parameter_refs`](@ref parameter_refs(::ParameterFunctionRef))
is intended as the preferred user function.
"""
function raw_parameter_refs(fref::ParameterFunctionRef)::VectorTuple{GeneralVariableRef}
    return _core_variable_object(fref).parameter_refs
end

"""
    parameter_refs(fref::ParameterFunctionRef)::Tuple

Return the parameter references associated with `fref`. This
is formatted as a Tuple of containing the parameter references as they inputted
to define `fref`.

**Example**
```julia-repl
julia> parameter_refs(p_func)
(t,)
```
"""
function parameter_refs(fref::ParameterFunctionRef)::Tuple
    return Tuple(raw_parameter_refs(fref))
end

"""
    parameter_list(fref::ParameterFunctionRef)::Vector{GeneralVariableRef}

Return a vector of the parameter references that `fref` depends on. This is
primarily an internal method where [`parameter_refs`](@ref parameter_refs(::ParameterFunctionRef))
is intended as the preferred user function.
"""
function parameter_list(fref::ParameterFunctionRef)::Vector{GeneralVariableRef}
    return raw_parameter_refs(fref).values
end

"""
    raw_function(fref::ParameterFunctionRef)::Function

Returns the raw function behind `fref` that takes a particular support of `fref`'s
infinite parameters as input.
"""
function raw_function(fref::ParameterFunctionRef)::Function
    return _core_variable_object(fref).func
end

# Extend _reduced_variable_dependencies
function _reduced_variable_dependencies(
    fref::ParameterFunctionRef
     )::Vector{ReducedVariableIndex}
    return _data_object(fref).reduced_var_indices
end

# Extend _derivative_dependencies
function _derivative_dependencies(
    fref::ParameterFunctionRef
     )::Vector{DerivativeIndex}
    return _data_object(fref).derivative_indices
end

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

# Extend _constraint_dependencies
function _constraint_dependencies(fref::ParameterFunctionRef)::Vector{ConstraintIndex}
    return _data_object(fref).constraint_indices
end

"""
    used_by_reduced_variable(fref::ParameterFunctionRef)::Bool

Return a `Bool` indicating if `fref` is used by a reduced infinite variable.

**Example**
```julia-repl
julia> used_by_reduced_variable(fref)
false
```
"""
function used_by_reduced_variable(fref::ParameterFunctionRef)::Bool
    return !isempty(_reduced_variable_dependencies(fref))
end

"""
    used_by_derivative(fref::ParameterFunctionRef)::Bool

Return a `Bool` indicating if `fref` is used by a derivative.

**Example**
```julia-repl
julia> used_by_derivative(vref)
true
```
"""
function used_by_derivative(fref::ParameterFunctionRef)::Bool
    return !isempty(_derivative_dependencies(fref))
end

"""
    used_by_measure(fref::ParameterFunctionRef)::Bool

Return a `Bool` indicating if `fref` is used by a measure.

**Example**
```julia-repl
julia> used_by_measure(fref)
true
```
"""
function used_by_measure(fref::ParameterFunctionRef)::Bool
    return !isempty(_measure_dependencies(fref))
end

"""
    used_by_constraint(fref::ParameterFunctionRef)::Bool

Return a `Bool` indicating if `fref` is used by a constraint.

**Example**
```julia-repl
julia> used_by_constraint(fref)
false
```
"""
function used_by_constraint(fref::ParameterFunctionRef)::Bool
    return !isempty(_constraint_dependencies(fref))
end

"""
    is_used(fref::ParameterFunctionRef)::Bool

Return a `Bool` indicating if `fref` is used in the model.

**Example**
```julia-repl
julia> is_used(fref)
true
```
"""
function is_used(fref::ParameterFunctionRef)::Bool
    return used_by_measure(fref) || used_by_constraint(fref) ||
           used_by_reduced_variable(fref) || used_by_derivative(fref)
end

"""
    JuMP.delete(model::InfiniteModel, fref::ParameterFunctionRef)::Nothing

Extend [`JuMP.delete`](@ref JuMP.delete(::JuMP.Model, ::JuMP.VariableRef)) to delete
infinite parameter functions and their dependencies. Errors if `fref` is invalid,
meaning it has already been deleted or it belongs to another model.
"""
function JuMP.delete(model::InfiniteModel, fref::ParameterFunctionRef)::Nothing
    @assert JuMP.is_valid(model, fref) "Infinite parameter function is invalid."
    # update the optimizer model status
    if is_used(fref)
        set_optimizer_model_ready(model, false)
    end
    # update parameter mapping
    all_prefs = parameter_list(fref)
    for pref in all_prefs
        filter!(e -> e != JuMP.index(fref), _parameter_function_dependencies(pref))
    end
    gvref = _make_variable_ref(model, JuMP.index(fref))
    # delete associated reduced variables and mapping
    for index in _reduced_variable_dependencies(fref)
        JuMP.delete(model, dispatch_variable_ref(model, index))
    end
    # delete associated derivative variables and mapping
    for index in _derivative_dependencies(fref)
        JuMP.delete(model, dispatch_variable_ref(model, index))
    end
    # remove from measures if used
    for mindex in _measure_dependencies(fref)
        mref = dispatch_variable_ref(model, mindex)
        func = measure_function(mref)
        data = measure_data(mref)
        if func isa GeneralVariableRef
            new_func = zero(JuMP.GenericAffExpr{Float64, GeneralVariableRef})
            new_meas = Measure(new_func, data, Int[], Int[], true)
        else
            _remove_variable(func, gvref)
            new_meas = build_measure(func, data)
        end
        _set_core_variable_object(mref, new_meas)
    end
    # remove from constraints if used
    for cindex in _constraint_dependencies(fref)
        cref = _temp_constraint_ref(model, cindex)
        func = JuMP.jump_function(JuMP.constraint_object(cref))
        if func isa GeneralVariableRef
            set = JuMP.moi_set(JuMP.constraint_object(cref))
            new_func = zero(JuMP.GenericAffExpr{Float64, GeneralVariableRef})
            new_constr = JuMP.ScalarConstraint(new_func, set)
            _set_core_constraint_object(cref, new_constr)
            empty!(_object_numbers(cref))
        else
            _remove_variable(func, gvref)
            _data_object(cref).object_nums = sort(_object_numbers(func))
        end
    end
    # delete the data object
    _delete_data_object(fref)
    return
end

################################################################################
#                              COMPARISON METHODS
################################################################################
## Extend for better comparisons than default
# GenericAffExpr
function Base.:(==)(aff1::JuMP.GenericAffExpr{C, V},
    aff2::JuMP.GenericAffExpr{C, V}
    )::Bool where {C, V <: GeneralVariableRef}
    return aff1.constant == aff2.constant && aff1.terms == aff2.terms
end

# GenericQuadExpr
function Base.:(==)(quad1::JuMP.GenericQuadExpr{C, V},
    quad2::JuMP.GenericQuadExpr{C, V}
    )::Bool where {C, V <: GeneralVariableRef}
    pairs1 = collect(pairs(quad1.terms))
    pairs2 = collect(pairs(quad2.terms))
    if length(pairs1) != length(pairs2)
        return false
    end
    for i in eachindex(pairs1)
        if pairs1[i][1].a != pairs2[i][1].a || pairs1[i][1].b != pairs2[i][1].b ||
           pairs1[i][2] != pairs2[i][2]
            return false
        end
    end
    return quad1.aff == quad2.aff
end

################################################################################
#                            VARIABLE LIST MAKING
################################################################################
## Determine which variables are present in a function
# GeneralVariableRef
function _all_function_variables(
    f::GeneralVariableRef
    )::Vector{GeneralVariableRef}
    return [f]
end

# GenericAffExpr
function _all_function_variables(
    f::JuMP.GenericAffExpr{C, V}
    )::Vector{V} where {C, V <: GeneralVariableRef}
    return collect(keys(f.terms))
end

# GenericQuadExpr
function _all_function_variables(
    f::JuMP.GenericQuadExpr{C, V}
    )::Vector{V} where {C, V <: GeneralVariableRef}
    vref_set = Set(keys(f.aff.terms))
    for pair in keys(f.terms)
        push!(vref_set, pair.a)
        push!(vref_set, pair.b)
    end
    return collect(vref_set)
end

# Fallback
function _all_function_variables(f)
    error("`_all_function_variables` not defined for expression of type $(typeof(f)).")
end

################################################################################
#                            OBJECT NUMBER METHODS
################################################################################
## Return the unique set of object numbers in an expression
# Dispatch fallback (--> should be defined for each non-empty variable type)
_object_numbers(expr::DispatchVariableRef)::Vector{Int} = Int[]

# GeneralVariableRef
function _object_numbers(expr::GeneralVariableRef)::Vector{Int}
    return _object_numbers(dispatch_variable_ref(expr))
end

# GenericAffExpr
function _object_numbers(expr::JuMP.GenericAffExpr)::Vector{Int}
    obj_nums = Set{Int}()
    for vref in keys(expr.terms)
        union!(obj_nums, _object_numbers(vref))
    end
    return collect(obj_nums)
end

# GenericQuadExpr
function _object_numbers(expr::JuMP.GenericQuadExpr)::Vector{Int}
    obj_nums = Set{Int}()
    for vref in keys(expr.aff.terms)
        union!(obj_nums, _object_numbers(vref))
    end
    for pair in keys(expr.terms)
        union!(obj_nums, _object_numbers(pair.a))
        union!(obj_nums, _object_numbers(pair.b))
    end
    return collect(obj_nums)
end

# Variable list
function _object_numbers(vrefs::Vector{GeneralVariableRef})::Vector{Int}
    obj_nums = Set{Int}()
    for vref in vrefs
        union!(obj_nums, _object_numbers(vref))
    end
    return collect(obj_nums)
end

################################################################################
#                             PARAMETER NUMBER METHODS
################################################################################
## Return the unique set of parameter numbers in an expression
# Dispatch fallback (--> should be defined for each non-empty variable type)
_parameter_numbers(expr::DispatchVariableRef)::Vector{Int} = Int[]

# GeneralVariableRef
function _parameter_numbers(expr::GeneralVariableRef)::Vector{Int}
    return _parameter_numbers(dispatch_variable_ref(expr))
end

# GenericAffExpr
function _parameter_numbers(expr::JuMP.GenericAffExpr)::Vector{Int}
    par_nums = Set{Int}()
    for vref in keys(expr.terms)
        union!(par_nums, _parameter_numbers(vref))
    end
    return collect(par_nums)
end

# GenericQuadExpr
function _parameter_numbers(expr::JuMP.GenericQuadExpr)::Vector{Int}
    par_nums = Set{Int}()
    for vref in keys(expr.aff.terms)
        union!(par_nums, _parameter_numbers(vref))
    end
    for pair in keys(expr.terms)
        union!(par_nums, _parameter_numbers(pair.a))
        union!(par_nums, _parameter_numbers(pair.b))
    end
    return collect(par_nums)
end

################################################################################
#                             MODEL EXTRACTION METHODS
################################################################################
## Get the model from an expression
# GeneralVariableRef
function _model_from_expr(expr::GeneralVariableRef)::InfiniteModel
    return JuMP.owner_model(expr)
end

# AffExpr
function _model_from_expr(expr::JuMP.GenericAffExpr)::Union{InfiniteModel, Nothing}
    if isempty(expr.terms)
        return
    else
        return JuMP.owner_model(first(keys(expr.terms)))
    end
end

# QuadExpr
function _model_from_expr(expr::JuMP.GenericQuadExpr)::Union{InfiniteModel, Nothing}
    result = _model_from_expr(expr.aff)
    if result !== nothing
        return result
    elseif isempty(expr.terms)
        return
    else
        return JuMP.owner_model(first(keys(expr.terms)).a)
    end
end

# Vector{GeneralVariableRef}
function _model_from_expr(vrefs::Vector{GeneralVariableRef})::Union{InfiniteModel, Nothing}
    if isempty(vrefs)
        return
    else
        return JuMP.owner_model(first(vrefs))
    end
end

# Fallback
function _model_from_expr(expr)
    error("`_model_from_expr` not defined for expr of type $(typeof(expr)).")
end

################################################################################
#                            VARIABLE REMOVAL BOUNDS
################################################################################
## Delete variables from an expression
# GenericAffExpr
function _remove_variable(f::JuMP.GenericAffExpr,
                          vref::GeneralVariableRef)::Nothing
    if haskey(f.terms, vref)
        delete!(f.terms, vref)
    end
    return
end

# GenericQuadExpr
function _remove_variable(f::JuMP.GenericQuadExpr,
                          vref::GeneralVariableRef)::Nothing
    _remove_variable(f.aff, vref)
    for (pair, coef) in f.terms # TODO should we preserve non-vref in pair? --> this would differ from JuMP
        if pair.a == vref || pair.b == vref
            delete!(f.terms, pair)
        end
    end
    return
end

################################################################################
#                         COEFFICIENT MODIFICATION METHODS
################################################################################
## Modify linear coefficient of variable in expression
# GeneralVariableRef
function _set_variable_coefficient!(expr::GeneralVariableRef,
    var::GeneralVariableRef,
    coeff::Real
    )::JuMP.GenericAffExpr{Float64, GeneralVariableRef}
    # Determine if variable is that of the expression and change accordingly
    if expr == var
        return Float64(coeff) * var
    else
        return expr + Float64(coeff) * var
    end
end

# GenericAffExpr
function _set_variable_coefficient!(expr::JuMP.GenericAffExpr{C, V},
    var::V,
    coeff::Real
    )::JuMP.GenericAffExpr{C, V} where {C, V <: GeneralVariableRef}
    # Determine if variable is in the expression and change accordingly
    if haskey(expr.terms, var)
        expr.terms[var] = coeff
        return expr
    else
        return expr + coeff * var
    end
end

# GenericQuadExpr
function _set_variable_coefficient!(expr::JuMP.GenericQuadExpr{C, V},
    var::V,
    coeff::Real
    )::JuMP.GenericQuadExpr{C, V} where {C, V <: GeneralVariableRef}
    # Determine if variable is in the expression and change accordingly
    if haskey(expr.aff.terms, var)
        expr.aff.terms[var] = coeff
        return expr
    else
        return expr + coeff * var
    end
end

# Fallback
function _set_variable_coefficient!(expr, var::GeneralVariableRef, coeff::Real)
    error("Unsupported function type for coefficient modification.")
end

################################################################################
#                         PARAMETER REFERNCE METHODS
################################################################################
## Return an element of a parameter reference tuple given the model and index
# IndependentParameterIndex
function _make_param_tuple_element(model::InfiniteModel,
    idx::IndependentParameterIndex,
    )::GeneralVariableRef
    return _make_parameter_ref(model, idx)
end

# DependentParametersIndex
function _make_param_tuple_element(model::InfiniteModel,
    idx::DependentParametersIndex,
    )::Vector{GeneralVariableRef}
    dpref = DependentParameterRef(model, DependentParameterIndex(idx, 1))
    num_params = _num_parameters(dpref)
    return [GeneralVariableRef(model, idx.value, DependentParameterIndex, i)
            for i in 1:num_params]
end

"""
    parameter_refs(expr)::Tuple

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

**Example**
```julia-repl
julia> parameter_refs(my_expr)
(t,)
```
"""
function parameter_refs(expr::Union{JuMP.GenericAffExpr, JuMP.GenericQuadExpr})::Tuple
    model = _model_from_expr(expr)
    if model === nothing
        return ()
    else
        obj_nums = sort!(_object_numbers(expr))
        obj_indices = _param_object_indices(model)[obj_nums]
        return Tuple(_make_param_tuple_element(model, idx) for idx in obj_indices)
    end
end

################################################################################
#                             EXPRESSION SEARCHING
################################################################################
# Check expression for a particular variable type via a recursive search
#=
function _has_variable(vrefs::Vector{GeneralVariableRef},
                       vref::GeneralVariableRef; prior=GeneralVariableRef[]
                       )::Bool
    if vrefs[1] == vref
        return true
    elseif _index_type(vrefs[1]) == MeasureIndex
        dvref = dispatch_variable_ref(vrefs[1])
        if length(vrefs) > 1
            return _has_variable(_all_function_variables(measure_function(dvref)),
                                 vref, prior = append!(prior, vrefs[2:end]))
        else
            return _has_variable(_all_function_variables(measure_function(dvref)),
                                 vref, prior = prior)
        end
    elseif length(vrefs) > 1
        return _has_variable(vrefs[2:end], vref, prior = prior)
    elseif length(prior) > 0
        return _has_variable(prior, vref)
    else
        return false
    end
end
=#


text:
"""
    WDError <: Exception
A Web Driver Exception.
"""
struct WDError <: Exception
    status::Int
    error::String
    message::String
end

function WDError(err::StatusError)
    json = JSON3.read(err.response.body)
    WDError(err.status, json.error, json.message)
end

ext: jl


text:
# Allow `UniformScaling`s as value for σ
_scale(w, σ) = σ*w
_scale(w::Number, σ::UniformScaling) = σ.λ*w

"""
    SDESolver

Abstract (super-)type for solving methods for stochastic differential equations.
"""
abstract type SDESolver
end

struct EulerMaruyama <: SDESolver
end
const Euler = EulerMaruyama

"""
    EulerMaruyama() <: SDESolver

Euler-Maruyama scheme. `Euler` is defined as alias.
"""
Euler, EulerMaruyama

"""
    StochasticHeun() <: SDESolver

Stochastic heun scheme.
"""
struct StochasticHeun <: SDESolver
end

"""
    Mdb() <: SDESolver

Euler scheme with the diffusion coefficient correction of the modified diffusion bridge.
"""
struct Mdb <: SDESolver
end

"""
    StochasticRungeKutta() <: SDESolver

Stochastic Runge-Kutta scheme for `T<:Number`-valued processes.
"""
struct StochasticRungeKutta <: SDESolver
end

struct EulerMaruyama! <: SDESolver
end

endpoint(y, P) = y

function solve!(::StochasticHeun, Y, u, W::SamplePath, P::ProcessOrCoefficients)

    N = length(W)
    N != length(Y) && error("Y and W differ in length.")

    ww = W.yy
    tt = Y.tt
    yy = Y.yy
    tt[:] = W.tt

    y = u

    for i in 1:N-2 # fix me
        yy[.., i] = y
        B = _b((i,tt[i]), y, P)
        y2 = y + B*(tt[i+1]-tt[i])
        y = y + 0.5*(_b((i+1,tt[i+1]), y2, P) + B)*(tt[i+1]-tt[i]) + σ(tt[i], y, P)*(ww[.., i+1]-ww[..,i])
    end
    yy[.., N-1] = endpoint(y, P)
    Y
end

"""
    solve(method::SDESolver, u, W::SamplePath, P) -> X
    solve(method::SDESolver, u, W::SamplePath, (b, σ)) -> X

Solve stochastic differential equation ``dX_t = b(t,X_t)dt + σ(t,X_t)dW_t``
using `method` in place.

# Example

```
solve(EulerMaruyama(), 1.0, sample(0:0.1:10, Wiener()), ((t,x)->-x, (t,x)->I))
```

```
struct OU <: ContinuousTimeProcess{Float64}
    μ::Float64
end
Bridge.b(s, x, P::OU) = -P.μ*x
Bridge.σ(s, x, P::OU) = I

solve(EulerMaruyama(), 1.0, sample(0:0.1:10, Wiener()), OU(1.4))
```
"""
solve(method::SDESolver, u::T, W::SamplePath, P::ProcessOrCoefficients) where {T} =
    solve!(method, samplepath(W.tt, zero(u)), u, W, P)

"""
    solve(method::SDESolver, u, W::VSamplePath, P) -> X

Solve stochastic differential equation ``dX_t = b(t,X_t)dt + σ(t,X_t)dW_t``
using `method`.
"""
solve(method::SDESolver, u, W::VSamplePath{T}, P::ProcessOrCoefficients) where {T} =
    solve!(method, VSamplePath(W.tt, zeros(T, size(u)..., length(W.tt))), u, W, P)

"""
    solve!(::EulerMaruyama, Y, u, W, P) -> X

Solve stochastic differential equation ``dX_t = b(t,X_t)dt + σ(t,X_t)dW_t``
using the Euler-Maruyama scheme in place.
"""
function solve!(::EulerMaruyama, Y, u::T, W::SamplePath, P::ProcessOrCoefficients) where {T}
    N = length(W)
    N != length(Y) && error("Y and W differ in length.")

    ww = W.yy
    tt = Y.tt
    yy = Y.yy
    tt[:] = W.tt

    y::T = u

    for i in 1:N-1
        yy[.., i] = y
        y = y + _b((i,tt[i]), y, P)*(tt[i+1]-tt[i]) + _scale((ww[.., i+1]-ww[..,i]), σ(tt[i], y, P))
    end
    yy[.., N] = endpoint(y, P)
    Y
end

# fallback method
function solve!(::EulerMaruyama, Y, u::T, W::AbstractPath, P::ProcessOrCoefficients) where {T}
    N = length(W)
    N != length(Y) && error("Y and W differ in length.")

    ww = W.yy
    tt = Y.tt
    yy = Y.yy
    tt[:] = W.tt

    y::T = u

    for (i, t, dt, dw) in increments(W)
        yy[.., i] = y
        y = y + _b((i,t), y, P)*dt + _scale(dw, σ(t, y, P))
    end
    yy[.., N] = endpoint(y, P)
    Y
end

"""
    solve(method, W, P) -> Y

Integrate with `method`, where `P` is a bridge proposal from
`startpoint(P)`` to `endpoint(P)`.

# Examples

```
cs = Bridge.CSpline(tt[1], tt[end], Bridge.b(tt[1], v[1], P),  Bridge.b(tt[end], v[2], P))
P° = BridgeProp(Pσ, v), Pσ.a, cs)
W = sample(tt, Wiener())
solve(Euler(), W, P°)
```
"""
solve(method, W::SamplePath, P) = let u = startpoint(P), Y = samplepath(W.tt, zero(u)); solve!(method, Y, u, W, P); Y end
solve!(method, Y::SamplePath, W::SamplePath, P::ContinuousTimeProcess) = solve!(method, Y, startpoint(P), W, P)

#### Guided Bridges
endpoint(y, P::GuidedBridge) =
 norm(P.H♢[end], 1) < eps() ? P.V[end] : y

solve(::Euler, u, W::SamplePath, P::Union{GuidedBridge,PartialBridge,PartialBridgeνH}) = let X = samplepath(W.tt, zero(u)); solve!(Euler(), X, u, W, P); X end
function solve!(::Euler, Y, u, W::SamplePath, P::Union{GuidedBridge,PartialBridge,PartialBridgeνH})
    W.tt === P.tt && error("Time axis mismatch between bridge P and driving W.") # not strictly an error

    N = length(W)
    N != length(Y) && error("Y and W differ in length.")

    ww = W.yy
    tt = Y.tt
    yy = Y.yy
    tt[:] = P.tt

    y::typeof(u) = u
    if typeof(u) != valtype(P)
    #    @warn "Starting point not of valtype." maxlog=10
    end
    for i in 1:N-1
        yy[.., i] = y
        y = y + _b((i, tt[i]), y, P)*(tt[i+1]-tt[i]) + _scale((ww[.., i+1]-ww[..,i]), σ(tt[i], y, P))
    end
    yy[.., N] = endpoint(y, P)
    yy[.., N]
end

####

"""
    solve!(method, Y, W, P) -> Y

Integrate with `method`, where `P is a bridge proposal overwriting `Y`.
"""
solve!

function solve!(::Mdb, Y, u, W::SamplePath, P::ContinuousTimeProcess{T}) where {T}
    W.tt === P.tt && error("Time axis mismatch between bridge P and driving W.") # not strictly an error

    N = length(W)
    N != length(Y) && error("Y and W differ in length.")

    ww = W.yy
    tt = Y.tt
    yy = Y.yy
    tt[:] = P.tt

    y::T = u

    for i in 1:N-1
        yy[.., i] = y
        y = y + _b((i, tt[i]), y, P)*(tt[i+1]-tt[i]) + _scale((ww[.., i+1]-ww[..,i]), σ(tt[i], y, P)*sqrt((tt[end]-tt[i+1])/(tt[end]-tt[i])))
    end
    yy[.., N] = endpoint(y, P)
    Y
end

function solve!(::StochasticRungeKutta, Y, u::T, W::SamplePath, P::ProcessOrCoefficients) where T<:Number

    N = length(W)
    N != length(Y) && error("Y and W differ in length.")

    ww = W.yy
    tt = Y.tt
    yy = Y.yy
    tt[:] = W.tt

    y = u

    for i in 1:N-1
        yy[.., i] = y
        delta = tt[i+1]-tt[i]
        sqdelta = sqrt(delta)
        B = _b((i,tt[i]), y, P)
        S = σ(tt[i], y, P)
        dw = ww[.., i+1]-ww[..,i]
        y = y + B*delta + S*dw
        ups = y + B*delta + S*sqdelta
        y = y + 0.5(σ(tt[i+1], ups, P) - S)*(dw^2 - delta)/sqdelta

    end
    yy[.., N] = y
    SamplePath{T}(tt, yy)
end

innovations(method, Y, P) = innovations!(method, copy(Y), Y, P)
function innovations!(::EulerMaruyama, W, Y::SamplePath, P)

    N = length(W)
    N != length(Y) && error("Y and W differ in length.")

    yy = Y.yy
    tt = Y.tt
    ww = W.yy
    W.tt[:] = Y.tt

    w = zero(ww[.., 1])

    for i in 1:N-1
        ww[.., i] = w
        w = w + inv(σ(tt[i], yy[.., i], P))*(yy[.., i+1] - yy[.., i] - _b((i,tt[i]), yy[.., i], P)*(tt[i+1]-tt[i]))
    end
    ww[.., N] = w
    W
end

function innovations!(::Mdb, W, Y::SamplePath, P)

    N = length(W)
    N != length(Y) && error("Y and W differ in length.")

    yy = Y.yy
    tt = Y.tt
    ww = W.yy
    W.tt[:] = Y.tt

    w = zero(ww[.., 1])

    for i in 1:N-1
        ww[.., i] = w
        w = w + sqrt((tt[end]-tt[i+1])/(tt[end]-tt[i]))\inv(σ(tt[i], yy[.., i], P))*(yy[.., i+1] - yy[.., i] - _b((i,tt[i]), yy[.., i], P)*(tt[i+1]-tt[i]))
    end
    ww[.., N] = w
    W
end


text:
#=
Examples of loops.
=#

foods = ["melon", "apple", "berry", "orange"]

# Loop through the items in an array
for f in foods
    println("Food item is $f")
end

# Loop through a sequence of numbers
for x in 1:10
    println("x is $x")
end

max_stars_repo_name: wigging/julia-computing

max_issues_repo_name: wigging/julia-computing

max_forks_repo_name: wigging/julia-computing


text:
@testset "prob_of_basepairs" begin
    sequence  = "GGGAACCUC"
    structure = "(((...)))"
    pt = Pairtable(structure)
    fc = FoldCompound(sequence)
    partfn(fc)
    bppm = bpp(fc)
    @test prob_of_basepairs(bppm, pt) ==
        prob_of_basepairs(bppm, structure) ==
        prob_of_basepairs(fc, pt) ==
        prob_of_basepairs(fc, structure) ==
        prob_of_basepairs(sequence, pt) ==
        prob_of_basepairs(sequence, structure)
end


text:
"""
    (A::KeyedArray)("a", 2.0, :γ) == A[1, 2, 3]
    A(:γ) == view(A, :, :, 3)

`KeyedArray`s are callable, and this behaves much like indexing,
except that it searches for the given keys in `axiskeys(A)`,
instead of `axes(A)` for indices.

A single key may be used to indicate a slice, provided that its type
only matches the eltype of one `axiskeys(A,d)`.
You can also slice explicitly with `A("a", :, :)`, both of these return a `view`.

An extra trailing colon (when all other indices are fixed) will return
a zero-dimensional `view`. This allows setting one value by
writing `A("a", 2.0, :γ, :) .= 100`.

Also accepts functions like `A(<=(2.0))` and selectors,
see `Nearest` and `Index`.
"""
@inline @propagate_inbounds (A::KeyedArray)(args...) = getkey(A, args...)

@inline function getkey(A, args...)
    if length(args) == ndims(A)
        inds_raw = map(findindex, args, axiskeys(A))
        inds = Base.to_indices(A, inds_raw)
        @boundscheck checkbounds(A, inds...)
        return @inbounds get_or_view(A, inds...)

    elseif length(args) > ndims(A) && all(args[ndims(A)+1:end] .== (:)) # trailing colons
        args_nd = args[1:ndims(A)]
        inds_raw = map(findindex, args_nd, axiskeys(A))
        inds = Base.to_indices(A, inds_raw)
        @boundscheck checkbounds(A, inds...)
        if inds isa NTuple{<:Any, Int}
            return @inbounds view(keyless(A), inds...) # zero-dim view of underlying
        else
            return @inbounds get_or_view(A, inds...)
        end

    elseif length(args)==1
        arg = first(args)
        d = inferdim(arg, axiskeys(A))
        i = findindex(arg, axiskeys(A,d))
        inds = ntuple(n -> n==d ? i : (:), ndims(A))
        @boundscheck checkbounds(A, inds...)
        return @inbounds view(A, inds...)

    end

    if length(args) != ndims(A)
        throw(ArgumentError(string("wrong number of keys: got ", length(args),
            " arguments, expected ndims(A) = ", ndims(A)," and perhaps a trailing colon.")))
    else
        throw(ArgumentError("can't understand what to do with $args, sorry"))
    end
end

@propagate_inbounds get_or_view(A, inds::Integer...) = getindex(A, inds...)
@propagate_inbounds get_or_view(A, inds...) = view(A, inds...)

@propagate_inbounds function setkey!(A, val, args...)
    length(args) == ndims(A) || error("wrong number of keys")
    inds = map((v,r) -> findindex(v,r), args, axiskeys(A))
    setindex!(A, val, inds...)
end

"""
    findindex(key, vec)

This is usually `findfirst(isequal(key), vec)`, and will error if it finds `nothing`.
But it also understands `findindex(:, vec) = (:)`,
and `findindex(array, vec) = vcat((findindex(x, vec) for x in array)...)`.

It also understands functions `findindex(<(4), vec) = findall(x -> x<4, vec)`,
and selectors like `Nearest(key)` and `Interval(lo,hi)`.
"""
@inline function findindex(a, r::AbstractArray)
    i = findfirst(isequal(a), r)
    i === nothing && throw(ArgumentError("could not find key $(repr(a)) in vector $r"))
    i
end

findindex(a::Colon, r::AbstractArray) = Colon()

findindex(a::Union{AbstractArray, Base.Generator}, r::AbstractArray) =
    [findindex(x, r) for x in a]

findindex(f::Function, r::AbstractArray) = findall(f, r)

# Faster than Base.findfirst(==(i), 1:10) etc:
# Now also https://github.com/JuliaLang/julia/pull/35256

findindex(i::Int, r::Base.OneTo{Int}) = 1 <= i <= r.stop ? i : nothing

findindex(i::Int, r::AbstractUnitRange) = first(r) <= i <= last(r) ? 1+Int(i - first(r)) : nothing

# Faster than Base.findall(==(i), 1:10) etc,
# but returning a range not an Array:

function findindex(eq::Base.Fix2{<:Union{typeof(==),typeof(isequal)},Int}, r::Base.OneTo{Int})
    1 <= eq.x <= r.stop ? (eq.x:eq.x) : (1:0)
end

function findindex(eq::Base.Fix2{<:Union{typeof(==),typeof(isequal)},Int}, r::AbstractUnitRange)
    val = 1 + Int(eq.x - first(r))
    first(r) <= eq.x <= last(r) ? (val:val) : (1:0)
end

# See also findrange.jl for findindex(<=(3), range) things.

"""
    inferdim(key, axiskeys::Tuple)

When you call `A(key)` for `ndims(A) > 1`, this returns which `d` you meant,
if unambigous, by comparing types & gradually widening.
"""
@generated inferdim(arg, tup) = _inferdim(arg, map(eltype, Tuple(tup.parameters)))

function _inferdim(argT, types, subtypes=())
    types == subtypes && throw(ArgumentError("key of type $argT doesn't match any dimensions"))

    # First look for direct match
    ds = findall(T -> argT <: T, types)

    if length(ds) == 1
        return first(ds)
    elseif length(ds) >= 2
        throw(ArgumentError("key of type $argT is ambiguous, matches dimensions $(Tuple(ds))"))
    end

    # If no direct match, look for a container whose eltype matches:
    if argT <: Selector || argT <: AbstractArray || argT <: Interval
        ds = findall(T -> eltype(argT) <: T, types)
    elseif argT <: Base.Fix2 # Base.Fix2{typeof(==),Int64}
        ds = findall(T -> argT.parameters[2] <: T, types)
    end

    if length(ds) == 1
        return first(ds)
    elseif length(ds) >= 2
        throw(ArgumentError("key of type $argT is ambiguous, matches dimensions $(Tuple(ds))"))
    end

    # Otherwise, widen the key types and try again.
    # This will recurse until types stop changing.
    supers = map(T -> supertype(T) == Any ? T : supertype(T), types)
    return _inferdim(argT, supers, types)
end


text:
__precompile__()

module Gadfly

using Colors
using Compat
using Compose
using DataArrays
using DataFrames
using DataStructures
using JSON
using Showoff
using IndirectArrays
using CategoricalArrays

import IterTools
import IterTools: distinct, drop, chain
import Compose: draw, hstack, vstack, gridstack, isinstalled, parse_colorant
import Base: +, -, /, *,
             copy, push!, start, next, done, show, getindex, cat,
             show, isfinite, display
import Distributions: Distribution

export Plot, Layer, Theme, Col, Row, Scale, Coord, Geom, Guide, Stat, Shape, render, plot,
       style, layer, spy, set_default_plot_size, set_default_plot_format, prepare_display

@deprecate circle Shape.circle
@deprecate square Shape.square
@deprecate diamond Shape.diamond
@deprecate cross Shape.cross
@deprecate xcross Shape.xcross
@deprecate utriangle Shape.utriangle
@deprecate dtriangle Shape.dtriangle
@deprecate star1 Shape.star1
@deprecate star2 Shape.star2
@deprecate hexagon Shape.hexagon
@deprecate octogon Shape.octogon
@deprecate hline Shape.hline
@deprecate vline Shape.vline

# Re-export some essentials from Compose
export SVGJS, SVG, PGF, PNG, PS, PDF, draw, inch, mm, cm, px, pt, color, @colorant_str, vstack, hstack, title, gridstack

function __init__()
    # Define an XML namespace for custom attributes
    Compose.xmlns["gadfly"] = "http://www.gadflyjl.org/ns"
    if haskey(ENV, "GADFLY_THEME")
        theme = ENV["GADFLY_THEME"]
        try
            push_theme(Symbol(strip(theme)))
        catch err
            warn("Error loading Gadfly theme $theme (set by GADFLY_THEME env variable)")
            show(err)
        end
    else
        push_theme(Juno.isactive() ? :dark : :default)
    end
end

const ColorOrNothing = Union{Colorant, (Void)}

element_aesthetics(::Any) = []
input_aesthetics(::Any) = []
output_aesthetics(::Any) = []
default_scales(::Any) = []
default_scales(x::Any, t) = default_scales(x)
default_statistic(::Any) = Stat.identity()
element_coordinate_type(::Any) = Coord.cartesian

abstract type Element end
abstract type ScaleElement       <: Element end
abstract type CoordinateElement  <: Element end
abstract type GeometryElement    <: Element end
abstract type GuideElement       <: Element end
abstract type StatisticElement   <: Element end

include("misc.jl")
include("ticks.jl")
include("color_misc.jl")
include("varset.jl")
include("shapes.jl")
include("data.jl")
include("aesthetics.jl")
include("mapping.jl")
include("scale.jl")
include("theme.jl")
include("open_file.jl")

# The layer and plot functions can also take functions that are evaluated with
# no arguments and are expected to produce an element.
const ElementOrFunction{T <: Element} = Union{Element, Base.Callable, Theme}

const gadflyjs = joinpath(dirname(Base.source_path()), "gadfly.js")

# Set prefereed canvas size when rendering a plot without an explicit call to
# `draw`.
"""
```
set_default_plot_size(width::Compose.MeasureOrNumber, height::Compose.MeasureOrNumber)
```
Sets preferred canvas size when rendering a plot without an explicit call to draw
"""
set_default_plot_size(width::Compose.MeasureOrNumber, height::Compose.MeasureOrNumber) =
        Compose.set_default_graphic_size(width, height)

"""
```
set_default_plot_format(fmt::Symbol)
```
Sets the default plot format
"""
set_default_plot_format(fmt::Symbol) = Compose.set_default_graphic_format(fmt)

# A plot has zero or more layers. Layers have a particular geometry and their
# own data, which is inherited from the plot if not given.
mutable struct Layer <: Element
    data_source::Union{(Void), MeltedData, AbstractArray, AbstractDataFrame}
    mapping::Dict
    statistics::Vector{StatisticElement}
    geom::GeometryElement
    theme::Union{(Void), Theme}
    order::Int
end
Layer() = Layer(nothing, Dict(), StatisticElement[], Geom.nil(), nothing, 0)
Layer(l::Layer) = Layer(l.data_source, l.mapping, l.statistics, l.geom, l.theme, l.order)
copy(l::Layer) = Layer(l)

"""
```
layer(data_source::@compat(Union{AbstractDataFrame, (@compat Void)}),
               elements::ElementOrFunction...; mapping...)
```
Creates layers based on elements

### Args
* data_source: The data source as a dataframe
* elements: The elements
* mapping: mapping

### Returns
An array of layers
"""
function layer(data_source::Union{AbstractDataFrame, (Void)},
               elements::ElementOrFunction...; mapping...)
    mapping = Dict{Symbol, Any}(mapping)
    lyr = Layer()
    lyr.data_source = data_source
    lyr.mapping = cleanmapping(mapping)
    if haskey(mapping, :order)
        lyr.order = mapping[:order]
    end
    lyrs = Layer[lyr]
    for element in elements
        add_plot_element!(lyrs, element)
    end
    lyrs
end

layer(elements::ElementOrFunction...; mapping...) =
        layer(nothing, elements...; mapping...)

add_plot_element!(lyrs::Vector{Layer}, arg::T) where {T<:Element} =
        error("Layers can't be used with elements of type $(typeof(arg))")

add_plot_element!(lyrs::Vector{Layer}, arg::ScaleElement) =
        error("Scales cannot be passed to layers, they must be specified at the plot level.")

function add_plot_element!(lyrs::Vector{Layer}, arg::GeometryElement)
    if lyrs[end].geom !== Geom.nil()
        push!(lyrs, copy(lyrs[end]))
    end
    lyrs[end].geom = arg
end

add_plot_element!(lyrs::Vector{Layer}, arg::Base.Callable) = add_plot_element!(lyrs, arg())

function add_plot_element!(lyrs::Vector{Layer}, arg::StatisticElement)
    for lyr in lyrs
        push!(lyr.statistics, arg)
    end
end

add_plot_element!(lyrs::Vector{Layer}, arg::Theme) = [lyr.theme = arg for lyr in lyrs]

# A full plot specification.
mutable struct Plot
    layers::Vector{Layer}
    data_source::Union{(Void), MeltedData, AbstractArray, AbstractDataFrame}
    data::Data
    scales::Vector{ScaleElement}
    statistics::Vector{StatisticElement}
    coord::Union{(Void), CoordinateElement}
    guides::Vector{GuideElement}
    theme::Theme
    mapping::Dict
end
Plot() = Plot(Layer[], nothing, Data(), ScaleElement[], StatisticElement[],
      nothing, GuideElement[], current_theme(), Dict())

layers(p::Plot) = p.layers

function add_plot_element!(p::Plot, arg::GeometryElement)
    if !isempty(p.layers) && isa(p.layers[end].geom, Geom.Nil)
        p.layers[end].geom = arg
    else
        layer = Layer()
        layer.geom = arg
        push!(p.layers, layer)
    end
end

function add_plot_element!(p::Plot, arg::StatisticElement)
    if isempty(p.layers)
        push!(p.layers, Layer())
    end

    push!(p.layers[end].statistics, arg)
end

add_plot_element!(p::Plot, arg::Function) = add_plot_element!(p, arg())
add_plot_element!(p::Plot, arg::ScaleElement) = push!(p.scales, arg)
add_plot_element!(p::Plot, arg::CoordinateElement) = p.coord = arg
add_plot_element!(p::Plot, arg::GuideElement) = push!(p.guides, arg)
add_plot_element!(p::Plot, arg::Layer) = push!(p.layers, arg)
add_plot_element!(p::Plot, arg::Vector{Layer}) = append!(p.layers, arg)
add_plot_element!(p::Plot, f::Type{T}) where {T <: Element} = add_plot_element!(p, f())
add_plot_element!(p::Plot, theme::Theme) = p.theme = theme

# Create a new plot.
#
# Grammar of graphics style plotting consists of specifying a dataset, one or
# more plot elements (scales, coordinates, geometries, etc), and binding of
# aesthetics to columns or expressions of the dataset.
#
# For example, a simple scatter plot would look something like:
#
#     plot(my_data, Geom.point, x="time", y="price")
#
# Where "time" and "price" are the names of columns in my_data.
#
# Args:
#   data_source: Data to be bound to aesthetics.
#   mapping: Aesthetics symbols (e.g. :x, :y, :color) mapped to
#            names of columns in the data frame or other expressions.
#   elements: Geometries, statistics, etc.

# because a call to layer() expands to a vector of layers (one for each Geom
# supplied), we need to allow Vector{Layer} to count as an Element for the
# purposes of plot().
const ElementOrFunctionOrLayers = Union{ElementOrFunction, Vector{Layer}}

"""
```
    function plot(data_source::@compat(Union{AbstractMatrix, AbstractDataFrame}),
              elements::ElementOrFunctionOrLayers...; mapping...)
```

Create a new plot.

Grammar of graphics style plotting consists of specifying a dataset, one or
more plot elements (scales, coordinates, geometries, etc), and binding of
aesthetics to columns or expressions of the dataset.

For example, a simple scatter plot would look something like:

plot(my_data, Geom.point, x="time", y="price")

Where "time" and "price" are the names of columns in my_data.

### Args:
* data_source: Data to be bound to aesthetics.
* elements: Geometries, statistics, etc.
* mapping: Aesthetics symbols (e.g. :x, :y, :color) mapped to names of columns in the data frame or other expressions.
"""
function plot(data_source::Union{AbstractArray, AbstractDataFrame},
              elements::ElementOrFunctionOrLayers...; mapping...)
    mappingdict = Dict{Symbol, Any}(mapping)
    return plot(data_source, mappingdict, elements...)
end

function plot(elements::ElementOrFunctionOrLayers...; mapping...)
    mappingdict = Dict{Symbol, Any}(mapping)
    plot(nothing, mappingdict, elements...)
end

# The old fashioned (pre named arguments) version of plot.
#
# This version takes an explicit mapping dictionary, mapping aesthetics symbols
# to expressions or columns in the data frame.
#
# Args:
#   data_source: Data to be bound to aesthetics.
#   mapping: Dictionary of aesthetics symbols (e.g. :x, :y, :color) to
#            names of columns in the data frame or other expressions.
#   elements: Geometries, statistics, etc.
#
# Returns:
#   A Plot object.
#
"""
```
function plot(data_source::@compat(Union{(@compat Void), AbstractMatrix, AbstractDataFrame}),
              mapping::Dict, elements::ElementOrFunctionOrLayers...)
```
The old fashioned (pre named arguments) version of plot.

This version takes an explicit mapping dictionary, mapping aesthetics symbols
to expressions or columns in the data frame.

### Args:
* data_source: Data to be bound to aesthetics.
* mapping: Dictionary of aesthetics symbols (e.g. :x, :y, :color) to
            names of columns in the data frame or other expressions.
*   elements: Geometries, statistics, etc.

### Returns:
A Plot object.
"""
function plot(data_source::Union{(Void), AbstractArray, AbstractDataFrame},
              mapping::Dict, elements::ElementOrFunctionOrLayers...)
    mapping = cleanmapping(mapping)
    p = Plot()
    for element in elements
        add_plot_element!(p, element)
    end

    p.data_source = evalmapping!(mapping, data_source, p.data)
    p.mapping = mapping

    return p
end

include("poetry.jl")

function Base.push!(p::Plot, element::ElementOrFunctionOrLayers)
    add_plot_element!(p, element)
    return p
end

# Turn a graph specification into a graphic.
#
# This is where magic happens (sausage is made). Processing all the parts of the
# plot is actually pretty simple. It's made complicated by trying to handle
# defaults. With that aside, plots are made in the following steps.
#
#    I. Apply scales to transform raw data to the form expected by the aesthetic.
#   II. Apply statistics to the scaled data. Statistics are essentially functions
#       that map one or more aesthetics to one or more aesthetics.
#  III. Apply coordinates. Currently all this does is figure out the coordinate
#       system used by the plot panel canvas.
#   IV. Render geometries. This gives us one or more compose forms suitable to be
#       composed with the plot's panel.
#    V. Render guides. Guides are conceptually very similar to geometries but with
#       the ability to be placed outside of the plot panel.
#
#  Finally there is a very important call to layout_guides which puts everything
#  together.
#
# Args:
#   plot: a plot to render.
#
# Returns:
#   A compose Canvas containing the graphic.
#
function render_prepare(plot::Plot)
    if isempty(plot.layers)
        layer = Layer()
        layer.geom = Geom.point()
        push!(plot.layers, layer)
    end

    # TODO: When subplots are given in multiple layers, we should rearrange,
    # putting the layers in one subplot instead.
    if sum([isa(layer.geom, Geom.SubplotGeometry) for layer in plot.layers]) > 1
        error("""
              Subplot geometries can not be used in multiple layers. Instead
              use layers within one subplot geometry.
              """)
    end

    # Process layers, filling inheriting mappings or data from the Plot where
    # they are missing.
    datas = Array{Data}(length(plot.layers))
    for (i, layer) in enumerate(plot.layers)
        if layer.data_source === nothing && isempty(layer.mapping)
            layer.data_source = plot.data_source
            layer.mapping = plot.mapping
            datas[i] = plot.data
        else
            datas[i] = Data()

            if layer.data_source === nothing
                layer.data_source = plot.data_source
            end

            if isempty(layer.mapping)
                layer.mapping = plot.mapping
            end

            evalmapping!(layer.mapping, layer.data_source, datas[i])
        end
        if isa(layer.geom, Geom.Nil); layer.geom = Geom.point(); end # see #1062
    end

    # We need to process subplot layers somewhat as though they were regular
    # plot layers. This is the only way scales, etc, can be consistently
    # applied.
    subplot_datas = Data[]
    for (layer, layer_data) in zip(plot.layers, datas)
        if isa(layer.geom, Geom.SubplotGeometry)
            for subplot_layer in layers(layer.geom)
                subplot_data = Data()
                if subplot_layer.data_source === nothing
                    subplot_layer.data_source = layer.data_source
                end

                if isempty(subplot_layer.mapping)
                    subplot_layer.mapping = layer.mapping
                end

                evalmapping!(subplot_layer.mapping, subplot_layer.data_source, subplot_data)
                push!(subplot_datas, subplot_data)
            end
        end
    end

    # Figure out the coordinates
    coord = plot.coord
    for layer in plot.layers
        coord_type = element_coordinate_type(layer.geom)
        if coord === nothing
            coord = coord_type()
        elseif typeof(coord) != coord_type
            error("Plot uses multiple coordinates: $(typeof(coord)) and $(coord_type)")
        end
    end

    # Add default statistics for geometries.
    layer_stats = Array{Vector{StatisticElement}}(length(plot.layers))
    for (i, layer) in enumerate(plot.layers)
        layer_stats[i] = isempty(layer.statistics) ? ( isa(layer.geom, Geom.SubplotGeometry) ?
                default_statistic(layer.geom) : [default_statistic(layer.geom)] ) : layer.statistics
    end

    # auto-enumeration: add Stat.x/y_enumerate when x and y is needed but only
    # one defined
    for (i, layer) in enumerate(plot.layers)
        layer_needed_aes = element_aesthetics(layer.geom)
        layer_defined_aes = Set{Symbol}()
        union!(layer_defined_aes, keys(layer.mapping))
        for stat in layer.statistics
            union!(layer_defined_aes, output_aesthetics(stat))
        end

        if mapreduce(x->in(x,layer_needed_aes),|,[:x,:xmax,:xmin]) &&
                mapreduce(y->in(y,layer_needed_aes),|,[:y,:ymax,:ymin])
            if !mapreduce(x->in(x,layer_defined_aes),|,[:x,:xmax,:xmin])
                unshift!(layer_stats[i], Stat.x_enumerate)
            elseif !mapreduce(y->in(y,layer_defined_aes),|,[:y,:ymax,:ymin])
                unshift!(layer_stats[i], Stat.y_enumerate)
            end
        end
    end

    used_aesthetics = Set{Symbol}()
    for layer in plot.layers
        union!(used_aesthetics, element_aesthetics(layer.geom))
    end

    for stats in layer_stats, stat in stats
        union!(used_aesthetics, input_aesthetics(stat))
    end

    mapped_aesthetics = Set(keys(plot.mapping))
    for layer in plot.layers
        union!(mapped_aesthetics, keys(layer.mapping))
    end

    defined_unused_aesthetics = setdiff(mapped_aesthetics, used_aesthetics)
    isempty(defined_unused_aesthetics) ||
            warn("The following aesthetics are mapped, but not used by any geometry:\n",
                join([string(a) for a in defined_unused_aesthetics], ", "))

    scaled_aesthetics = Set{Symbol}()
    for scale in plot.scales
        union!(scaled_aesthetics, element_aesthetics(scale))
    end

    # Only one scale can be applied to an aesthetic (without getting some weird
    # and incorrect results), so we organize scales into a dict.
    scales = Dict{Symbol, ScaleElement}()
    for scale in plot.scales
        for var in element_aesthetics(scale)
            scales[var] = scale
        end
    end

    unscaled_aesthetics = setdiff(used_aesthetics, scaled_aesthetics)

    _theme(plt, lyr) = lyr.theme == nothing ? plt.theme : lyr.theme

    # Add default scales for statistics.
    layer_stats_with_theme = map(plot.layers, layer_stats) do l, stats
        map(s->(s, _theme(l, plot)), collect(stats))
    end

    for element in chain([(s, plot.theme) for s in plot.statistics],
                         [(l.geom, _theme(plot, l)) for l in plot.layers],
                         layer_stats_with_theme...)

        for scale in default_scales(element...)
            # Use the statistics default scale only when it covers some
            # aesthetic that is not already scaled.
            scale_aes = Set(element_aesthetics(scale))
            if !isempty(intersect(scale_aes, unscaled_aesthetics))
                for var in scale_aes
                    scales[var] = scale
                end
                setdiff!(unscaled_aesthetics, scale_aes)
            end
        end
    end

    # Assign scales to mapped aesthetics first.
    for var in unscaled_aesthetics
        in(var, mapped_aesthetics) || continue

        var_data = getfield(plot.data, var)
        if var_data == nothing
            for data in datas
                var_layer_data = getfield(data, var)
                if var_layer_data != nothing
                    var_data = var_layer_data
                    break
                end
            end
        end

        var_data == nothing && continue

        t = classify_data(var_data)
        if scale_exists(t, var)
            scale = get_scale(t, var, plot.theme)
            scale_aes = Set(element_aesthetics(scale))
            for var in scale_aes
                scales[var] = scale
            end
        end
    end

    for var in unscaled_aesthetics
        (haskey(plot.mapping, var) || haskey(scales, var)) && continue

        t = :categorical
        for data in chain(datas, subplot_datas)
            val = getfield(data, var)
            if val != nothing && val != :categorical
                t = classify_data(val)
            end
        end

        if scale_exists(t, var)
            scale = get_scale(t, var, plot.theme)
            scale_aes = Set(element_aesthetics(scale))
            for var in scale_aes
                scales[var] = scale
            end
        end
    end

    # Avoid clobbering user-defined guides with default guides (e.g.
    # in the case of labels.)
    guides = copy(plot.guides)
    explicit_guide_types = Set()
    for guide in guides
        push!(explicit_guide_types, typeof(guide))
    end

    statistics = Set{StatisticElement}()
    for statistic in plot.statistics
        push!(statistics, statistic)
    end

    # Default guides and statistics
    facet_plot = true
    for layer in plot.layers
        if typeof(layer.geom) != Geom.subplot_grid
            facet_plot = false
            break
        end
    end

    if !facet_plot
        in(Guide.PanelBackground, explicit_guide_types) || push!(guides, Guide.background())
        in(Guide.QuestionMark, explicit_guide_types) || push!(guides, Guide.questionmark())
        in(Guide.HelpScreen, explicit_guide_types) || push!(guides, Guide.helpscreen())
        in(Guide.CrossHair, explicit_guide_types) || push!(guides, Guide.crosshair())
        in(Guide.XTicks, explicit_guide_types) || push!(guides, Guide.xticks())
        in(Guide.YTicks, explicit_guide_types) || push!(guides, Guide.yticks())
    end

    for guide in guides
        push!(statistics, default_statistic(guide))
    end

    mapped_and_used = function(vs)
        any(Bool[in(v, mapped_aesthetics) && in(v, used_aesthetics) for v in vs])
    end

    choose_name = function(vs, fallback)
        for v in vs
            if haskey(plot.data.titles, v)
                return plot.data.titles[v]
            end
        end

        for v in vs
            for data in datas
                if haskey(data.titles, v)
                    return data.titles[v]
                end
            end
        end

        fallback
    end

    if mapped_and_used(x_axis_label_aesthetics) &&
        !in(Guide.XLabel, explicit_guide_types)
        label = choose_name(x_axis_label_aesthetics, "x")
        if facet_plot && haskey(plot.data.titles, :xgroup)
            label = string(label, " <i><b>by</b></i> ", plot.data.titles[:xgroup])
        end

        push!(guides, Guide.xlabel(label))
    end

    if mapped_and_used(y_axis_label_aesthetics) &&
       !in(Guide.YLabel, explicit_guide_types)
        label = choose_name(y_axis_label_aesthetics, "y")
        if facet_plot && haskey(plot.data.titles, :ygroup)
            label = string(label, " <i><b>by</b></i> ", plot.data.titles[:ygroup])
        end

        push!(guides, Guide.ylabel(label))
    end

    # I. Scales
    layer_aess = Scale.apply_scales(IterTools.distinct(values(scales)),
                                    datas..., subplot_datas...)

    # set default labels
    for (i, layer) in enumerate(plot.layers)
        if layer_aess[i].color_key_title == nothing &&
           haskey(layer.mapping, :color) &&
           !isa(layer.mapping[:color], AbstractArray)
           layer_aess[i].color_key_title = string(layer.mapping[:color])
       end
    end

    if layer_aess[1].color_key_title == nothing &&
       haskey(plot.mapping, :color) && !isa(plot.mapping[:color], AbstractArray)
        layer_aess[1].color_key_title = string(plot.mapping[:color])
    end

    # IIa. Layer-wise statistics
    for (stats, aes) in zip(layer_stats, layer_aess)
        Stat.apply_statistics(stats, scales, coord, aes)
    end

    # IIb. Plot-wise Statistics
    plot_aes = concat(layer_aess...)
    statistics = collect(statistics)
    Stat.apply_statistics(statistics, scales, coord, plot_aes)

    # Add some default guides determined by defined aesthetics
    supress_colorkey = false
    for layer in plot.layers
        if isa(layer.geom, Geom.SubplotGeometry) &&
                haskey(layer.geom.guides, Guide.ColorKey)
            supress_colorkey = true
            break
        end
    end

    if !supress_colorkey &&
       !all([aes.color === nothing for aes in [plot_aes, layer_aess...]]) &&
       !in(Guide.ColorKey, explicit_guide_types) &&
       !in(Guide.ManualColorKey, explicit_guide_types)
        push!(guides, Guide.colorkey())
    end

    # build arrays of scaled aesthetics for layers within subplots
    layer_subplot_aess = Array{Vector{Aesthetics}}(length(plot.layers))
    layer_subplot_datas = Array{Vector{Data}}(length(plot.layers))
    j = 1
    for (i, layer) in enumerate(plot.layers)
        layer_subplot_aess[i] = Aesthetics[]
        layer_subplot_datas[i] = Data[]
        if isa(layer.geom, Geom.SubplotGeometry)
            for subplot_layer in layers(layer.geom)
                push!(layer_subplot_aess[i], layer_aess[length(datas) + j])
                push!(layer_subplot_datas[i], subplot_datas[j])
                j += 1
            end
        end
    end

    return (plot, coord, plot_aes,
            layer_aess, layer_stats, layer_subplot_aess, layer_subplot_datas,
            scales, guides)
end

"""
```
render(plot::Plot)
```
Render a plot based on the Plot object

### Args
*   plot: Plot to be rendered.

### Returns
A Compose context containing the rendered plot.
"""
function render(plot::Plot)
    (plot, coord, plot_aes,
     layer_aess, layer_stats, layer_subplot_aess, layer_subplot_datas,
     scales, guides) = render_prepare(plot)

    root_context = render_prepared(plot, coord, plot_aes, layer_aess,
                                   layer_stats, layer_subplot_aess,
                                   layer_subplot_datas,
                                   scales, guides)

    ctx =  pad_inner(root_context, plot.theme.plot_padding...)

    if plot.theme.background_color != nothing
        compose!(ctx, (context(order=-1000000),
                        fill(plot.theme.background_color),
                        stroke(nothing), rectangle()))
    end

    return ctx
end

# Render a plot given a precomputed Aesthetics object for each layer.
#
# Additionally, without all the work to choose reasonable defaults performed by
# `render`. This is a separate function from `render` to facilitate rendering
# subplots.
#
# Args:
#   plot: Plot to be rendered.
#   aess: A vector of precomputed Aesthetics objects of the same length
#       as plot.layers.
#   layer_stats: A vector of statistic elements of the same length as
#       plot.layers.
#   layer_subplot_aesthetics: An array of aesthetics for each layer in the plot.
#       If the layer in a subplot geometry, the array is scaled data for each
#       sub-geometry, otherwise it's empty. I just melted your brain, didn't I?
#   scales: Dictionary mapping an aesthetics symbol to the scale applied to it.
#   statistics: Statistic elements applied plot-wise.
#   guides: Guide elements indexed by type. (Only one type of each guide may
#       be in the same plot.)
#   preserve_plot_context_size: Don't squish the plot to fit the guides.
#       Guides will be drawn outside the context
#
# Returns:
#   A Compose context containing the rendered plot.
#
function render_prepared(plot::Plot,
                         coord::CoordinateElement,
                         plot_aes::Aesthetics,
                         layer_aess::Vector{Aesthetics},
                         layer_stats::Vector{Vector{StatisticElement}},
                         layer_subplot_aess::Vector{Vector{Aesthetics}},
                         layer_subplot_datas::Vector{Vector{Data}},
                         scales::Dict{Symbol, ScaleElement},
                         guides::Vector{GuideElement};
                         table_only=false,
                         dynamic=true)
    # III. Coordinates
    plot_context = Coord.apply_coordinate(coord, vcat(plot_aes,
                                          layer_aess), scales)

    # IV. Geometries
    themes = Theme[layer.theme === nothing ? plot.theme : layer.theme
                   for layer in plot.layers]
    zips = trim_zip(plot.layers, layer_aess,
                                                   layer_subplot_aess,
                                                   layer_subplot_datas,
               themes)

    compose!(plot_context,
             [compose(context(order=layer.order), render(layer.geom, theme, aes,
                                                         subplot_aes, subplot_data,
                                                         scales))
              for (layer, aes, subplot_aes, subplot_data, theme) in zips]...)

    # V. Guides
    guide_contexts = Any[]
    for guide in guides
        guide_context = render(guide, plot.theme, plot_aes, dynamic)
        if guide_context != nothing
            append!(guide_contexts, guide_context)
        end
    end

    tbl = Guide.layout_guides(plot_context, coord,
                              plot.theme, guide_contexts...)
    if table_only
        return tbl
    end

    c = compose!(context(), tbl)
    class = "plotroot"
    if haskey(scales, :x) && isa(scales[:x], Scale.ContinuousScale) && scales[:x].scalable
        class = string(class, " xscalable")
    end
    if haskey(scales, :y) && isa(scales[:y], Scale.ContinuousScale) && scales[:y].scalable
        class = string(class, " yscalable")
    end

    compose(c, svgclass(class), jsinclude(gadflyjs, ("Gadfly", "Gadfly")))
end

# A convenience version of Compose.draw that let's you skip the call to render.
"""
```
draw(backend::Compose.Backend, p::Plot)
```
A convenience version of Compose.draw without having to call render

### Args
* backend: The Compose.Backend object
* p: The Plot object
"""
draw(backend::Compose.Backend, p::Plot) = draw(backend, render(p))

"""
    title(ctx::Context, str::String, props::Property...) -> Context

Add a title string to a group of plots, typically created with `vstack`, `hstack`, or `gridstack`.

# Examples

```
p1 = plot(x=[1,2], y=[3,4], Geom.line);
p2 = plot(x=[1,2], y=[4,3], Geom.line);
title(hstack(p1,p2), "my latest data", Compose.fontsize(18pt), fill(colorant"red"))
```
"""
title(ctx::Context, str::String, props::Compose.Property...) = vstack(
    compose(context(0, 0, 1, 0.1), text(0.5, 1.0, str, hcenter, vbottom), props...),
    compose(context(0, 0, 1, 0.9), ctx))

# Convenience stacking functions
"""
    vstack(ps::Union{Plot,Context}...)
    vstack(ps::Vector)

Arrange plots into a vertical column.  Use `context()` as a placeholder for an empty panel.  Heterogeneous vectors must be typed.  See also `hstack`, `gridstack`, `subplot_grid`.

# Examples

```
p1 = plot(x=[1,2], y=[3,4], Geom.line);
p2 = Compose.context();
vstack(p1, p2)
vstack(Union{Plot,Compose.Context}[p1, p2])
"""
vstack(ps::Union{Plot,Context}...) = vstack([typeof(p)==Plot ? render(p) : p for p in ps]...)
vstack(ps::Vector{Plot}) = vstack(ps...)
vstack(ps::Vector{Union{Plot,Context}}) = vstack(ps...)

"""
    hstack(ps::Union{Plot,Context}...)
    hstack(ps::Vector)

Arrange plots into a horizontal row.  Use `context()` as a placeholder for an empty panel.  Heterogeneous vectors must be typed.  See also `vstack`, `gridstack`, `subplot_grid`.

# Examples

```
p1 = plot(x=[1,2], y=[3,4], Geom.line);
p2 = Compose.context();
hstack(p1, p2)
hstack(Union{Plot,Compose.Context}[p1, p2])
"""
hstack(ps::Union{Plot,Context}...) = hstack([typeof(p)==Plot ? render(p) : p for p in ps]...)
hstack(ps::Vector{Plot}) = hstack(ps...)
hstack(ps::Vector{Union{Plot,Context}}) = hstack(ps...)

"""
    gridstack(ps::Matrix{Union{Plot,Context}})

Arrange plots into a rectangular array.  Use `context()` as a placeholder for an empty panel.  Heterogeneous matrices must be typed.  See also `hstack`, `vstack`.

# Examples

```
p1 = plot(x=[1,2], y=[3,4], Geom.line);
p2 = Compose.context();
gridstack([p1 p1; p1 p1])
gridstack(Union{Plot,Compose.Context}[p1 p2; p2 p1])
```
"""
_gridstack(ps::Matrix) = gridstack(map(p->typeof(p)==Plot ? render(p) : p, ps))
gridstack(ps::Matrix{Plot}) = _gridstack(ps)
gridstack(ps::Matrix{Union{Plot,Context}}) = _gridstack(ps)

# show functions for all supported compose backends.

function show(io::IO, m::MIME"text/html", p::Plot)
    buf = IOBuffer()
    svg = SVGJS(buf, Compose.default_graphic_width,
                Compose.default_graphic_height, false)
    draw(svg, p)
    show(io, m, svg)
end

function show(io::IO, m::MIME"image/svg+xml", p::Plot)
    buf = IOBuffer()
    svg = SVG(buf, Compose.default_graphic_width,
              Compose.default_graphic_height, false)
    draw(svg, p)
    show(io, m, svg)
end

try
    getfield(Compose, :Cairo) # throws if Cairo isn't being used
    global show
    function show(io::IO, ::MIME"image/png", p::Plot)
        draw(PNG(io, Compose.default_graphic_width,
                 Compose.default_graphic_height), p)
    end
end

try
    getfield(Compose, :Cairo) # throws if Cairo isn't being used
    global show
    function show(io::IO, ::MIME"application/postscript", p::Plot)
        draw(PS(io, Compose.default_graphic_width,
             Compose.default_graphic_height), p);
    end
end

function show(io::IO, ::MIME"text/plain", p::Plot)
    write(io, "Plot(...)")
end

function default_mime()
    if Compose.default_graphic_format == :png
        "image/png"
    elseif Compose.default_graphic_format == :svg
        "image/svg+xml"
    elseif Compose.default_graphic_format == :html
        "text/html"
    elseif Compose.default_graphic_format == :ps
        "application/postscript"
    elseif Compose.default_graphic_format == :pdf
        "application/pdf"
    else
        ""
    end
end

import Base.REPL: REPLDisplay

"""
    display(p::Plot)

Render `p` to a multimedia display, typically an internet browser.

This function is handy when rendering by `plot` has been suppressed
with either trailing semi-colon or by calling it within a function.
"""
function display(d::REPLDisplay, p::Union{Plot,Compose.Context})
    if mimewritable("text/html", p)
        display(d,"text/html", p)
        return
    elseif mimewritable("image/png", p)
        display(d,"image/png", p)
        return
    elseif mimewritable("application/pdf", p)
        display(d,"application/pdf", p)
        return
    elseif mimewritable("image/svg+xml", p)
        display(d,"image/svg+xml", p)
        return
    elseif mimewritable("application/postscript", p)
        display(d,"application/postscript", p)
        return
    end
    throw(MethodError)
end

# Fallback display method. When there isn't a better option, we write to a
# temporary file and try to open it.
function display(d::REPLDisplay, ::MIME"image/png", p::Union{Plot,Compose.Context})
    filename = string(tempname(), ".png")
    output = open(filename, "w")
    draw(PNG(output, Compose.default_graphic_width,
             Compose.default_graphic_height), p)
    close(output)
    open_file(filename)
end

function display(d::REPLDisplay, ::MIME"image/svg+xml", p::Union{Plot,Compose.Context})
    filename = string(tempname(), ".svg")
    output = open(filename, "w")
    draw(SVG(output, Compose.default_graphic_width,
             Compose.default_graphic_height), p)
    close(output)
    open_file(filename)
end

function display(d::REPLDisplay, ::MIME"text/html", p::Union{Plot,Compose.Context})
    filename = string(tempname(), ".html")
    output = open(filename, "w")

    plot_output = IOBuffer()
    draw(SVGJS(plot_output, Compose.default_graphic_width,
               Compose.default_graphic_height, false), p)
    plotsvg = String(take!(plot_output))

    write(output,
        """
        <!DOCTYPE html>
        <html>
          <head>
            <title>Gadfly Plot</title>
            <meta charset="utf-8">
          </head>
            <body style="margin:0">
            <script charset="utf-8">
                $(readstring(Compose.snapsvgjs))
            </script>
            <script charset="utf-8">
                $(readstring(gadflyjs))
            </script>

            $(plotsvg)
          </body>
        </html>
        """)
    close(output)
    open_file(filename)
end

function display(d::REPLDisplay, ::MIME"application/postscript", p::Union{Plot,Compose.Context})
    filename = string(tempname(), ".ps")
    output = open(filename, "w")
    draw(PS(output, Compose.default_graphic_width,
            Compose.default_graphic_height), p)
    close(output)
    open_file(filename)
end

function display(d::REPLDisplay, ::MIME"application/pdf", p::Union{Plot,Compose.Context})
    filename = string(tempname(), ".pdf")
    output = open(filename, "w")
    draw(PDF(output, Compose.default_graphic_width,
             Compose.default_graphic_height), p)
    close(output)
    open_file(filename)
end

# Display in Juno

import Juno: Juno, @render, media, Media

media(Plot, Media.Plot)

@render Juno.PlotPane p::Plot begin
    x, y = Juno.plotsize()
    set_default_plot_size(x*Gadfly.px, y*Gadfly.px)
    HTML(stringmime("text/html", p))
end

@render Juno.Editor p::Gadfly.Plot begin
    Juno.icon("graph")
end

include("coord.jl")
include("geometry.jl")
include("guide.jl")
include("statistics.jl")

# All aesthetics must have a scale. If none is given, we use a default.
# The default depends on whether the input is discrete or continuous (i.e.,
# PooledDataVector or DataVector, respectively).
const default_aes_scales = Dict{Symbol, Dict}(

    :distribution => Dict{Symbol, Any}(
        :x => Scale.x_distribution(),
        :y => Scale.y_distribution()
    ),

    :functional => Dict{Symbol, Any}(
        :z      => Scale.z_func(),
        :y      => Scale.y_func(),
        :shape  => Scale.shape_identity(),
        :size   => Scale.size_identity(),
        :color  => Scale.color_identity(),
    ),

    :numerical => Dict{Symbol, Any}(
        :x           => Scale.x_continuous(),
        :xmin        => Scale.x_continuous(),
        :xmax        => Scale.x_continuous(),
        :xintercept  => Scale.x_continuous(),
        :xslope      => Scale.x_continuous(),
        :xend        => Scale.x_continuous(),
        :yend        => Scale.y_continuous(),
        :y           => Scale.y_continuous(),
        :ymin        => Scale.y_continuous(),
        :ymax        => Scale.y_continuous(),
        :yintercept  => Scale.y_continuous(),
        :yslope      => Scale.y_continuous(),
        :middle      => Scale.y_continuous(),
        :upper_fence => Scale.y_continuous(),
        :lower_fence => Scale.y_continuous(),
        :upper_hinge => Scale.y_continuous(),
        :lower_hinge => Scale.y_continuous(),
        :xgroup      => Scale.xgroup(),
        :ygroup      => Scale.ygroup(),
        :shape       => Scale.shape_discrete(),
        :size        => Scale.size_continuous(),
        :group       => Scale.group_discrete(),
        :label       => Scale.label()
    ),

    :categorical => Dict{Symbol, Any}(
        :x          => Scale.x_discrete(),
        :xmin       => Scale.x_discrete(),
        :xmax       => Scale.x_discrete(),
        :xintercept => Scale.x_discrete(),
        :xend       => Scale.x_discrete(),
        :yend       => Scale.y_discrete(),
        :y          => Scale.y_discrete(),
        :ymin       => Scale.y_discrete(),
        :ymax       => Scale.y_discrete(),
        :yintercept => Scale.y_discrete(),
        :xgroup     => Scale.xgroup(),
        :ygroup     => Scale.ygroup(),
        :shape      => Scale.shape_discrete(),
        :size       => Scale.size_discrete(),
        :group      => Scale.group_discrete(),
        :label      => Scale.label()
    )
)

get_scale(::Val{t}, ::Val{var}, theme::Theme) where {t,var} = default_aes_scales[t][var]
get_scale(t::Symbol, var::Symbol, theme::Theme) = get_scale(Val{t}(), Val{var}(), theme)

### Override default getters for color scales
get_scale(::Val{:categorical}, ::Val{:color}, theme::Theme=current_theme()) =
        theme.discrete_color_scale
get_scale(::Val{:numerical}, ::Val{:color}, theme::Theme=current_theme()) =
        theme.continuous_color_scale

function scale_exists(t::Symbol, var::Symbol)
    if !haskey(default_aes_scales, t) || !haskey(default_aes_scales[t], var)
        method = methods(get_scale, (Val{t}, Val{var}, Theme))
        catchall = methods(get_scale, (Val{1}, Val{nothing}, Theme))
        first(method) !== first(catchall)
    else
        true
    end
end

# Determine whether the input is categorical or numerical

const CategoricalType = Union{AbstractString, Bool, Symbol}

classify_data(data::AbstractArray{T, N}) where {N, T <: Union{CategoricalType,Missing}}        = :categorical
classify_data(data::AbstractArray{T, N}) where {N, T <: Union{Base.Callable,Measure,Colorant}} = :functional
classify_data(data::CategoricalArray) = :categorical
classify_data(data::T) where {T <: Base.Callable} = :functional
classify_data(data::AbstractArray) = :numerical
classify_data(data::Distribution) = :distribution

function classify_data(data::AbstractArray{Any})
    for val in data
        if isa(val, CategoricalType)
            return :categorical
        end
    end
    :numerical
end

# Axis labels are taken whatever is mapped to these aesthetics, in order of
# preference.
const x_axis_label_aesthetics = [:x, :xmin, :xmax]
const y_axis_label_aesthetics = [:y, :ymin, :ymax]

end # module Gadfly


text:
using DrWatson
@quickactivate "StatReth"

# %%
using CSV
using DataFrames
using Turing
using Plots

include(srcdir("quap.jl"))
include(srcdir("tools.jl"))

# %%
# This is from the chapter
d = DataFrame(CSV.File(datadir("exp_raw/Howell_1.csv")))
sort!(d, :weight)  # to make the plots look nice

z_weight, unz_weight = zscore_transform(d.weight)
d.weight_s = z_weight(d.weight)

f_parabola(weight_s, a, b1, b2) = a + b1 * weight_s + b2 * weight_s^2

@model function parabola(weight_s, height)
    a ~ Normal(178, 20)
    b1 ~ LogNormal(0, 1)
    b2 ~ Normal(0, 1)
    σ ~ Uniform(0, 50)
    μ = f_parabola.(weight_s, a, b1, b2)
    height ~ MvNormal(μ, σ)
end

mparab = parabola(d.weight_s, d.height)
prior = sample(mparab, Prior(), 1_000) |> DataFrame

# %%
plot(legend = false)
for c in eachrow(prior[1:100, :])
    p = f_parabola.(d.weight_s, c.a, c.b1, c.b2)
    plot!(d.weight_s, p, color = :black, alpha = 0.2)
end
plot!()

# %%
"""
Thinking about the cbrt explanation, b2 should be negative so I get a right turn here as
well. Since I don't want the extremum to be in the range of my data I need the maximum to
be at like 3 or 4 or something like that. Over the range of ~4 stddev I want the height to
change 100 cm or so.
f(x) = b₂ * x² + b₁ * x + a
f'(x) = 2b₂ * x + b₁

f(2) - f(-2) = 100
f'(3) = 0
b₂ < 0

f(2) - f(-2) = (b₂ * 4 + b₁ * 2 + a) - (b₂ * 4 + b₁ * -2 + a)
             = 4b₁ = 100
            =>  b₁ =  25
f'(3) = 2b₂ * 3 + b₁ = 0 => -b₁ = 6b₂
     => b₂ = -25 / 6 ≈ -4

I'll move a by hand until it seams reasonable.
"""

@model function parabola(weight_s, height)
    a ~ Normal(130, 20)
    b1 ~ Normal(25, 5)
    b2 ~ truncated(Normal(-4, 2), -Inf, 0)
    σ ~ Exponential(5)
    μ = f_parabola.(weight_s, a, b1, b2)
    height ~ MvNormal(μ, σ)
end

mparab = parabola(d.weight_s, d.height)
prior = sample(mparab, Prior(), 100) |> DataFrame

plot(legend = false)
for c in eachrow(prior)
    p = f_parabola.(d.weight_s, c.a, c.b1, c.b2)
    plot!(d.weight_s, p, color = :black, alpha = 0.2)
end
plot!(ylims = (0, 200))


text:
function main(inputmatrix)

  println("Initial matrix looks like this:")
  println(inputmatrix)

  # an array of tuples representing the positions
  # dummy initialization so that it's never empty
  top3positions = [(1,1),(1,2),(1,3)]

  numrows = size(inputmatrix)[1]
  numcols = size(inputmatrix)[2]

  for i=1:numrows
    for j=1:numcols
      elem = (i,j)
      top3positions = updatepositions(inputmatrix,top3positions,elem)
    end
  end

  position1 = top3positions[1]
  position2 = top3positions[2]
  position3 = top3positions[3]

  outputmatrix = deepcopy(inputmatrix)

  outputmatrix[position1[1],position1[2]] = 0
  outputmatrix[position2[1],position2[2]] = 0
  outputmatrix[position3[1],position3[2]] = 0

  println("At the end, it looks like this:")
  println(outputmatrix)

  outputmatrix
end

function updatepositions(matrix,oldpositions,newposition)

  # from largest to smallest
  positions = orderpositions(matrix,oldpositions)

  position1 = positions[1]
  position2 = positions[2]
  position3 = positions[3]

  # value1 >= value2 >= value3
  value1 = matrix[position1[1],:][position1[2]]
  value2 = matrix[position2[1],:][position2[2]]
  value3 = matrix[position3[1],:][position3[2]]

  newvalue = matrix[newposition[1],:][newposition[2]]

  if newvalue > value1
    positions[3] = positions[2]
    positions[2] = positions[1]
    positions[1] = newposition
  elseif newvalue > value2
    positions[3]= positions[2]
    positions[2] = newposition
    # no change in positions[1]
  elseif newvalue > value3
    positions[3] = newposition
    # no change in positions[1] or positions[2]
  else
    # nothing to do
  end

  positions

end

function orderpositions(matrix,positions)

  @assert size(positions)[1] == 3

  position1 = positions[1]
  position2 = positions[2]
  position3 = positions[3]

  value1 = matrix[position1[1],:][position1[2]]
  value2 = matrix[position2[1],:][position2[2]]
  value3 = matrix[position3[1],:][position3[2]]

  if(value1 >= value2 >= value3)
    orderedpositions = [position1,position2,position3]
  elseif(value1 >= value3 >= value2)
    orderedpositions = [position1,position3,position2]
  elseif(value2 >= value1 >= value3)
    orderedpositions = [position2,position1,position3]
  elseif(value2 >= value3 >= value1)
    orderedpositions =  [position2,position3,position1]
  elseif(value3 >= value1 >= value2)
    orderedpositions = [position3,position1,position2]
  elseif(value3 >= value2 >= value1)
    orderedpositions = [position3,position2,position1]
  end

  @assert size(orderedpositions)[1] == 3

  orderedpositions
end

M = [5 10 -5 22; 1 33 15 3; 8 29 12 1;3 11 39 20]

main(M)


text:
############
#Observação#
############

# Para usar essa função é necessario mandar uma matriz AA que seja do tipo float
# tentei fazer a conversão de tipo porem os resultados saiam como se fossem um número inteiro.

function eliminacao00(A)

	(m, n) = size(A);

	for i = 1:n-1
		pivot = A[i,i];
		for j = i+1:m

			mult = A[j,i]/pivot;

			A[j,:] = A[j,:] - mult*A[i,:];

		end

	end

	return A;

end


text:
using JuMP
using LightGraphs, SimpleWeightedGraphs
using Random

lazy_cec_constraints = 0

function cycle_elimination_constraints!(
    model,
    graph::SimpleWeightedGraph,
    k::Int)

    global lazy_cec_constraints = 0
    @constraint(model, sum(variable_by_name(model, "y[1,$j]") for j in 2:nv(graph)) == 1) # Constraint 32
    @constraint(model, sum(variable_by_name(model, "y[$j,1]") for j in 2:nv(graph)) >= 1) # Constraint 33
    @constraint(model, sum(variable_by_name(model, "y[$j,1]") for j in 2:nv(graph)) <= k-1) # Constraint 34
    for i in 2:nv(graph)
        @constraint(model, sum(variable_by_name(model, "y[$i,$j]")
                               for j in 1:nv(graph) if has_edge(graph, i, j))
                    >= sum(variable_by_name(model, "y[$j,$i]")
                           for j in 1:nv(graph) if has_edge(graph, j, i))) # Constraint 35

        @constraint(model, sum(variable_by_name(model, "y[$i,$j]")
                               for j in 1:nv(graph) if has_edge(graph, i, j))
                    <= sum(variable_by_name(model, "y[$j,$i]")
                           for j in 1:nv(graph) if has_edge(graph, j, i)) * (k-1)) # Constraint 36
        @constraint(model, sum(variable_by_name(model, "y[$j,$i]") for j in 1:nv(graph) if has_edge(graph, j, i)) <= 1) # Constraint 37
    end

    MOI.set(model,
            MOI.LazyConstraintCallback(),
            cec_lazy_clause_generator(model, graph, k))
end

function cec_lazy_clause_generator(model, graph::SimpleWeightedGraph, k::Int)
    es = collect(SimpleWeightedEdge, edges(graph))
    sort!(es, by=dst)
    sort!(es, by=src, alg=MergeSort)
    y::Dict{Tuple{Int,Int},VariableRef} = Dict{Tuple{Int,Int},VariableRef}()
    for e in es
        i::Int = src(e)
        j::Int = dst(e)
        y[i,j] = variable_by_name(model, "y[$i,$j]")
        y[j,i] = variable_by_name(model, "y[$j,$i]")
    end
    return function (cb_data)
        status = callback_node_status(cb_data, model)
        if status == MOI.CALLBACK_NODE_STATUS_FRACTIONAL
            @debug "Found fractional solution"
            return
        elseif status == MOI.CALLBACK_NODE_STATUS_INTEGER
            @debug "Found integer solution"
            # x_val::Dict{Int, } = Dict{Int, Float64}()
            y_val::Dict{Tuple{Int,Int},Int} = Dict{Tuple{Int,Int},Int}()
            for e in es
                i::Int = src(e)
                j::Int = dst(e)
                y_val[i,j] = round(callback_value(cb_data, y[i,j]))
                y_val[j,i] = round(callback_value(cb_data, y[j,i]))
            end
            cycle = find_cycle(graph, y_val,k)
            if !isnothing(cycle)
                cyclelength::Int = length(cycle)
                @debug "Detected cycle with length $cyclelength: $cycle."
                constraint = @build_constraint(sum(y[i,j]
                                                    for i in cycle
                                                        for j in cycle
                                                            if has_edge(graph, i, j))
                                               <= cyclelength - 1) # Constraint 38

                # for debugging:
                # arcs::Vector{String} = Vector{String}()
                # for i in cycle
                #    for j in cycle
                #        if has_edge(graph, i, j)
                #            push!(arcs, "y[$i,$j]")
                #        end
                #    end
                # end
                # @debug join(arcs, " + ") * " <= $(cyclelength - 1)"

                global lazy_cec_constraints += 1
                MOI.submit(model, MOI.LazyConstraint(cb_data), constraint)
            end
        else
            @assert status == MOI.CALLBACK_NODE_STATUS_UNKNOWN
            @debug "Found either fractional or integer solution"
        end
    end
end

function find_cycle(
    graph::SimpleWeightedGraph,
    y_val::Dict{Tuple{Int,Int},Int64},
    k::Int)::Union{Nothing,Set{Int}}
    visited::Set{Int} = Set{Int}()
    queue::Set{Int} = Set{Int}([1])
    while !isempty(queue)
        @assert !isempty(queue) "queue should not be empty."
        c::Int = pop!(queue)
        ns = neighbors(graph, c)
        for n::Int in ns
            if y_val[c, n] == 1 && n != 1
                push!(queue, n)
            end
        end
        push!(visited, c)
    end
    if length(visited) == k + 1
        return # no cycles
    end
    @assert length(visited) <= k "Solution should be a forest, but only has one tree."

    to_check::Set{Int} = Set{Int}(2:nv(graph))

    filter!(v -> !(v in visited), to_check)
    @assert length(to_check) == (nv(graph) - 1) - (length(visited) - 1) "Should have pruned $(length(visited) - 1) of $(nv(graph) - 1) items. However only $(length(to_check)) left."

    while !isempty(to_check)
        @assert isempty(queue) "queue should be empty."
        push!(queue, pop!(to_check))
        visited = Set{Int}()
        @assert isempty(visited) "visited should be empty."
        while !isempty(queue)
            @assert !isempty(queue) "queue should not be empty."
            c::Int = pop!(queue)
            # @debug "Visiting node $c"
            push!(visited, c)
            ns = neighbors(graph, c)
            for n::Int in ns
                if y_val[c, n] == 1 && n != 1
                    if n in visited
                        return visited
                    end
                    push!(queue, n)
                end
            end
        end
        filter!(v -> !(v in visited), to_check)
    end

    @warn "No cycle detected. Something is wrong with the implementation."

end

max_stars_repo_name: Entze/jkMST

max_issues_repo_name: Entze/jkMST

max_forks_repo_name: Entze/jkMST


text:
export MixedMatHRep, MixedMatVRep

# H-Representation
"""
    hrep(A::AbstractMatrix, b::AbstractVector, linset::BitSet=BitSet())

Creates an H-representation for the polyhedron defined by the inequalities ``\\langle A_i, x \\rangle = b_i`` if `i in linset`
and ``\\langle A_i, x \\rangle \\le b_i`` otherwise where ``A_i`` is the ``i``th row of `A`, i.e. `A[i,:]` and ``b_i`` is `b[i]`.
"""
hrep(A::AbstractMatrix, b::AbstractVector, linset::BitSet=BitSet()) = MixedMatHRep(A, b, linset)

# No copy since I do not modify anything and a copy is done when building a polyhedron
mutable struct MixedMatHRep{T, MT<:AbstractMatrix{T}} <: MixedHRep{T}
    # Ax <= b
    A::MT
    b::Vector{T}
    linset::BitSet

    function MixedMatHRep{T, MT}(A::MT, b::AbstractVector, linset::BitSet) where {T, MT}
        if size(A, 1) != length(b)
            error("The length of b must be equal to the number of rows of A")
        end
        if !isempty(linset) && last(linset) > length(b)
            error("The elements of linset should be between 1 and the number of rows of A/length of b")
        end
        new{T, typeof(A)}(A, b, linset)
    end
end
FullDim(rep::MixedMatHRep) = size(rep.A, 2)

similar_type(::Type{<:MixedMatHRep{S, MT}}, ::FullDim, ::Type{T}) where {S, T, MT} = MixedMatHRep{T, similar_type(MT, T)}
hvectortype(p::Type{MixedMatHRep{T, MT}}) where {T, MT} = vectortype(MT)
fulltype(::Type{MixedMatHRep{T, MT}}) where {T, MT} = MixedMatHRep{T, MT}

MixedMatHRep{T}(A::AbstractMatrix{T}, b::AbstractVector{T}, linset::BitSet) where {T} = MixedMatHRep{T, typeof(A)}(A, b, linset)
MixedMatHRep(A::AbstractMatrix{T}, b::AbstractVector{T}, linset::BitSet) where T = MixedMatHRep{T, typeof(A)}(A, b, linset)
MixedMatHRep{U}(A::AbstractMatrix{S}, b::AbstractVector{T}, linset::BitSet) where {S, T, U} = MixedMatHRep{U}(AbstractMatrix{U}(A), AbstractVector{U}(b), linset)
function MixedMatHRep(A::AbstractMatrix{S}, b::AbstractVector{T}, linset::BitSet) where {S, T}
    U = promote_type(S, T)
    MixedMatHRep(AbstractMatrix{U}(A), AbstractVector{U}(b), linset)
end
MixedMatHRep(A::AbstractMatrix{T}, b::AbstractVector{T}, linset::BitSet) where T <: Real = MixedMatHRep{T}(A, b, linset)

MixedMatHRep(h::HRep{T}) where {T} = MixedMatHRep{T}(h)
MixedMatHRep{T}(h::HRep) where {T} = convert(MixedMatHRep{T, hmatrixtype(typeof(h), T)}, h)

function MixedMatHRep{T, MT}(d::FullDim, hyperplanes::HyperPlaneIt{T},
                             halfspaces::HalfSpaceIt{T}) where {T, MT}
    nhyperplane = length(hyperplanes)
    nhrep = nhyperplane + length(halfspaces)
    A = emptymatrix(MT, nhrep, fulldim(d))
    b = Vector{T}(undef, nhrep)
    linset = BitSet(1:nhyperplane)
    for (i, h) in enumerate(hyperplanes)
        A[i,:] = h.a
        b[i] = h.β
    end
    for (i, h) in enumerate(halfspaces)
        A[nhyperplane+i,:] = h.a
        b[nhyperplane+i] = h.β
    end
    MixedMatHRep{T, MT}(A, b, linset)
end

Base.copy(ine::MixedMatHRep{T, MT}) where {T, MT} = MixedMatHRep{T, MT}(copy(ine.A), copy(ine.b), copy(ine.linset))

Base.isvalid(hrep::MixedMatHRep{T}, idx::HIndex{T}) where {T} = 0 < idx.value <= size(hrep.A, 1) && (idx.value in hrep.linset) == islin(idx)
done(idxs::HIndices{T, <:MixedMatHRep{T}}, idx::HIndex{T}) where {T} = idx.value > size(idxs.rep.A, 1)
Base.get(hrep::MixedMatHRep{T}, idx::HIndex{T}) where {T} = valuetype(idx)(hrep.A[idx.value,:], hrep.b[idx.value])

# V-Representation
"""
    vrep(V::AbstractMatrix, R::AbstractMatrix, Rlinset::BitSet=BitSet())

Creates a V-representation for the polyhedron defined by the points ``V_i``,
lines ``R_i`` if `i in Rlinset` and rays ``R_i`` otherwise
where ``V_i`` (resp. ``R_i``) is the ``i``th row of `V` (resp. `R`), i.e. `V[i,:]` (resp. `R[i,:]`).
"""
vrep(V::AbstractMatrix, R::AbstractMatrix, Rlinset::BitSet=BitSet()) = MixedMatVRep(V, R, Rlinset)
vrep(V::AbstractMatrix) = vrep(V, similar(V, 0, size(V, 2)))

mutable struct MixedMatVRep{T, MT<:AbstractMatrix{T}} <: MixedVRep{T}
    V::MT # each row is a vertex
    R::MT # each row is a ray
    Rlinset::BitSet

    function MixedMatVRep{T, MT}(V::MT, R::MT, Rlinset::BitSet) where {T, MT}
        if iszero(size(V, 1)) && !iszero(size(R, 1))
            vconsistencyerror()
        end
        if (length(R) > 0 && length(V) > 0 && size(R, 2) != size(V, 2))
            error("Dimension of rays does not match and dimension of points")
        end
        if !isempty(Rlinset) && last(Rlinset) > size(R, 1)
            error("The elements of Rlinset should be between 1 and the number of rows of R")
        end
        new{T, MT}(V, R, Rlinset)
    end
end
FullDim(rep::MixedMatVRep) = size(rep.V, 2)

similar_type(::Type{<:MixedMatVRep{S, MT}}, ::FullDim, ::Type{T}) where {S, T, MT} = MixedMatVRep{T, similar_type(MT, T)}
vvectortype(::Type{MixedMatVRep{T, MT}}) where {T, MT} = vectortype(MT)
fulltype(::Type{MixedMatVRep{T, MT}}) where {T, MT} = MixedMatVRep{T, MT}

MixedMatVRep{T}(V::MT, R::MT, Rlinset::BitSet) where {T, MT<:AbstractMatrix{T}} = MixedMatVRep{T, MT}(V, R, Rlinset)
MixedMatVRep(V::MT, R::MT, Rlinset::BitSet) where {T, MT<:AbstractMatrix{T}} = MixedMatVRep{T, MT}(V, R, Rlinset)
MixedMatVRep{U}(V::AbstractMatrix{S}, R::AbstractMatrix{T}, Rlinset::BitSet) where {S, T, U} = MixedMatVRep{U}(AbstractMatrix{U}(V), AbstractMatrix{U}(R), Rlinset)
function MixedMatVRep(V::AbstractMatrix{S}, R::AbstractMatrix{T}, Rlinset::BitSet) where {S, T}
    U = promote_type(S, T)
    MixedMatVRep(AbstractMatrix{U}(V), AbstractMatrix{U}(R), Rlinset)
end

MixedMatVRep(v::VRep{T}) where {T} = MixedMatVRep{T}(v)
MixedMatVRep{T}(v::VRep) where {T} = convert(MixedMatVRep{T, vmatrixtype(typeof(v), T)}, v)

function MixedMatVRep{T, MT}(d::FullDim, vits::VIt{T}...) where {T, MT}
    points, lines, rays = fillvits(d, vits...)
    npoint = length(points)
    nline = length(lines)
    nray = length(rays)
    N = fulldim(d)
    V = emptymatrix(MT, npoint, N)
    R = emptymatrix(MT, nline + nray, N)
    Rlinset = BitSet()
    function _fill!(M, linset, offset, ps)
        for (i, p) in enumerate(ps)
            M[offset + i, :] = coord(p)
            if islin(p)
                push!(linset, offset + i)
            end
        end
    end
    _fill!(V, nothing, 0, points)
    _fill!(R, Rlinset, 0, lines)
    _fill!(R, Rlinset, nline, rays)
    MixedMatVRep{T, MT}(V, R, Rlinset)
end

Base.copy(ext::MixedMatVRep{T, MT}) where {T, MT} = MixedMatVRep{T, MT}(copy(ext.V), copy(ext.R), copy(ext.Rlinset))

_mat(hrep::MixedMatVRep, ::PIndex) = hrep.V
_mat(hrep::MixedMatVRep, ::RIndex) = hrep.R
_islin(hrep::MixedMatVRep, ::PIndex) = false
_islin(hrep::MixedMatVRep, idx::RIndex) = idx.value in hrep.Rlinset

Base.isvalid(vrep::MixedMatVRep{T}, idx::VIndex{T}) where {T} = 0 < idx.value <= size(_mat(vrep, idx), 1) && _islin(vrep, idx) == islin(idx)
done(idxs::VIndices{T, <:MixedMatVRep{T}}, idx::VIndex{T}) where {T} = idx.value > size(_mat(idxs.rep, idx), 1)
Base.get(vrep::MixedMatVRep{T}, idx::VIndex{T}) where {T} = valuetype(idx)(_mat(vrep, idx)[idx.value, :])

# sparse halfspaces does not mean sparse points
dualtype(::Type{MixedMatHRep{T, MT}}) where {T, MT} = MixedMatVRep{T, Matrix{T}}
dualtype(::Type{MixedMatVRep{T, MT}}) where {T, MT} = MixedMatHRep{T, Matrix{T}}
function dualfullspace(h::MixedMatHRep, d::FullDim, ::Type{T}) where {T}
    N = fulldim(d)
    MixedMatVRep{T}(zeros(T, 1, N), Matrix(one(T) * I, N, N), BitSet(1:N))
end


text:
# Tests for README & Documentation Examples

# Manual control
money = 1USD
magn = money.val
symb = currency(money)
a = π^2 * magn

# Give change
COINS = [500EUR, 200EUR, 100EUR, 50EUR, 20EUR, 10EUR, 5EUR, 2EUR, 1EUR, 0.5EUR,
    0.2EUR, 0.1EUR, 0.05EUR, 0.02EUR, 0.01EUR]
function change(amount::Monetary{:EUR,Int})
    coins = Dict{Monetary{:EUR,Int}, Int}()
    for denomination in COINS
        coins[denomination], amount = divrem(amount, denomination)
    end
    coins
end

@testset "Doc Examples" begin
    @test Monetary(symb, round(Int, 100a)) == 9.87USD
    @test sum([k*v for (k, v) in change(167.25EUR)]) == 167.25EUR
end