text:
module Integers

##############################################################################
##
## TODO
##
##############################################################################

# - Should remove module and rely on package-level module as per convention
# - Methods e.g., promote, convert, macro +, -, *, /

##############################################################################
##
## Exports
##
##############################################################################

export IntegerMod, makeModular, extendedEuclideanAlgorithm,
       gcd

##############################################################################
##
## Implementation
##
##############################################################################

import Base.show,
       Base.length,
       Base.convert,
       Base.divrem,
       Base.promote_rule,
       Base.promote,
       Base.abs,
       Base.div,
       Base.rem,
       Base.zero,
       Base.zeros

type IntegerMod{P} <: Number
    n::Integer
    p::Integer

    function IntegerMod()
      new(0, P)
    end
    function IntegerMod(N::Integer)
        new(mod(N,P), P)
    end
end

function makeModular(n::Integer, p::Integer)
    return IntegerMod{p}(n)
end

convert{P}(::Type{IntegerMod{P}}, n::Integer) = IntegerMod{P}(n)
convert{P}(::Type{IntegerMod{P}}, a::Array) = [IntegerMod{P}(n) for n in a]

promote_rule{P}(::Type{IntegerMod{P}}, ::Type{Integer}) = IntegerMod{P}

function show{P}(io::IO, a::IntegerMod{P})
    print(io, string(a.n, " (mod ", P, ")"))
end

function =={P}(a::IntegerMod{P}, b::IntegerMod{P})
  return a.n == b.n
end

function =={P}(a::IntegerMod{P}, b::Integer)
  return a.n == mod(b, P)
end

function =={P}(a::Integer, b::IntegerMod{P})
  return mod(a, P) == b.n
end

function +{P}(a::IntegerMod{P}, b::IntegerMod{P})
  return IntegerMod{P}(a.n + b.n)
end

function +{P}(a::IntegerMod{P}, b::Integer)
  return IntegerMod{P}(a.n + b)
end

function +{P}(a::Integer, b::IntegerMod{P})
  return IntegerMod{P}(a + b.n)
end

function -{P}(a::IntegerMod{P})
  return IntegerMod{P}(-a.n)
end

function -{P}(a::IntegerMod{P}, b::IntegerMod{P})
  return IntegerMod{P}(a.n - b.n)
end

function -{P}(a::IntegerMod{P}, b::Integer)
  return IntegerMod{P}(a.n - b)
end

function -{P}(a::IntegerMod{P}, b::IntegerMod{P})
  return return IntegerMod{P}(a - b.n)
end

function *{P}(a::IntegerMod{P}, b::IntegerMod{P})
  return IntegerMod{P}(a.n * b.n)
end

function *{P}(a::IntegerMod{P}, b::Bool)
  return
end

function *{P}(a::Bool, b::IntegerMod{P})
  return
end

function *{P}(a::IntegerMod{P}, b::Integer)
  return IntegerMod{P}(a.n * b)
end

function *{P}(a::Integer, b::IntegerMod{P})
  return IntegerMod{P}(a * b.n)
end

function zero{P}(a::IntegerMod{P})
  return IntegerMod{P}(0)
end

function zero{P}(a::Type{IntegerMod{P}})
  return IntegerMod{P}(0)
end

function zeros{P}(T::Type{IntegerMod{P}}, n = 0)
  result = Array(IntegerMod{P}, n)
  for i in 1:n
    result[i] = 0
  end
  return result
end

function divrem{P}(a::IntegerMod{P}, b::IntegerMod{P})
    q, r = divrem(a.n, b.n)
  return (IntegerMod{P}(q), IntegerMod{P}(r))
end

function /{P}(a::IntegerMod{P}, b::IntegerMod{P})
    return a * inverse(b)
end

function abs{P}(a::IntegerMod{P})
  return abs(a.n)
end

function inverse{P}(a::IntegerMod{P})
    x, y, d = extendedEuclideanAlgorithm(a.n, P)
    return IntegerMod{P}(x)
end

function gcd(a, b)
  if abs(b) > abs(a)
    return gcd(b, a)
  end

  while abs(b) > 0
    _, r = divrem(a,b)
    a, b = b, r
  end

  return a
end

function extendedEuclideanAlgorithm{T<:Number}(a::T, b::T)
    if abs(b) > abs(a)
        x, y, d = extendedEuclideanAlgorithm(b,a)
        return (y,x,d)
    end

    if abs(b) == 0
        return (1, 0, a)
    end

    x1, x2, y1, y2 = convert(T, 0), convert(T, 1), convert(T, 1), convert(T, 0)
    while abs(b) > 0
        q, r = divrem(a, b)
        x = x2 - q*x1
        y = y2 - q*y1
        a, b, x2, x1, y2, y1 = b, r, x1, x, y1, y
    end

    return (x2, y2, a)
end

end # module Integer


text:
################################################################################
#                           SUPPORT ITERATION METHODS
################################################################################
## Form a placeholder parameter reference given the object index
# IndependentParameterIndex
function _temp_parameter_ref(model::InfiniteOpt.InfiniteModel,
    index::InfiniteOpt.IndependentParameterIndex
    )::InfiniteOpt.IndependentParameterRef
    return InfiniteOpt.IndependentParameterRef(model, index)
end

# DependentParametersIndex
function _temp_parameter_ref(model::InfiniteOpt.InfiniteModel,
    index::InfiniteOpt.DependentParametersIndex
    )::InfiniteOpt.DependentParameterRef
    idx = InfiniteOpt.DependentParameterIndex(index, 1)
    return InfiniteOpt.DependentParameterRef(model, idx)
end

# Return the collected supports of an infinite parameter
function _collected_supports(
    pref::Union{InfiniteOpt.IndependentParameterRef, InfiniteOpt.DependentParameterRef}
    )::Vector
    supp_dict = InfiniteOpt._parameter_supports(pref)
    supp_list = collect(keys(supp_dict))
    # append a placeholder NaN support at the end to be used for efficient combinatorics
    return push!(supp_list, map(i -> NaN, first(supp_list)))
end

# Return the collected support labels of an infinite parameter
function _collected_support_labels(
    pref::Union{InfiniteOpt.IndependentParameterRef, InfiniteOpt.DependentParameterRef},
    supports::Vector
    )::Vector{Set{DataType}}
    supp_dict = InfiniteOpt._parameter_supports(pref)
    default = Set{DataType}()
    return map(k -> get(supp_dict, k, default), supports)
end

"""
    set_parameter_supports(trans_model::JuMP.Model,
                           inf_model::InfiniteOpt.InfiniteModel)::Nothing

Collect the infinite parameter supports stored in their respective dictionaries
form `inf_model` and process them into a tuple of vectors where each vector
contains the collected supports of a particular infinite parameter. These support
collections are ordered in accordance with the definition order of the
parameters (i.e., their object numbers). A support collection assocciated with
an independent will be a `Vector{Float64}` and a support collection associated
with a group of dependent parameters will be a `Vector{Vector{Float64}}`. Note
that each collection vector will include an extra final placeholder element
comprised of `NaN`s for convenience in generating support indices via
[`support_index_iterator`](@ref). This also gathers the associated support labels.

Before this is all done, `InfiniteOpt.add_generative_supports` is invoked as needed.
"""
function set_parameter_supports(trans_model::JuMP.Model,
    inf_model::InfiniteOpt.InfiniteModel
    )::Nothing
    # gather the basic information
    param_indices = InfiniteOpt._param_object_indices(inf_model)
    prefs = map(idx -> _temp_parameter_ref(inf_model, idx), param_indices)
    data = transcription_data(trans_model)
    # check and add supports to prefs as needed
    for pref in prefs
        InfiniteOpt.add_generative_supports(pref)
        if InfiniteOpt.has_internal_supports(pref)
            data.has_internal_supports = true
        end
    end
    # build and add the support/label tuples
    supps = Tuple(_collected_supports(pref) for pref in prefs)
    labels = Tuple(_collected_support_labels(pref, supps[i])
                   for (i, pref) in enumerate(prefs))
    data.supports = supps
    data.support_labels = labels
    return
end

################################################################################
#                        VARIABLE INITIALIZATION METHODS
################################################################################
"""
    transcribe_hold_variables!(trans_model::JuMP.Model,
                               inf_model::InfiniteOpt.InfiniteModel)::Nothing

Create a transcription variable (i.e., a JuMP variable) for each `HoldVariable`
stored in `inf_model` and add it to `trans_model`. The variable mapping is
also stored in `TranscriptionData.finvar_mappings` which enables
[`transcription_variable`](@ref) and [`lookup_by_support`](@ref).
"""
function transcribe_hold_variables!(trans_model::JuMP.Model,
    inf_model::InfiniteOpt.InfiniteModel
    )::Nothing
    for (index, object) in InfiniteOpt._data_dictionary(inf_model, InfiniteOpt.HoldVariable)
        hvref = InfiniteOpt._make_variable_ref(inf_model, index)
        vref = JuMP.add_variable(trans_model,
                                 JuMP.ScalarVariable(object.variable.info),
                                 object.name)
        transcription_data(trans_model).finvar_mappings[hvref] = vref
    end
    return
end

# Return a proper scalar variable info object given on with a start value function
function _format_infinite_info(var::InfiniteOpt.InfiniteVariable,
    support::Vector{Float64}
    )::JuMP.VariableInfo
    # generate the start value
    info = var.info
    if var.is_vector_start
        start = info.start(support)
    else
        start = info.start(Tuple(support, var.parameter_refs)...)
    end
    # make the info and return
    return JuMP.VariableInfo(info.has_lb, info.lower_bound, info.has_ub,
                             info.upper_bound, info.has_fix, info.fixed_value,
                             info.has_start, start, info.binary, info.integer)
end

"""
    transcribe_infinite_variables!(trans_model::JuMP.Model,
                                   inf_model::InfiniteOpt.InfiniteModel)::Nothing

Create transcription variables (i.e., JuMP variables) for each `InfiniteVariable`
stored in `inf_model` and add them to `trans_model`. The variable mappings are
also stored in `TranscriptionData.infvar_mappings` in accordance with
`TranscriptionData.infvar_lookup` which enable [`transcription_variable`](@ref)
and [`lookup_by_support`](@ref). Note that the supports will not be generated
until `InfiniteOpt.variable_supports` is invoked via `InfiniteOpt.supports`.
Note that `TranscriptionData.infvar_support_labels` is also populated.
"""
function transcribe_infinite_variables!(trans_model::JuMP.Model,
    inf_model::InfiniteOpt.InfiniteModel
    )::Nothing
    for (index, object) in InfiniteOpt._data_dictionary(inf_model, InfiniteOpt.InfiniteVariable)
        # get the basic variable information
        var = object.variable
        base_name = object.name
        param_nums = var.parameter_nums
        # prepare for iterating over its supports
        supp_indices = support_index_iterator(trans_model, var.object_nums)
        vrefs = Vector{JuMP.VariableRef}(undef, length(supp_indices))
        labels = Vector{Set{DataType}}(undef, length(supp_indices))
        lookup_dict = Dict{Vector{Float64}, Int}()
        # create a variable for each support
        for (counter, i) in enumerate(supp_indices)
            raw_supp = index_to_support(trans_model, i)
            supp = raw_supp[param_nums]
            info = _format_infinite_info(var, supp)
            var_name = string(base_name, "(support: ", counter, ")")
            @inbounds vrefs[counter] = JuMP.add_variable(trans_model,
                                                         JuMP.ScalarVariable(info),
                                                         var_name)
            lookup_dict[supp] = counter
            @inbounds labels[counter] = index_to_labels(trans_model, i)
        end
        # save the transcription information
        ivref = InfiniteOpt._make_variable_ref(inf_model, index)
        data = transcription_data(trans_model)
        data.infvar_lookup[ivref] = lookup_dict
        data.infvar_mappings[ivref] = vrefs
        data.infvar_support_labels[ivref] = labels
    end
    return
end

# Return a proper scalar variable info object given on with a start value function
function _format_derivative_info(d::InfiniteOpt.Derivative,
    support::Vector{Float64}
    )::JuMP.VariableInfo
    # generate the start value
    info = d.info
    if d.is_vector_start
        start = info.start(support)
    else
        prefs = InfiniteOpt.raw_parameter_refs(d.variable_ref)
        start = info.start(Tuple(support, prefs)...)
    end
    # make the info and return
    return JuMP.VariableInfo(info.has_lb, info.lower_bound, info.has_ub,
                             info.upper_bound, info.has_fix, info.fixed_value,
                             info.has_start, start, info.binary, info.integer)
end

"""
    transcribe_derivative_variables!(trans_model::JuMP.Model,
                                     inf_model::InfiniteOpt.InfiniteModel)::Nothing

Create transcription variables (i.e., JuMP variables) for each `Derivative`
stored in `inf_model` and add them to `trans_model`. The variable mappings are
also stored in `TranscriptionData.infvar_mappings` in accordance with
`TranscriptionData.infvar_lookup` which enable [`transcription_variable`](@ref)
and [`lookup_by_support`](@ref). Note that the supports will not be generated
until `InfiniteOpt.variable_supports` is invoked via `InfiniteOpt.supports`. The
futher derivative evaluation constraints are added when
`transcribe_derivative_evaluations!` is invoked. Note that
`TranscriptionData.infvar_support_labels` is also populated.
"""
function transcribe_derivative_variables!(trans_model::JuMP.Model,
    inf_model::InfiniteOpt.InfiniteModel
    )::Nothing
    for (index, object) in InfiniteOpt._data_dictionary(inf_model, InfiniteOpt.Derivative)
        # get the basic variable information
        dref = InfiniteOpt._make_variable_ref(inf_model, index)
        d = object.variable
        base_name = InfiniteOpt.variable_string(JuMP.REPLMode, dispatch_variable_ref(dref))
        param_nums = InfiniteOpt._parameter_numbers(d.variable_ref)
        obj_nums = InfiniteOpt._object_numbers(d.variable_ref)
        # prepare for iterating over its supports
        supp_indices = support_index_iterator(trans_model, obj_nums)
        vrefs = Vector{JuMP.VariableRef}(undef, length(supp_indices))
        labels = Vector{Set{DataType}}(undef, length(supp_indices))
        lookup_dict = Dict{Vector{Float64}, Int}()
        # create a variable for each support
        for (counter, i) in enumerate(supp_indices)
            raw_supp = index_to_support(trans_model, i)
            supp = raw_supp[param_nums]
            info = _format_derivative_info(d, supp)
            deriv_name = string(base_name, "(support: ", counter, ")")
            @inbounds vrefs[counter] = JuMP.add_variable(trans_model,
                                                         JuMP.ScalarVariable(info),
                                                         deriv_name)
            lookup_dict[supp] = counter
            @inbounds labels[counter] = index_to_labels(trans_model, i)
        end
        # save the transcription information
        data = transcription_data(trans_model)
        data.infvar_lookup[dref] = lookup_dict
        data.infvar_mappings[dref] = vrefs
        data.infvar_support_labels[dref] = labels
    end
    return
end

# Setup the mapping for a given reduced variable
function _set_reduced_variable_mapping(trans_model::JuMP.Model,
    var::InfiniteOpt.ReducedVariable,
    rvref::InfiniteOpt.GeneralVariableRef,
    index_type
    )::Nothing
    param_nums = var.parameter_nums
    ivref = var.infinite_variable_ref
    ivref_param_nums = InfiniteOpt._parameter_numbers(ivref)
    eval_supps = var.eval_supports
    # prepare for iterating over its supports
    supp_indices = support_index_iterator(trans_model, var.object_nums)
    vrefs = Vector{JuMP.VariableRef}(undef, length(supp_indices))
    labels = Vector{Set{DataType}}(undef, length(supp_indices))
    lookup_dict = Dict{Vector{Float64}, Int}()
    counter = 1
    # map a variable for each support
    for i in supp_indices
        raw_supp = index_to_support(trans_model, i)
        # ensure this support is valid with the reduced restriction
        if any(!isnan(raw_supp[ivref_param_nums[k]]) && raw_supp[ivref_param_nums[k]] != v
               for (k, v) in eval_supps)
            continue
        end
        # map to the current transcription variable
        supp = raw_supp[param_nums]
        ivref_supp = [haskey(eval_supps, j) ? eval_supps[j] : raw_supp[k]
                      for (j, k) in enumerate(ivref_param_nums)]
        @inbounds vrefs[counter] = lookup_by_support(trans_model, ivref, ivref_supp)
        lookup_dict[supp] = counter
        @inbounds labels[counter] = index_to_labels(trans_model, i)
        counter += 1
    end
    # truncate vrefs if any supports were skipped because of dependent parameter supps
    deleteat!(vrefs, counter:length(vrefs))
    deleteat!(labels, counter:length(vrefs))
    # save the transcription information
    data = transcription_data(trans_model)
    data.infvar_lookup[rvref] = lookup_dict
    data.infvar_mappings[rvref] = vrefs
    data.infvar_support_labels[rvref] = labels
    return
end

# Empty mapping dispatch for infinite parameter functions
function _set_reduced_variable_mapping(trans_model::JuMP.Model,
    var::InfiniteOpt.ReducedVariable,
    rvref::InfiniteOpt.GeneralVariableRef,
    index_type::Type{InfiniteOpt.ParameterFunctionIndex}
    )::Nothing
    return
end

"""
    transcribe_reduced_variables!(trans_model::JuMP.Model,
                                  inf_model::InfiniteOpt.InfiniteModel)::Nothing

Map each `ReducedVariable` in `inf_model` to transcription variables stored in
`trans_model`. The variable mappings are also stored in
`TranscriptionData.infvar_mappings` in accordance with
`TranscriptionData.infvar_lookup` which enable [`transcription_variable`](@ref)
and [`lookup_by_support`](@ref). Note that [`transcribe_infinite_variables!`](@ref)
must be called first. Note that the supports will not be generated
until `InfiniteOpt.variable_supports` is invoked via `InfiniteOpt.supports`.
Note that `TranscriptionData.infvar_support_labels` is also populated.
"""
function transcribe_reduced_variables!(trans_model::JuMP.Model,
    inf_model::InfiniteOpt.InfiniteModel
    )::Nothing
    for (index, object) in InfiniteOpt._data_dictionary(inf_model, InfiniteOpt.ReducedVariable)
        # get the basic variable information
        var = object.variable
        rvref = InfiniteOpt._make_variable_ref(inf_model, index)
        # setup the mappings
        ivref = InfiniteOpt.infinite_variable_ref(rvref)
        _set_reduced_variable_mapping(trans_model, var, rvref,
                                      InfiniteOpt._index_type(ivref))
    end
    return
end

# Override the info of the jump variable with the point variable's if it is different
function _update_point_info(gvref::InfiniteOpt.GeneralVariableRef,
    vref::JuMP.VariableRef
    )::Nothing
    pvref = InfiniteOpt.dispatch_variable_ref(gvref)
    if JuMP.has_lower_bound(pvref)
        if JuMP.is_fixed(vref)
            JuMP.unfix(vref)
        end
        JuMP.set_lower_bound(vref, JuMP.lower_bound(pvref))
    end
    if JuMP.has_upper_bound(pvref)
        if JuMP.is_fixed(vref)
            JuMP.unfix(vref)
        end
        JuMP.set_upper_bound(vref, JuMP.upper_bound(pvref))
    end
    if JuMP.is_fixed(pvref)
        JuMP.fix(vref, JuMP.fix_value(pvref), force = true)
    end
    if JuMP.is_binary(pvref) && !JuMP.is_binary(vref)
        if JuMP.is_integer(vref)
            JuMP.unset_integer(vref)
        end
        JuMP.set_binary(vref)
    elseif JuMP.is_integer(pvref) && !JuMP.is_integer(vref)
        if JuMP.is_binary(vref)
            JuMP.unset_binary(vref)
        end
        JuMP.set_integer(vref)
    end
    if JuMP.start_value(pvref) != JuMP.start_value(vref)
        JuMP.set_start_value(vref, JuMP.start_value(pvref))
    end
    return
end

"""
    transcribe_point_variables!(trans_model::JuMP.Model,
                               inf_model::InfiniteOpt.InfiniteModel)::Nothing

Map each `PointVariable` in `inf_model` to a transcription variable stored in
`trans_model`. The variable mapping is also stored in
`TranscriptionData.finvar_mappings` which enables
[`transcription_variable`](@ref) and [`lookup_by_support`](@ref). Note that
[`transcribe_infinite_variables!`](@ref) must be called first and that the
info constraints associated with the transcription variable will be updated in
accordance with the point variable.
"""
function transcribe_point_variables!(trans_model::JuMP.Model,
    inf_model::InfiniteOpt.InfiniteModel
    )::Nothing
    for (index, object) in InfiniteOpt._data_dictionary(inf_model, InfiniteOpt.PointVariable)
        # get the basic variable information
        var = object.variable
        ivref = var.infinite_variable_ref
        param_nums = InfiniteOpt._parameter_numbers(ivref)
        supp = var.parameter_values
        # find the corresponding variable record the mapping
        vref = lookup_by_support(trans_model, ivref, supp)
        pvref = InfiniteOpt._make_variable_ref(inf_model, index)
        transcription_data(trans_model).finvar_mappings[pvref] = vref
        # update the info constraints as needed
        _update_point_info(pvref, vref)
    end
    return
end

################################################################################
#                       TRANSCRIPTION EXPRESSION METHODS
################################################################################
"""
    transcription_expression(trans_model::JuMP.Model, expr, support::Vector{Float64})

Given the `expr` from an `InfiniteModel`, form its transcripted version in
accordance with the variable mappings available in `trans_model` defined at
`support`. This should only be used once all variables and measures have been
transcribed (e.g., via [`transcribe_hold_variables!`](@ref)).
"""
function transcription_expression(trans_model::JuMP.Model, expr, support)
    error("Unsupported expression type `$(typeof(expr))` for automated " *
          "transcription.")
end

# GeneralVariableRef
function transcription_expression(trans_model::JuMP.Model,
    vref::InfiniteOpt.GeneralVariableRef,
    support::Vector{Float64}
    )
    return transcription_expression(trans_model, vref,
                                    InfiniteOpt._index_type(vref), support)
end

# Infinite variables, infinite parameter functions, and measures
function transcription_expression(trans_model::JuMP.Model,
    vref::InfiniteOpt.GeneralVariableRef,
    index_type::Type{V},
    support::Vector{Float64}
    ) where {V <: Union{InfVarIndex, InfiniteOpt.ParameterFunctionIndex, InfiniteOpt.MeasureIndex}}
    param_nums = InfiniteOpt._parameter_numbers(vref)
    return lookup_by_support(trans_model, vref, index_type, support[param_nums])
end

# Reduced variables
function transcription_expression(trans_model::JuMP.Model,
    vref::InfiniteOpt.GeneralVariableRef,
    index_type::Type{InfiniteOpt.ReducedVariableIndex},
    support::Vector{Float64}
    )
    ivref = InfiniteOpt.infinite_variable_ref(vref)
    if InfiniteOpt._index_type(ivref) == InfiniteOpt.ParameterFunctionIndex
        prefs = InfiniteOpt.raw_parameter_refs(ivref)
        param_nums = InfiniteOpt._parameter_numbers(ivref)
        func = InfiniteOpt.raw_function(ivref)
        return func(Tuple(support[param_nums], prefs)...)
    else
        param_nums = InfiniteOpt._parameter_numbers(vref)
        return lookup_by_support(trans_model, vref, index_type, support[param_nums])
    end
end

# Point variables and hold variables
function transcription_expression(trans_model::JuMP.Model,
    vref::InfiniteOpt.GeneralVariableRef,
    index_type::Type{V},
    support::Vector{Float64}
    )::JuMP.VariableRef where {V <: FinVarIndex}
    return lookup_by_support(trans_model, vref, index_type, support)
end

# Infinite parameters
function transcription_expression(trans_model::JuMP.Model,
    vref::InfiniteOpt.GeneralVariableRef,
    index_type::Type{V},
    support::Vector{Float64}
    )::Float64 where {V <: InfiniteOpt.InfiniteParameterIndex}
    param_num = InfiniteOpt._parameter_number(vref)
    return support[param_num]
end

# Finite parameters
function transcription_expression(trans_model::JuMP.Model,
    vref::InfiniteOpt.GeneralVariableRef,
    index_type::Type{InfiniteOpt.FiniteParameterIndex},
    support::Vector{Float64}
    )::Float64
    return InfiniteOpt.parameter_value(vref)
end

# AffExpr
function transcription_expression(trans_model::JuMP.Model,
    aff::JuMP.GenericAffExpr,
    support::Vector{Float64}
    )
    return JuMP.@expression(trans_model,
                 sum(coef * transcription_expression(trans_model, vref, support)
                     for (coef, vref) in JuMP.linear_terms(aff)) + JuMP.constant(aff))
end

# QuadExpr
function transcription_expression(trans_model::JuMP.Model,
    quad::JuMP.GenericQuadExpr,
    support::Vector{Float64}
    )
    return JuMP.@expression(trans_model,
                 sum(coef * transcription_expression(trans_model, v1, support) *
                     transcription_expression(trans_model, v2, support)
                     for (coef, v1, v2) in JuMP.quad_terms(quad)) +
                 transcription_expression(trans_model, quad.aff, support))
end

# Real Number
function transcription_expression(trans_model::JuMP.Model,
    num::Real,
    support::Vector{Float64}
    )
    return zero(JuMP.AffExpr) + num
end

################################################################################
#                         MEASURE TRANSCRIPTION METHODS
################################################################################
"""
    transcribe_measures!(trans_model::JuMP.Model,
                         inf_model::InfiniteOpt.InfiniteModel)::Nothing

For each `Measure` in `inf_model` expand it via `InfiniteOpt.expand_measure` or
`analytic_expansion` as appropriate and transcribe the expanded expression via
[`transcription_expression`](@ref). Then store the measure to transcripted
expression mappings in `TranscriptionData.measure_mappings` and
`TranscriptionData.measure_lookup` to enable [`transcription_variable`](@ref)
and [`lookup_by_support`](@ref). Note that the supports will not be generated
until `InfiniteOpt.variable_supports` is invoked via `InfiniteOpt.supports`.
Note that `TranscriptionData.measure_support_labels` is also populated.
"""
function transcribe_measures!(trans_model::JuMP.Model,
    inf_model::InfiniteOpt.InfiniteModel
    )::Nothing
    for (index, object) in InfiniteOpt._data_dictionary(inf_model, InfiniteOpt.Measure)
        # get the basic information
        meas = object.measure
        # expand the measure
        if meas.constant_func
            new_expr = InfiniteOpt.analytic_expansion(meas.func, meas.data, trans_model)
        else
            new_expr = InfiniteOpt.expand_measure(meas.func, meas.data, trans_model)
        end
        # prepare to transcribe over the supports
        supp_indices = support_index_iterator(trans_model, meas.object_nums)
        exprs = Vector{JuMP.AbstractJuMPScalar}(undef, length(supp_indices))
        labels = Vector{Set{DataType}}(undef, length(supp_indices))
        lookup_dict = Dict{Vector{Float64}, Int}()
        # map a variable for each support
        for (counter, i) in enumerate(supp_indices)
            raw_supp = index_to_support(trans_model, i)
            @inbounds exprs[counter] = transcription_expression(trans_model,
                                                                new_expr, raw_supp)
            supp = raw_supp[meas.parameter_nums]
            lookup_dict[supp] = counter
            @inbounds labels[counter] = index_to_labels(trans_model, i)
        end
        # save the transcription information
        mref = InfiniteOpt._make_variable_ref(inf_model, index)
        data = transcription_data(trans_model)
        data.measure_lookup[mref] = lookup_dict
        data.measure_mappings[mref] = exprs
        data.measure_support_labels[mref] = labels
    end
    return
end

################################################################################
#                        OBJECTIVE TRANSCRIPTION METHODS
################################################################################
"""
    transcribe_objective!(trans_model::JuMP.Model,
                          inf_model::InfiniteOpt.InfiniteModel)::Nothing

Form the transcripted version of the objective stored in `inf_model` and add it
to `trans_model`. Note that all the variables and measures in `inf_model` must
by transcripted first (e.g., via [`transcribe_infinite_variables!`](@ref)).
"""
function transcribe_objective!(trans_model::JuMP.Model,
    inf_model::InfiniteOpt.InfiniteModel
    )::Nothing
    expr = JuMP.objective_function(inf_model)
    sense = JuMP.objective_sense(inf_model)
    trans_expr = transcription_expression(trans_model, expr, Float64[])
    JuMP.set_objective(trans_model, sense, trans_expr)
    return
end

################################################################################
#                        CONSTRAINT TRANSCRIPTION METHODS
################################################################################
## Given a variable and its set from an info constraint, get the transcribed version
# Lower bound constraint
function _get_info_constr_from_var(trans_model::JuMP.Model,
    vref::InfiniteOpt.GeneralVariableRef,
    set::MOI.GreaterThan,
    support::Vector{Float64}
    )::Union{JuMP.ConstraintRef, Nothing}
    trans_vref = transcription_expression(trans_model, vref, support)
    return JuMP.has_lower_bound(trans_vref) ? JuMP.LowerBoundRef(trans_vref) : nothing
end

# Upper bound constraint
function _get_info_constr_from_var(trans_model::JuMP.Model,
    vref::InfiniteOpt.GeneralVariableRef,
    set::MOI.LessThan,
    support::Vector{Float64}
    )::Union{JuMP.ConstraintRef, Nothing}
    trans_vref = transcription_expression(trans_model, vref, support)
    return JuMP.has_upper_bound(trans_vref) ? JuMP.UpperBoundRef(trans_vref) : nothing
end

# Fix constraint
function _get_info_constr_from_var(trans_model::JuMP.Model,
    vref::InfiniteOpt.GeneralVariableRef,
    set::MOI.EqualTo,
    support::Vector{Float64}
    )::Union{JuMP.ConstraintRef, Nothing}
    trans_vref = transcription_expression(trans_model, vref, support)
    return JuMP.is_fixed(trans_vref) ? JuMP.FixRef(trans_vref) : nothing
end

# Binary constraint
function _get_info_constr_from_var(trans_model::JuMP.Model,
    vref::InfiniteOpt.GeneralVariableRef,
    set::MOI.ZeroOne,
    support::Vector{Float64}
    )::Union{JuMP.ConstraintRef, Nothing}
    trans_vref = transcription_expression(trans_model, vref, support)
    return JuMP.is_binary(trans_vref) ? JuMP.BinaryRef(trans_vref) : nothing
end

# Integer constraint
function _get_info_constr_from_var(trans_model::JuMP.Model,
    vref::InfiniteOpt.GeneralVariableRef,
    set::MOI.Integer,
    support::Vector{Float64}
    )::Union{JuMP.ConstraintRef, Nothing}
    trans_vref = transcription_expression(trans_model, vref, support)
    return JuMP.is_integer(trans_vref) ? JuMP.IntegerRef(trans_vref) : nothing
end

## Extract the parameter bounds from a constr if it has them and collect them
## to return them in a convenient form for iteration
# ScalarConstraint
function _get_constr_bounds(constr::JuMP.ScalarConstraint
    )::Tuple{Vector{Int}, Vector{InfiniteOpt.IntervalSet}}
    return Int[], InfiniteOpt.IntervalSet[]
end

# BoundedScalarConstraint
function _get_constr_bounds(constr::InfiniteOpt.BoundedScalarConstraint
    )::Tuple{Vector{Int}, Vector{InfiniteOpt.IntervalSet}}
    bounds = InfiniteOpt.parameter_bounds(constr)
    prefs = collect(keys(bounds))
    sets = map(p -> bounds[p], prefs)
    indices = map(p -> InfiniteOpt._parameter_number(p), prefs)
    return indices, sets
end

## Determine if a given raw support satisfies constraint parameter bounds
function _support_in_bounds(support::Vector{Float64},
    indices::Vector{Int},
    sets::Vector{InfiniteOpt.IntervalSet}
    )::Bool
    for i in eachindex(indices)
        s = support[indices[i]]
        if !isnan(s) && (s < JuMP.lower_bound(sets[i]) || s > JuMP.upper_bound(sets[i]))
            return false
        end
    end
    return true
end

"""
    transcribe_constraints!(trans_model::JuMP.Model,
                            inf_model::InfiniteOpt.InfiniteModel)::Nothing

For each constraint in `inf_model` form its transcripted version(s) and add them
to `trans_model`. The mappings are stored in `TranscriptionData.constr_mappings`
and the associated supports are stored in `TranscriptionData.constr_supports`
to enable [`transcription_constraint`](@ref) and `InfiniteOpt.constraint_supports`.
Note that variable info constraints are simply mapped to the existing info
constraints already generated along with the transcription variables. Note that
the variables and measures must all first be transcripted (e.g., via
[`transcribe_measures!`](@ref)). Note that
`TranscriptionData.constr_support_labels` is also populated.
"""
function transcribe_constraints!(trans_model::JuMP.Model,
    inf_model::InfiniteOpt.InfiniteModel
    )::Nothing
    param_supps = parameter_supports(trans_model)
    for (index, object) in inf_model.constraints
        # get the basic information
        constr = object.constraint
        func = JuMP.jump_function(constr)
        set = JuMP.moi_set(constr)
        obj_nums = object.object_nums
        # prepare the iteration helpers
        supp_indices = support_index_iterator(trans_model, obj_nums)
        crefs = Vector{JuMP.ConstraintRef}(undef, length(supp_indices))
        supps = Vector{Tuple}(undef, length(supp_indices))
        labels = Vector{Set{DataType}}(undef, length(supp_indices))
        counter = 1
        # iterate over the support indices for the info constraints
        if object.is_info_constraint
            for i in supp_indices
                raw_supp = index_to_support(trans_model, i)
                info_ref = _get_info_constr_from_var(trans_model, func, set,
                                                     raw_supp)
                # not all supports may be defined if overwritten by a point variable
                if info_ref !== nothing
                    @inbounds crefs[counter] = info_ref
                    @inbounds supps[counter] = Tuple(param_supps[j][i[j]]
                                                     for j in obj_nums)
                    @inbounds labels[counter] = index_to_labels(trans_model, i)
                    counter += 1
                end
            end
        # iterate over the supports for regular constraints
        else
            # get basic setup information
            bound_info = _get_constr_bounds(constr)
            name = object.name
            for i in supp_indices
                raw_supp = index_to_support(trans_model, i)
                # ensure the support satisfies parameter bounds and then add it
                if _support_in_bounds(raw_supp, bound_info[1], bound_info[2])
                    new_func = transcription_expression(trans_model, func, raw_supp)
                    trans_constr = JuMP.build_constraint(error, new_func, set)
                    new_name = isempty(name) ? "" : string(name, "(support: ", counter, ")")
                    @inbounds crefs[counter] = JuMP.add_constraint(trans_model,
                                                         trans_constr, new_name)
                    @inbounds supps[counter] = Tuple(param_supps[j][i[j]]
                                                     for j in obj_nums)
                    @inbounds labels[counter] = index_to_labels(trans_model, i)
                    counter += 1
                end
            end
        end
        # truncate the arrays in case not all the supports satisfied the bounds
        deleteat!(crefs, counter:length(crefs))
        deleteat!(supps, counter:length(supps))
        deleteat!(labels, counter:length(supps))
        # add the constraint mappings to the trans model
        cref = InfiniteOpt._make_constraint_ref(inf_model, index)
        data = transcription_data(trans_model)
        data.constr_mappings[cref] = crefs
        data.constr_supports[cref] = supps
        data.constr_support_labels[cref] = labels
    end
    return
end

################################################################################
#                      DERIVATIVE CONSTRAINT TRANSCRIPTION METHODS
################################################################################
"""
    transcribe_derivative_evaluations!(trans_model::JuMP.Model,
                                       inf_model::InfiniteOpt.InfiniteModel)::Nothing

Generate the auxiliary derivative evaluation equations and transcribe them
appropriately for all the derivatives in `inf_model`. These are in turn added to
`trans_model`. Note that no mapping information is recorded since the InfiniteModel
won't have any constraints that correspond to these equations. Also Note that
the variables and measures must all first be transcripted (e.g., via
[`transcribe_derivative_variables!`](@ref)).
"""
function transcribe_derivative_evaluations!(trans_model::JuMP.Model,
    inf_model::InfiniteOpt.InfiniteModel
    )::Nothing
    for (index, object) in InfiniteOpt._data_dictionary(inf_model, InfiniteOpt.Derivative)
        # get the basic variable information
        dref = InfiniteOpt._make_variable_ref(inf_model, index)
        pref = dispatch_variable_ref(object.variable.parameter_ref)
        method = InfiniteOpt.derivative_method(pref)
        # generate the transcription constraints as needed
        if !InfiniteOpt.has_derivative_constraints(dref)
            # generate the evaluation expressions
            exprs = InfiniteOpt.evaluate_derivative(dref, method, trans_model)
            # prepare the iteration helpers
            param_obj_num = InfiniteOpt._object_number(pref)
            obj_nums = filter(!isequal(param_obj_num), InfiniteOpt._object_numbers(dref))
            supp_indices = support_index_iterator(trans_model, obj_nums)
            # transcribe the constraints
            set = MOI.EqualTo(0.0)
            for i in supp_indices
                raw_supp = index_to_support(trans_model, i)
                for expr in exprs
                    new_expr = transcription_expression(trans_model, expr, raw_supp)
                    trans_constr = JuMP.build_constraint(error, new_expr, set)
                    JuMP.add_constraint(trans_model, trans_constr) # TODO maybe add name?
                end
            end
        end
    end
    return
end

################################################################################
#                      INFINITEMODEL TRANSCRIPTION METHODS
################################################################################
"""
    build_transcription_model!(trans_model::JuMP.Model,
                               inf_model::InfiniteOpt.InfiniteModel;
                               [check_support_dims::Bool = true])::Nothing

Given an empty `trans_model` build it using the information stored in `inf_model`.
This is intended for a `TranscriptionModel` that serves as a internal optimizer model
of `inf_model`. This detail is important to correctly enable internally generated
reduced variables during the transcription process such that `inf_model` is not
modified. Note that this will add supports to `inf_model` via
[`InfiniteOpt.fill_in_supports!`](@ref) for infinite parameters that contain
no supports. Also a warning is thrown when the transcription model contains
more than 15,000 support points to alert users when they may naively have
a few independent supports whose product quickly yields a very large grid.
For example having 3 independent parameters with 100 supports each would result
in 1,000,000 supports if all three are together in at least 1 constraint. This
behavior can be overcome using dependent parameters. The warning can be turned off
via `check_support_dims = false`.
"""
function build_transcription_model!(trans_model::JuMP.Model,
    inf_model::InfiniteOpt.InfiniteModel; check_support_dims::Bool = true
    )::Nothing
    # ensure there are supports to add and add them to the trans model
    InfiniteOpt.fill_in_supports!(inf_model, modify = false)
    set_parameter_supports(trans_model, inf_model)
    # check that there isn't a crazy amount of supports from taking the product
    supps = parameter_supports(trans_model)
    num_supps = prod(length, supps)
    if check_support_dims && length(supps) > 1 && num_supps > 15000 # NOTE this is an arbitrary cutoff
        @warn("Due to necessarily considering the combinatorics of independent " *
              "parameter supports, the model will be transcripted over up to $(num_supps) " *
              "supports (if any constraint uses all the infinite parameters) "  *
              "and thus naive solution of the discretized problem may be slow. " *
              "This warning can be turned off via `check_support_dims = false`.")
    end
    # define the variables
    transcribe_hold_variables!(trans_model, inf_model)
    transcribe_infinite_variables!(trans_model, inf_model)
    transcribe_derivative_variables!(trans_model, inf_model)
    transcribe_reduced_variables!(trans_model, inf_model)
    transcribe_point_variables!(trans_model, inf_model)
    transcribe_measures!(trans_model, inf_model)
    # define the objective
    transcribe_objective!(trans_model, inf_model)
    # define the constraints
    transcribe_constraints!(trans_model, inf_model)
    # define the derivative evaluation constraints
    transcribe_derivative_evaluations!(trans_model, inf_model)
    return
end

# TODO make TranscriptionModel(::InfiniteModel) --> have to deal with accessing reduced vars

ext: jl


text:
"""
    SolidMaterial(density=3000.0, YoungsModulus=2e11)

Generates a SolidMaterial object by providing the material properties of a [Solid] with keyword arguments.

# Arguments
Arguments that are not provided have value `NaN`.
- `density` in [kg/m^3]: Density, see [Wikipedia](https://en.wikipedia.org/wiki/Density).
- `YoungsModulus` in [Pa]: Youngs's modulus, see [Wikipedia](https://en.wikipedia.org/wiki/Young%27s_modulus).
- `PoissonsRatio`: Poisson's ratio, see [Wikipedia](https://en.wikipedia.org/wiki/Poisson%27s_ratio).
- `meltingPoint` in [K]: Melting point, see [Wikipedia](https://en.wikipedia.org/wiki/Melting_point).
                         If the material is destroyed before its melting point (e.g. wood that is burning)
                        then `meltingPoint` is the temperature when destruction of the solid starts.
- `specificHeatCapacity` in [J/(kg.K)]: Specific heat capacity, see [Wikipedia](https://en.wikipedia.org/wiki/Heat_capacity).
- `thermalConductivity` in [W/(m.K)]: Thermal conductivity, see [Wikipedia](https://en.wikipedia.org/wiki/Thermal_conductivity) and
                                [List of thermal conductivities](https://en.wikipedia.org/wiki/List_of_thermal_conductivities)
- `linearThermalExpansionCoefficient::Float64` in [1/K]: Linear thermal expansion coefficient, see [Wikipedia](https://en.wikipedia.org/wiki/Thermal_expansion).

"""
mutable struct SolidMaterial
   density::Float64              # [kg/m^3], https://en.wikipedia.org/wiki/Density
   YoungsModulus::Float64        # [Pa], https://en.wikipedia.org/wiki/Young%27s_modulus
   PoissonsRatio::Float64        # [],   https://en.wikipedia.org/wiki/Poisson%27s_ratio
   meltingPoint::Float64         # [K],  https://en.wikipedia.org/wiki/Melting_point
                                 #       If the material is destroyed before its melting point
                                 #       (e.g. wood that is burning) then meltingPoint is the
                                 #       temperature when destruction starts.
   specificHeatCapacity::Float64 # [J/(kg.K)], https://en.wikipedia.org/wiki/Heat_capacity
   thermalConductivity::Float64  # [W/(m.K)],  https://en.wikipedia.org/wiki/Thermal_conductivity
                                 #             https://en.wikipedia.org/wiki/List_of_thermal_conductivities
   linearThermalExpansionCoefficient::Float64  # [1/K], https://en.wikipedia.org/wiki/Thermal_expansion
end
SolidMaterial(; density=NaN,
                 YoungsModulus=NaN,
                 PoissonsRatio=NaN,
                 meltingPoint=NaN,
                 specificHeatCapacity=NaN,
                 thermalConductivity=NaN,
                 linearThermalExpansionCoefficient=NaN) =
              SolidMaterial(density, YoungsModulus, PoissonsRatio, meltingPoint, specificHeatCapacity,
                             thermalConductivity, linearThermalExpansionCoefficient)

"""
    readSolidMaterialFromJSON(fileName)

Read a JSON file consisting of a dictionary of SolidMaterial instances from `fileName` and
return a `Dict{String, SolidMaterial}` dictionary.
"""
readSolidMaterialFromJSON(fileName::AbstractString) =
         Basics.readDictOfStructsFromJSON(fileName, SolidMaterial)

"""
    const solidMaterialPalette

Dictionary of solid material data, see [`SolidMaterial`](@ref)
"""
const solidMaterialPalette =
         readSolidMaterialFromJSON( joinpath(Modia3D.path, "palettes", "solidMaterials.json") )


text:
using NNlib: DenseConvDims
import NNlib: conv!, ∇conv_filter!, ∇conv_data!, conv_bias_act!

using CUDA.CUDNN: scalingParameter, CUDNN_CONVOLUTION, convdims,
                  cudnnConvolutionDescriptor, cudnnConvolutionBwdDataAlgoPerf,
                  cudnnConvolutionForward!, cudnnConvolutionBwdFilterAlgoPerf,
                  cudnnConvolutionBackwardData, cudnnConvolutionBackwardFilter,
                  cudnnConvolutionBackwardBias

const CUDNNFloat = Union{Float16,Float32,Float64}

function cudnnConvolutionDescriptorAndPaddedInput(cdims::DenseConvDims, x::DenseCuArray{T}) where T
    # The main purpose of this function is to catch asymmetric padding which cudnn does not support
    # If we find asymmetric padding we'll make a copy of x which is manually padded so that we can
    # call cudnn with symmetric padding.
    pad = NNlib.padding(cdims)
    sdims = NNlib.spatial_dims(cdims)
    all(i -> pad[i] .== pad[i+1], 1:2:2sdims) && return (cudnnConvolutionDescriptor(cdims, x), x, identity)

    # Naive implementation, is there a faster way?
    # How much we need to pad x manually: The absolute difference between pad_left and pad_right, pad_top
    # and pad_bottom etc. respectively. We keep the sign here though because we use it below to figure out
    # which side of x to pad. Oh, and we use a CartesianIndex as we will mainly use this to index in x
    pad_manual = CartesianIndex(ntuple(i -> i > sdims ? 0 : pad[2(i-1)+1] - pad[2(i-1)+2], ndims(x)))

    # How much we can let cudnn pad: The smallest padding amount between pad_left and pad_right, pad_top
    # and pad_bottom etc. respectively
    pad_cudnn = ntuple(i -> min(pad[2(i-1)+1], pad[2(i-1)+2]), sdims)

    x_padded_size = ntuple(i -> i <= sdims ? size(x, i) + abs(pad_manual[i]) : size(x ,i), ndims(x))
    x_padded = similar(x, x_padded_size)
    fill!(x_padded, 0)
    # This is a bit yucky, but we are basically figuring out where in x_padded we shall insert x
    # Haven't benchmarked if this has any advantages over a more readable solution, e.g. writing dim
    # by dim to an array in a loop
    xIs = CartesianIndices(x)
    xI_first = first(xIs)
    xI_last = last(xIs)
    xIs_pad = max(xI_first, xI_first + pad_manual) : max(xI_last, xI_last + pad_manual)
    x_padded[xIs_pad] = x

    return cudnnConvolutionDescriptor(cdims, x_padded, pad_cudnn), x_padded, _x -> _x[xIs_pad]
end

function cudnnConvolutionDescriptor(cdims::DenseConvDims, x::DenseCuArray{T}, pad = nnlibPadding(cdims)) where T
    mode=(NNlib.flipkernel(cdims) ? CUDNN_CROSS_CORRELATION : CUDNN_CONVOLUTION)
    cudnnConvolutionDescriptor(convdims(pad, size(x),0),
                               convdims(NNlib.stride(cdims),size(x),1),
                               convdims(NNlib.dilation(cdims),size(x),1),
                               mode,
                               cudnnDataType(T),
                               math_mode(),
                               CUDNN_DEFAULT_REORDER,
                               Cint(NNlib.groupcount(cdims)))
end

function conv!(y::DenseCuArray{T}, x::DenseCuArray{T}, w::DenseCuArray{T}, cdims::DenseConvDims;
               alpha=1, beta=0, algo=-1) where T<:CUDNNFloat
    if cudnnversion() < v"6"
        all(x -> x == 1, dilation(cdims)) || error("Only dilation = 1 is supported in cuDNN version < 6")
    end
    if algo != -1
        @warn "algo option has been deprecated, the fastest algo is computed automatically" maxlog=1
    end
    d, x, _ = cudnnConvolutionDescriptorAndPaddedInput(cdims, x)
    cudnnConvolutionForward!(y, w, x, d; alpha, beta, z=y)
end

function conv_bias_act!(y::DenseCuArray{T}, x::DenseCuArray{T}, w::DenseCuArray{T},
                        cdims::DenseConvDims, bias::DenseCuArray{T}, σ=identity;
                        z::DenseCuArray{T}=y, alpha=1, beta=0, algo=-1) where T<:CUDNNFloat
    if cudnnversion() < v"6"
        all(x -> x == 1, dilation(cdims)) || error("Only dilation = 1 is supported in cuDNN version < 6")
    end
    if algo != -1
        @warn "The algo option has been deprecated, the fastest algo is computed automatically" maxlog=1
    end
    d, x, _ = cudnnConvolutionDescriptorAndPaddedInput(cdims, x)
    # only relu and identity are supported by cudnnConvolutionForward!
    activation = (σ == NNlib.relu ? CUDNN_ACTIVATION_RELU : CUDNN_ACTIVATION_IDENTITY)
    cudnnConvolutionForward!(y, w, x, d; z, bias, activation, alpha, beta)
    if activation === CUDNN_ACTIVATION_IDENTITY && σ ∉ (nothing, identity)
        y = σ.(y)
    end
    return y
end

function ∇conv_data!(dx::DenseCuArray{T}, dy::DenseCuArray{T}, w::DenseCuArray{T},
                     cdims::DenseConvDims; alpha=1, beta=0, algo=-1) where T<:CUDNNFloat
    if cudnnversion() < v"6"
        all(x -> x == 1, dilation(cdims)) || error("Only dilation = 1 is supported in cuDNN version < 6")
    end
    if algo != -1
        @warn "The algo option has been deprecated, the fastest algo is computed automatically" maxlog=1
    end
    alpha, beta = scalingParameter(T,alpha), scalingParameter(T,beta);
    convDesc, dx, depad = cudnnConvolutionDescriptorAndPaddedInput(cdims, dx)
    xDesc, yDesc, wDesc = cudnnTensorDescriptor(dx), cudnnTensorDescriptor(dy), cudnnFilterDescriptor(w)
    p = cudnnConvolutionBwdDataAlgoPerf(wDesc, w, yDesc, dy, convDesc, xDesc, dx)
    with_workspace(p.memory) do workspace
        cudnnConvolutionBackwardData(handle(), alpha, wDesc, w, yDesc, dy, convDesc, p.algo, workspace, sizeof(workspace), beta, xDesc, dx)
    end
    return depad(dx)
end

function ∇conv_filter!(dw::DenseCuArray{T}, x::DenseCuArray{T}, dy::DenseCuArray{T},
                       cdims::DenseConvDims; alpha=1, beta=0, algo=-1) where T<:CUDNNFloat
    if cudnnversion() < v"6"
        all(x -> x == 1, dilation(cdims)) || error("Only dilation = 1 is supported in cuDNN version < 6")
    end
    if algo != -1
        @warn "The algo option has been deprecated, the fastest algo is computed automatically" maxlog=1
    end
    alpha, beta = scalingParameter(T,alpha), scalingParameter(T,beta);
    convDesc, x, _ = cudnnConvolutionDescriptorAndPaddedInput(cdims, x)
    xDesc, yDesc, wDesc = cudnnTensorDescriptor(x), cudnnTensorDescriptor(dy), cudnnFilterDescriptor(dw)
    p = cudnnConvolutionBwdFilterAlgoPerf(xDesc, x, yDesc, dy, convDesc, wDesc, dw);
    with_workspace(p.memory) do workspace
        cudnnConvolutionBackwardFilter(handle(), alpha, xDesc, x, yDesc, dy, convDesc, p.algo, workspace, sizeof(workspace), beta, wDesc, dw);
    end
    return dw
end

function ∇conv_bias!(db::DenseCuArray{T}, dy::DenseCuArray{T}; alpha=1, beta=0) where T<:CUDNNFloat
    alpha,beta = scalingParameter(T,alpha), scalingParameter(T,beta)
    bDesc, yDesc = cudnnTensorDescriptor.((db,dy))
    cudnnConvolutionBackwardBias(handle(), alpha, yDesc, dy, beta, bDesc, db)
    return db
end


text:
struct Spambase <: Name end

# mandatory methods
format(::Type{Spambase}) = TabularData
task(::Type{Spambase}) = TwoClass
source(::Type{Spambase}) = "http://archive.ics.uci.edu/ml/datasets/spambase"

function downloadlink(::Type{Spambase})
    return "https://archive.ics.uci.edu/ml/machine-learning-databases/spambase/spambase.data"
end

function preprocess(N::Type{Spambase})
    return path -> preprocess_csv(N, path, :train; istarget = 58)
end

nattributes(::Type{Spambase}) = (57, )
ninstances(::Type{Spambase}) = (4601, 0, 0)
positivelabel(::Type{Spambase}) = 1

# optional methods
function checksum(::Type{Spambase})
    return "b1ef93de71f97714d3d7d4f58fc9f718da7bbc8ac8a150eff2778616a8097b12"
end

name(::Type{Spambase}) = "Spam e-mail database"
function author(::Type{Spambase})
    return ["Mark Hopkins", "Erik Reeber", "George Forman", "Jaap Suermondt"]
end

function licence(::Type{Spambase})
    return """
    This data set is licensed under a Creative Commons Attribution 4.0 International (CC BY 4.0) license. This allows for the sharing and adaptation of the datasets for any purpose, provided that the appropriate credit is given.

    URL: https://creativecommons.org/licenses/by/4.0/legalcode
    """
end

function citation(::Type{Spambase})
    return """
    @misc{misc_spambase_94,
        author       = {Hopkins, Mark, Reeber, Erik, Forman, George & Suermondt, Jaap},
        title        = {{Spambase}},
        year         = {1999},
        howpublished = {UCI Machine Learning Repository}
    }
    """
end

max_stars_repo_name: VaclavMacha/DatasetProvider

max_issues_repo_name: VaclavMacha/DatasetProvider

max_forks_repo_name: VaclavMacha/DatasetProvider


text:
#!/bin/julia
# What Would Hamilton Do? (WWHD?)
using Quaternions

# Origin of MA in Fract Coords
origin=[0.5 0.5 0.5]' # origin of MA in unit cell, fractional coords

POSCARlines = readlines(open("POSCAR"))

shift=[0.05*randn() 0.05*randn() 0.05*randn()]'

# Generate a normalized Quaternion with normally distributed random 4 component
qr=normalize(Quaternion(randn(),randn(),randn(),randn()))
# This should be an evenly distributed rotation matrix
rotate=rotationmatrix(qr)

#rotate=eye(3) #3-dimensional identity matrix; null rotation operation

#println("Determinant (should be 1): ",det(rotate))
#println("Does Q^TQ = I? \n",rotate*rotate')

for POSCARline in enumerate(POSCARlines)
    # Start line 10 if Selective dynamics; 9 otherwise
    if (POSCARline[1]>=9 && POSCARline[1]<=9+8) #if MA; for Cubic Perov structs.
        r=float(split(POSCARline[2])) # "1.2 3.4 5.6" --> [1.2; 3.4; 5.6]
        r=shift+origin+(rotate*(r-origin)) # apply shift + rotate, shift back
        r=mod.(r,1) # Modulo arithmatic so that all {x,y,z} values are on [0,1]
        @printf "  %.16f  %.16f  %.16f \n" r[1] r[2] r[3]
        # 16 digits of precision seems to match VASP output
    else # otherwise echo line unchanged to output
        println(POSCARline[2])
    end
end

max_stars_repo_name: jarvist/Julia-Phonons

max_issues_repo_name: lucydot/Julia-Phonons

max_forks_repo_name: lucydot/Julia-Phonons


text:
using BinaryBuilder

# Collection of sources required to build SuiteSparse
sources = [
    "http://faculty.cse.tamu.edu/davis/SuiteSparse/SuiteSparse-5.2.0.tar.gz" =>
    "3c46c035ea8217649958a0f73360e825b0c9dcca4e32a9349d2c7678c0d48813",
    "./src",
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/SuiteSparse

UNAME=`uname`
LDFLAGS="-L$prefix/lib"

if [[ ${UNAME} == MSYS_NT-6.3 ]]; then
    UNAME="Windows"
    LDFLAGS="-L$prefix/lib -shared"
    CFOPENMP=
fi

if [[ ${nbits} == 64 ]]; then
    BLAS="-lopenblas64_"
    LAPACK="-lopenblas64_"
    UMFPACK_CONFIG="-DSUN64 -DLONGBLAS='long long'"
    CHOLMOD_CONFIG="-DSUN64 -DLONGBLAS='long long' -DNPARTITION"
    SPQR_CONFIG="-DSUN64 -DLONGBLAS='long long'"
else
    BLAS="-lopenblas"
    LAPACK="-lopenblas"
    UMFPACK_CONFIG=
    CHOLMOD_CONFIG="-DNPARTITION"
    SPQR_CONFIG=
fi

make -j -C SuiteSparse_config library INSTALL="$prefix" BLAS="$BLAS" LAPACK="$LAPACK" UMFPACK_CONFIG="$UMFPACK_CONFIG" CHOLMOD_CONFIG="$CHOLMOD_CONFIG" SPQR_CONFIG="$SPQR_CONFIG" UNAME="$UNAME" LDFLAGS="$LDFLAGS" CFOPENMP="$CFOPENMP" config

for proj in SuiteSparse_config AMD BTF CAMD CCOLAMD COLAMD CHOLMOD LDL KLU UMFPACK RBio SPQR; do
    make -j -C $proj library INSTALL="$prefix" BLAS="$BLAS" LAPACK="$LAPACK" UMFPACK_CONFIG="$UMFPACK_CONFIG" CHOLMOD_CONFIG="$CHOLMOD_CONFIG" SPQR_CONFIG="$SPQR_CONFIG" UNAME="$UNAME" LDFLAGS="$LDFLAGS" CFOPENMP="$CFOPENMP"
    make -j -C $proj install INSTALL="$prefix" BLAS="$BLAS" LAPACK="$LAPACK" UMFPACK_CONFIG="$UMFPACK_CONFIG" CHOLMOD_CONFIG="$CHOLMOD_CONFIG" SPQR_CONFIG="$SPQR_CONFIG" UNAME="$UNAME" LDFLAGS="$LDFLAGS" CFOPENMP="$CFOPENMP"
echo make "$MAKE_OPTS"
done

# For now, we'll have to adjust the name of the OpenBLAS library on macOS.
# Eventually, this should be fixed upstream
if [[ ${target} == "x86_64-apple-darwin14" ]]; then
    echo "-- Modifying library name for OpenBLAS"

    for nm in libcholmod.3.0.12 libspqr.2.0.8 libumfpack.5.7.6; do
        install_name_tool -change libopenblas64_.0.2.20.dylib @rpath/libopenblas64_.dylib ${prefix}/lib/${nm}.dylib
    done
fi

# Compile suitesparse_wrapper shim
cd $WORKSPACE/srcdir/SuiteSparse_wrapper
make install

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    BinaryProvider.Windows(:i686, :blank_libc, :blank_abi),
    BinaryProvider.Windows(:x86_64, :blank_libc, :blank_abi),
    BinaryProvider.Linux(:x86_64, :glibc, :blank_abi),
    BinaryProvider.Linux(:i686, :glibc, :blank_abi),
    BinaryProvider.Linux(:aarch64, :glibc, :blank_abi),
    BinaryProvider.Linux(:armv7l, :glibc, :eabihf),
    BinaryProvider.MacOS(:x86_64, :blank_libc, :blank_abi),
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libsuitesparseconfig",   :libsuitesparseconfig),
    LibraryProduct(prefix, "libamd",                 :libamd),
    LibraryProduct(prefix, "libbtf",                 :libbtf),
    LibraryProduct(prefix, "libcamd",                :libcamd),
    LibraryProduct(prefix, "libccolamd",             :libccolamd),
    LibraryProduct(prefix, "libcolamd",              :libcolamd),
    LibraryProduct(prefix, "libcholmod",             :libcholmod),
    LibraryProduct(prefix, "libldl",                 :libldl),
    LibraryProduct(prefix, "libklu",                 :libklu),
    LibraryProduct(prefix, "libumfpack",             :libumfpack),
    LibraryProduct(prefix, "librbio",                :librbio),
    LibraryProduct(prefix, "libspqr",                :libspqr),
    LibraryProduct(prefix, "libsuitesparse_wrapper", :libsuitesparse_wrapper),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    "https://github.com/staticfloat/OpenBLASBuilder/releases/download/v0.2.20-7/build.jl"
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, "SuiteSparse", sources, script, platforms, products, dependencies)


text:
# Unit tests for Distances

@testset "individual metrics" begin

a = 1
b = 2
@test sqeuclidean(a, a) == 0.
@test sqeuclidean(a, b) == 1.

@test euclidean(a, a) == 0.
@test euclidean(a, b) == 1.

@test cityblock(a, a) == 0.
@test cityblock(a, b) == 1.

@test chebyshev(a, a) == 0.
@test chebyshev(a, b) == 1.

@test chebyshev(a, a) == 0.

@test minkowski(a, a, 2) == 0.
@test minkowski(a, b, 2) == 1.

@test hamming(a, a) == 0
@test hamming(a, b) == 1

bt = [true, false, true]
bf = [false, true, true]
@test rogerstanimoto(bt, bt) == 0
@test rogerstanimoto(bt, bf) == 4./5

p = rand(12)
p[p .< 0.3] = 0.
q = rand(12)
a = [1., 2., 1., 3., 2., 1.]
b = [1., 3., 0., 2., 2., 0.]
for (x, y) in (([4., 5., 6., 7.], [3., 9., 8., 1.]),
                ([4., 5., 6., 7.], [3. 8.; 9. 1.]))
    @test sqeuclidean(x, x) == 0.
    @test sqeuclidean(x, y) == 57.

    @test euclidean(x, x) == 0.
    @test euclidean(x, y) == sqrt(57.)

    @test jaccard(x, x) == 0
    @test jaccard(x, y) == 13./28

    @test cityblock(x, x) == 0.
    @test cityblock(x, y) == 13.

    @test chebyshev(x, x) == 0.
    @test chebyshev(x, y) == 6.

    @test minkowski(x, x, 2) == 0.
    @test minkowski(x, y, 2) == sqrt(57.)

    @test cosine_dist(x, x) ≈ 0.0

    @test_throws DimensionMismatch cosine_dist(1.:2, 1.:3)
    @test cosine_dist(x, y) ≈ (1.0 - 112. / sqrt(19530.))

    @test corr_dist(x, x) < 1.0e-14
    @test corr_dist(x, y) ≈ cosine_dist(x .- mean(x), vec(y) .- mean(y))

    @test chisq_dist(x, x) == 0.
    @test chisq_dist(x, y) == sum((x - vec(y)).^2 ./ (x + vec(y)))

    klv = 0.
    for i = 1 : length(p)
        if p[i] > 0
            klv += p[i] * log(p[i] / q[i])
        end
    end
    @test kl_divergence(p, q) ≈ klv

    pm = (p + q) / 2
    jsv = kl_divergence(p, pm) / 2 + kl_divergence(q, pm) / 2
    @test js_divergence(p, p) ≈ 0.0
    @test js_divergence(p, q) ≈ jsv

    @test spannorm_dist(x, x) == 0.
    @test spannorm_dist(x, y) == maximum(x - vec(y)) - minimum(x - vec(y))

    w = ones(4)
    @test sqeuclidean(x, y) ≈ wsqeuclidean(x, y, w)

    w = rand(size(x))

    @test wsqeuclidean(x, x, w) == 0.
    @test wsqeuclidean(x, y, w) ≈ dot((x - vec(y)).^2, w)

    @test weuclidean(x, x, w) == 0.
    @test weuclidean(x, y, w) == sqrt(wsqeuclidean(x, y, w))

    @test wcityblock(x, x, w) == 0.
    @test wcityblock(x, y, w) ≈ dot(abs(x - vec(y)), w)

    @test wminkowski(x, x, w, 2) == 0.
    @test wminkowski(x, y, w, 2) ≈ weuclidean(x, y, w)

    w = rand(size(a))

    @test whamming(a, a, w) == 0.
    @test whamming(a, b, w) == sum((a .!= b) .* w)
end

end # testset

@testset "NaN behavior" begin

a = [NaN, 0]; b = [0, NaN]
@test isnan(chebyshev(a, b)) == isnan(maximum(a-b))
a = [NaN, 0]; b = [0, 1]
@test isnan(chebyshev(a, b)) == isnan(maximum(a-b))

end #testset

@testset "empty vector" begin

a = Float64[]
b = Float64[]
@test sqeuclidean(a, b) == 0.
@test isa(sqeuclidean(a, b), Float64)
@test euclidean(a, b) == 0.
@test isa(euclidean(a, b), Float64)
@test cityblock(a, b) == 0.
@test isa(cityblock(a, b), Float64)
@test chebyshev(a, b) == 0
@test isa(chebyshev(a, b), Float64)
@test minkowski(a, b, 2) == 0.
@test isa(minkowski(a, b, 2), Float64)
@test hamming(a, b) == 0.0
@test isa(hamming(a, b), Int)

w = Float64[]
@test isa(whamming(a, b, w), Float64)

end # testset

@testset "DimensionMismatch throwing" begin

a = [1, 0]; b = [2]
@test_throws DimensionMismatch sqeuclidean(a, b)
a = [1, 0]; b = [2.0] ; w = [3.0]
@test_throws DimensionMismatch wsqeuclidean(a, b, w)
a = [1, 0]; b = [2.0, 4.0] ; w = [3.0]
@test_throws DimensionMismatch wsqeuclidean(a, b, w)

end # testset

@testset "mahalanobis" begin

x, y = [4., 5., 6., 7.], [3., 9., 8., 1.]
a = [1., 2., 1., 3., 2., 1.]
b = [1., 3., 0., 2., 2., 0.]

Q = rand(length(x), length(x))
Q = Q * Q'  # make sure Q is positive-definite
@test sqmahalanobis(x, x, Q) == 0.
@test sqmahalanobis(x, y, Q) ≈ dot(x - y, Q * (x - y))

@test mahalanobis(x, x, Q) == 0.
@test mahalanobis(x, y, Q) == sqrt(sqmahalanobis(x, y, Q))

end #testset

@testset "bhattacharyya / hellinger" begin

x, y = [4., 5., 6., 7.], [3., 9., 8., 1.]
a = [1., 2., 1., 3., 2., 1.]
b = [1., 3., 0., 2., 2., 0.]
p = rand(12)
p[p .< 0.3] = 0.
q = rand(12)

# Bhattacharyya and Hellinger distances are defined for discrete
# probability distributions so to calculate the expected values
# we need to normalize vectors.
px = x ./ sum(x)
py = y ./ sum(y)
expected_bc_x_y = sum(sqrt(px .* py))
@test Distances.bhattacharyya_coeff(x, y) ≈ expected_bc_x_y
@test bhattacharyya(x, y) ≈ (-log(expected_bc_x_y))
@test hellinger(x, y) ≈ sqrt(1 - expected_bc_x_y)

pa = a ./ sum(a)
pb = b ./ sum(b)
expected_bc_a_b = sum(sqrt(pa .* pb))
@test Distances.bhattacharyya_coeff(a, b) ≈ expected_bc_a_b
@test bhattacharyya(a, b) ≈ (-log(expected_bc_a_b))
@test hellinger(a, b) ≈ sqrt(1 - expected_bc_a_b)

pp = p ./ sum(p)
pq = q ./ sum(q)
expected_bc_p_q = sum(sqrt(pp .* pq))
@test Distances.bhattacharyya_coeff(p, q) ≈ expected_bc_p_q
@test bhattacharyya(p, q) ≈ (-log(expected_bc_p_q))
@test hellinger(p, q) ≈ sqrt(1 - expected_bc_p_q)

# Ensure it is semimetric
@test bhattacharyya(x, y) ≈ bhattacharyya(y, x)

end #testset

macro test_colwise(dist, x, y)
    quote
        local n = size($x, 2)
        r1 = zeros(n)
        r2 = zeros(n)
        r3 = zeros(n)
        for j = 1 : n
            r1[j] = evaluate($dist, ($x)[:,j], ($y)[:,j])
            r2[j] = evaluate($dist, ($x)[:,1], ($y)[:,j])
            r3[j] = evaluate($dist, ($x)[:,j], ($y)[:,1])
        end
        @test colwise($dist, $x, $y) ≈ r1
        @test colwise($dist, ($x)[:,1], $y) ≈ r2
        @test colwise($dist, $x, ($y)[:,1]) ≈ r3
    end
end

@testset "column-wise metrics" begin

m = 5
n = 8
X = rand(m, n)
Y = rand(m, n)
A = rand(1:3, m, n)
B = rand(1:3, m, n)

P = rand(m, n)
Q = rand(m, n)
P[P .< 0.3] = 0.

@test_colwise SqEuclidean() X Y
@test_colwise Euclidean() X Y
@test_colwise Cityblock() X Y
@test_colwise Chebyshev() X Y
@test_colwise Minkowski(2.5) X Y
@test_colwise Hamming() A B

@test_colwise CosineDist() X Y
@test_colwise CorrDist() X Y

@test_colwise ChiSqDist() X Y
@test_colwise KLDivergence() P Q
@test_colwise JSDivergence() P Q
@test_colwise SpanNormDist() X Y

@test_colwise BhattacharyyaDist() X Y
@test_colwise HellingerDist() X Y

w = rand(m)

@test_colwise WeightedSqEuclidean(w) X Y
@test_colwise WeightedEuclidean(w) X Y
@test_colwise WeightedCityblock(w) X Y
@test_colwise WeightedMinkowski(w, 2.5) X Y
@test_colwise WeightedHamming(w) A B

Q = rand(m, m)
Q = Q * Q'  # make sure Q is positive-definite

@test_colwise SqMahalanobis(Q) X Y
@test_colwise Mahalanobis(Q) X Y

end # testset

macro test_pairwise(dist, x, y)
    quote
        local nx = size($x, 2)
        local ny = size($y, 2)
        rxy = zeros(nx, ny)
        rxx = zeros(nx, nx)
        for j = 1 : ny, i = 1 : nx
            rxy[i, j] = evaluate($dist, ($x)[:,i], ($y)[:,j])
        end
        for j = 1 : nx, i = 1 : nx
            rxx[i, j] = evaluate($dist, ($x)[:,i], ($x)[:,j])
        end
        @test pairwise($dist, $x, $y) ≈ rxy
        @test pairwise($dist, $x) ≈ rxx
    end
end

@testset "pairwise metrics" begin

m = 5
n = 8
nx = 6
ny = 8

X = rand(m, nx)
Y = rand(m, ny)
A = rand(1:3, m, nx)
B = rand(1:3, m, ny)

P = rand(m, nx)
Q = rand(m, ny)

@test_pairwise SqEuclidean() X Y
@test_pairwise Euclidean() X Y
@test_pairwise Cityblock() X Y
@test_pairwise Chebyshev() X Y
@test_pairwise Minkowski(2.5) X Y
@test_pairwise Hamming() A B

@test_pairwise CosineDist() X Y
@test_pairwise CorrDist() X Y

@test_pairwise ChiSqDist() X Y
@test_pairwise KLDivergence() P Q
@test_pairwise JSDivergence() P Q

@test_pairwise BhattacharyyaDist() X Y
@test_pairwise HellingerDist() X Y

w = rand(m)

@test_pairwise WeightedSqEuclidean(w) X Y
@test_pairwise WeightedEuclidean(w) X Y
@test_pairwise WeightedCityblock(w) X Y
@test_pairwise WeightedMinkowski(w, 2.5) X Y
@test_pairwise WeightedHamming(w) A B

Q = rand(m, m)
Q = Q * Q'  # make sure Q is positive-definite

@test_pairwise SqMahalanobis(Q) X Y
@test_pairwise Mahalanobis(Q) X Y

end #testset


text:
module Delaunay

using PyCall
using SparseArrays

export Triangulation, delaunay

const spatial = PyNULL()

function __init__()
    return copy!(spatial, pyimport_conda("scipy.spatial", "scipy"))
end

"""
   Triangulation

Result of Delaunay triangulation.
"""
struct Triangulation
    "Coordinates of input points (npoints, ndim)"
    points::Array{Float64,2}
    "(One-based) indices of the points forming the simplices (nsimplex, ndim+1)"
    simplices::Array{Int,2}
    "Indices of neighbor simplices for each simplex; 0 at boundary (nsimplex, ndim+1)"
    neighbors::Array{Int,2}
    "[normal, offset] forming the hyperplane equation of the facet on the paraboloid (nsimplex, ndim+2)"
    equations::Array{Float64,2}
    "Affine transform from x to the barycentric coordinates c (nsimplex, ndim+1,ndim)"
    transform::Array{Float64,3}
    "Lookup array, from a vertex, to some simplex which it is a part of (npoints)"
    vertex_to_simplex::Array{Int,1}
    "Vertices of facets forming the convex hull of the point set (nfaces, ndim)"
    convex_hull::Array{Int,2}
    "Indices of coplanar points and the corresponding indices of the nearest facet and the nearest vertex (ncoplanar, 3)"
    coplanar::Array{Int,2}
    "Neighboring vertices of vertices [neighbors, vertex]"
    vertex_neighbor_vertices::Union{Nothing,SparseMatrixCSC{Nothing,Int}}
end

function Base.show(io::IO, tri::Triangulation)
    println(io,
            "Triangulation of $(size(tri.points, 1)) vertices in $(size(tri.points, 2)) dimensions")
    println(io, "Coordinates of vertices:")
    println(io, "    ", tri.points)
    println(io, "Indices of simplices:")
    println(io, "    ", tri.simplices)
    println(io, "Lookup array from vertices to some containing simplex:")
    println(io, "    ", tri.vertex_to_simplex)
    println(io, "Indices of facets of convex hull of the vertices:")
    println(io, "    ", tri.convex_hull)
    println(io,
            "Indices of coplanar points, nearest facet, and nearest vertex:")
    println(io, "    ", isempty(tri.coplanar) ? "[]" : tri.coplanar)
    println(io, "Neighboring vertices of vertices:")
    vnv = tri.vertex_neighbor_vertices
    if vnv === nothing
        println(io, "    $vnv")
    else
        for j in 1:size(vnv, 2)
            print(io, "    [$j]: [")
            comma = ""
            for i0 in nzrange(vnv, j)
                i = rowvals(vnv)[i0]
                print(io, comma, i)
                comma = ", "
            end
            println(io, "]")
        end
    end
end

"""
    delaunay(vertices::Array{Float64, 2})::Triangulation

Calculate the Delaunay triangulation of the given vertices.

Input format: `size(vertices) = N, D` where `N` is the number of
vertices and `D` the spatial dimension.

Algorithm: Uses SciPy's `Delaunay` function. See
<https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.Delaunay.html>.
"""
function delaunay(vertices::Array{Float64,2},
                  options::Union{Nothing,Vector{Symbol}}=nothing)::Triangulation
    nvertices, dim = size(vertices)
    @assert dim > 0
    dim == 1 && return delaunay_1d(vertices)
    if options === nothing
        if dim <= 4
            options = [:Qbb, :Qc, :Qz, :Q12]
        else
            options = [:Qbb, :Qc, :Qz, :Qx, :Q12]
        end
    end
    options = join(options, " ")
    py = spatial.Delaunay(vertices; qhull_options=options)
    points = convert(Array{Float64,2}, py."points")
    simplices = inc!(convert(Array{Int,2}, py."simplices"))
    neighbors = inc!(convert(Array{Int,2}, py."neighbors"))
    equations = convert(Array{Float64,2}, py."equations")
    transform = convert(Array{Float64,3}, py."transform")
    vertex_to_simplex = inc!(convert(Array{Int,1}, py."vertex_to_simplex"))
    convex_hull = inc!(convert(Array{Int,2}, py."convex_hull"))
    coplanar = inc!(convert(Array{Int,2}, py."coplanar"))
    if hasproperty(py, "vertex_neighbor_vertices")
        indptr = inc!(convert(Array{Int,1},
                              get(py."vertex_neighbor_vertices", 0)))
        indices = inc!(convert(Array{Int,1},
                               get(py."vertex_neighbor_vertices", 1)))
        vertex_neighbor_vertices = SparseMatrixCSC{Nothing,Int}(nvertices,
                                                                nvertices,
                                                                indptr, indices,
                                                                fill(nothing,
                                                                     length(indices)))
    else
        vertex_neighbor_vertices = nothing
    end
    return Triangulation(points, simplices, neighbors, equations, transform,
                         vertex_to_simplex, convex_hull, coplanar,
                         vertex_neighbor_vertices)
end

function delaunay_1d(vertices::Array{Float64,2})::Triangulation
    nvertices, dim = size(vertices)
    @assert dim == 1

    points = vertices

    # Sort vertices
    perm = sortperm(reshape(vertices, nvertices))

    nsimplices = nvertices - 1
    simplices = Array{Int}(undef, nsimplices, 2)
    for i in 1:nsimplices
        pi0 = perm[i]
        pi1 = perm[i + 1]
        simplices[i, 1] = min(pi0, pi1)
        simplices[i, 2] = max(pi0, pi1)
    end

    neighbors = Array{Int}(undef, nsimplices, dim + 1)
    for i in 1:nsimplices
        neighbors[i, 1] = i == 1 ? 0 : i - 1
        neighbors[i, 2] = i == nsimplices ? 0 : i + 1
    end

    equations = Array{Int}(undef, 0, dim + 2)          # TODO
    transform = Array{Float64}(undef, 0, dim + 1, dim) # TODO

    vertex_to_simplex = Array{Int}(undef, nvertices)
    for i in 1:nvertices
        vertex_to_simplex[perm[i]] = i == nvertices ? nvertices - 1 : i
    end

    convex_hull = Array{Int}(undef, 2, dim)
    convex_hull[1] = perm[1]
    convex_hull[2] = perm[nvertices]

    coplanar = Array{Int}(undef, 0, 3)

    # Convert to sparse matrix format
    I = Int[]
    J = Int[]
    V = Nothing[]
    for i in 1:nvertices
        if i > 1
            push!(I, perm[i])
            push!(J, perm[i - 1])
            push!(V, nothing)
        end
        if i < nvertices
            push!(I, perm[i])
            push!(J, perm[i + 1])
            push!(V, nothing)
        end
    end
    # We transpose I and J here
    vertex_neighbor_vertices = sparse(J, I, V, nvertices, nvertices)
    return Triangulation(points, simplices, neighbors, equations, transform,
                         vertex_to_simplex, convex_hull, coplanar,
                         vertex_neighbor_vertices)
end

"""
   inc!(A::AbstractArray)

Convert Python zero-based to Julia 1-based array indices
"""
function inc!(A)
    for i in eachindex(A)
        A[i] += 1
    end
    return A
end

end


text:
#=
Toeplitz matrix examples from Trefethen & Embree, Spectra & Pseudospectra
mainly ch. 7

Instructions
============

1. Establish a plotter for Pseudospectra:
```
using Pseudospectra
setpsplotter(:PyPlot) -or- setpsplotter(:Plots)
```
If using Plots, you should probably set a backend first.

2. `include()` this file.
3. Invoke `M.runme()` with arguments of interest.
=#
module M

using Pseudospectra, ToeplitzMatrices, Polynomials, LinearAlgebra

p = getpsplotter()
if p == :undef
    error("plotter must be set first")
elseif p == :PyPlot
    using PyPlot
else
    using Plots
end

"""
Draws a spectral portrait of a Toeplitz matrix, also showing
the locations of eigenvalues of perturbations of said matrix.
"""
function runme(casename="limaçon",n=64,niter=100)

    plotter = getpsplotter()

    gs = setgs()

    ax = zeros(0) # override as needed

    # conventions: x^0 terms must be equal
    # make life easier for those of us w/o French keyboards
    if casename[1:3] == "lim" # limaçon
        ax = [-2,2,-2,2]
        vc = [0.0,1.0,1.0]
        vr = [0.0,0.0]
    elseif casename[1:3] == "tri" # triangle
        vc = [0.0,1.0]
        vr = [0.0,0.0,1/4]
    elseif casename[1:3] == "wha" # whale
        vr = [0.0,1.0,1im,-3-2im,-1]
        vc = [0.0,10,3+1im,4,1im]
    elseif casename[1:3] == "but" # butterfly
        vc = [0.0,-1im,1.0]
        vr = [0.0,1im,-1]
    elseif casename[1:3] == "grc" # Grcar
        vc = ones(4)
        vr = [1.0,-1.0]
    elseif casename[1:3] == "ch3" # chapter 3 example
        ax = [-1.5,1.5,-1,1]
        vc = [0.0,1/4]
        vr = [0.0,1.0]
    elseif casename[1:3] == "ch7" # chapter 7 example
        vc = [0.0,0,-4,-2im]
        vr = [0.0,2im,-1,2]
    else
        throw(ArgumentError("unrecognized casename"))
    end

    vrp = Poly(vr)
    vcp = Poly(vcat([0.0],vc[2:end]))

    symbolfunc(z) = polyval(vrp,z) + polyval(vcp,1/z)

    At = Toeplitz(vcat(vc,zeros(n-length(vc))),vcat(vr,zeros(n-length(vr))))
    # some methods are missing for Toeplitz, so ...
    A = Matrix(At)

    opts=Dict{Symbol,Any}()
    isempty(ax) || (opts[:ax] = ax)
    ps_data = new_matrix(A,opts)
    driver!(ps_data,opts,gs)

    ϵ = 0.001
    for it=1:niter
        pert = randn(size(A)) + im*randn(size(A))
        pn = norm(pert)
        pert .*= ϵ/pn
        A1 = A + pert
        ews = eigvals(A1)
        if plotter == :PyPlot
            plot(real(ews),imag(ews),"r.")
        else
            scatter!(real(ews),imag(ews),markersize=2,c=:red,leg=false,
                     markerstrokealpha=0.0)
        end
    end
    θv = range(0, stop=2π, length=201)

    symcurve = [symbolfunc(cis(θ)) for θ in θv]
    if plotter == :PyPlot
        plot(real(symcurve),imag(symcurve),"k")
        p = nothing
    else
        p = plot!(real(symcurve),imag(symcurve))
    end

    #=
    # Sometimes perturbed ews are inside this curve.
    # Theory seems confused.

    efac = ϵ^(-1/n)
    symcurve = [symbolfunc(efac*cis(θ)) for θ in θv]
    plot(real(symcurve),imag(symcurve),"g")
    =#

    p
end

end


text:
# Generate a 100 grid of squares

V, EV = Lar.randomcuboids(100, .4);
V = GL.normalize2(V, flag = true);

model = CAGD.Model(V);
CAGD.addModelCells!(model, 1, EV);

# Planar Arrangment

A_model = CAGD.planar_arrangement(model);

A_model.G
Matrix(A_model.T[1])
Matrix(A_model.T[2])

# Display Models

if todisplay
    VV = [[k] for k=1:size(A_model.G, 2)];
    triangulated_faces = CAGD.triangulate2D(A_model);
    FVs = convert(Array{Lar.Cells}, triangulated_faces);

    GL.VIEW(GL.GLExplode(A_model.G, FVs, 1.0, 1.0, 1.0, 99, 1));
    GL.VIEW(GL.GLExplode(A_model.G, FVs, 1.2, 1.2, 1.2, 99, 1));
end;

# Generate 50 S1

store = []
mycircle(r,n) = Lar.circle(r)([n])
for k=1:50
    r = rand()
    n = 8
    while true
        n = abs(rand(Int8)+1)
        if n>2 break end
    end
    scale = Lar.s(0.25,0.25)
    transl = Lar.t(rand(2)...)
    s = Lar.Struct([ transl, scale, mycircle(r,n) ])
    push!(store, Lar.struct2lar(s))
end

V,EV = Lar.struct2lar(Lar.Struct(store))
V = GL.normalize2(V)

if todisplay  GL.VIEW([ GL.GLLines(V,EV) ])  end

model = CAGD.Model(V);
CAGD.addModelCells!(model, 1, EV);

# Planar Arrangment

A_model = CAGD.planar_arrangement(model);

A_model.G
Matrix(A_model.T[1])
Matrix(A_model.T[2])

# Display Models

if todisplay
    VV = [[k] for k=1:size(A_model.G, 2)];
    triangulated_faces = CAGD.triangulate2D(A_model);
    FVs = convert(Array{Lar.Cells}, triangulated_faces);

    GL.VIEW(GL.GLExplode(A_model.G, FVs, 1.0, 1.0, 1.0, 99, 1));
    GL.VIEW(GL.GLExplode(A_model.G, FVs, 1.2, 1.2, 1.2, 99, 1));
end;


text:
module NewPkgEval

using BinaryBuilder
import Pkg.TOML
using Pkg
import Base: UUID
using Dates

downloads_dir(name) = joinpath(@__DIR__, "..", "deps", "downloads", name)
julia_path(ver) = joinpath(@__DIR__, "..", "deps", "julia-$ver")
versions_file() = joinpath(@__DIR__, "..", "deps", "Versions.toml")
registry_path() = joinpath(first(DEPOT_PATH), "registries", "General")

include("build_julia.jl")

"""
    get_registry()

Download the default registry, or if it already exists, update it.
"""
function get_registry()
    Pkg.Types.clone_default_registries()
    Pkg.Types.update_registries(Pkg.Types.Context())
end

"""
    read_versions() -> Dict

Parse the `deps/Versions.toml` file containing version and download information for
various versions of Julia.
"""
function read_versions()
    vers = TOML.parse(read(versions_file(), String))
end

"""
    obtain_julia(the_ver)

Download the specified version of Julia using the information provided in `Versions.toml`.
"""
function obtain_julia(the_ver::VersionNumber)
    vers = read_versions()
    for (ver, data) in vers
        ver = VersionNumber(ver)
        ver == the_ver || continue
        if haskey(data, "url")
            file = get(data, "file", "julia-$ver.tar.gz")
            @assert !isabspath(file)
            Pkg.PlatformEngines.download_verify_unpack(
                data["url"],
                data["sha"],
                julia_path(ver);
                tarball_path=downloads_dir(file),
                force=true
            )
        else
            file = data["file"]
            !isabspath(file) && (file = downloads_dir(file))
            Pkg.PlatformEngines.verify(file, data["sha"])
            isdir(julia_path(ver)) || Pkg.PlatformEngines.unpack(file, julia_path(ver))
        end
        return
    end
    error("Requested Julia version not found")
end

function installed_julia_dir(ver)
     jp = julia_path(ver)
     jp_contents = readdir(jp)
     # Allow the unpacked directory to either be insider another directory (as produced by
     # the buildbots) or directly inside the mapped directory (as produced by the BB script)
     if length(jp_contents) == 1
         jp = joinpath(jp, first(jp_contents))
     end
     jp
end

function with_mounted_shards(f, runner)
    BinaryBuilder.mount_shards(runner; verbose=true)
    try
        f()
    finally
        BinaryBuilder.unmount_shards(runner; verbose=true)
    end
end

"""
    run_sandboxed_julia(args=``; ver::VersionNumber, do_obtain=true, kwargs...)

Run Julia inside of a sandbox, passing the given arguments `args` to it. The keyword
argument `ver` specifies the version of Julia to use, and `do_obtain` dictates whether
the specified version should first be downloaded. If `do_obtain` is `false`, it must
already be installed.
"""
function run_sandboxed_julia(args=``; ver::VersionNumber, do_obtain=true, stdout=stdout, stdin=stdin)
    runner, cmd = runner_sandboxed_julia(args; ver=ver, do_obtain=do_obtain)
    with_mounted_shards(runner) do
        Base.run(pipeline(cmd, stdout=stdout, stderr=stderr))
    end
end

function runner_sandboxed_julia(args=``; ver::VersionNumber, do_obtain=true)
    if do_obtain
        obtain_julia(ver)
    else
        @assert ispath(julia_path(ver))
    end
    ispath(registry_path()) || error("Please run `NewPkgEval.get_registry()` first")
    tmpdir = joinpath(tempdir(), "NewPkgEval")
    mkpath(tmpdir)
    tmpdir = mktempdir(tmpdir)
    runner = BinaryBuilder.UserNSRunner(tmpdir,
        workspaces=[
            installed_julia_dir(ver) => "/maps/julia",
            registry_path() => "/maps/registries/General"
        ])
    cmd = `/maps/julia/bin/julia --color=yes $args`
    cmd = setenv(`$(runner.sandbox_cmd) -- $(cmd)`, runner.env) # extracted from run_interactive in BinaryBuilder
    return runner, cmd
end

log_path(ver) = joinpath(@__DIR__, "..", "logs/logs-$ver")

"""
    run_sandboxed_test(pkg; ver::VersionNumber, do_depwarns=false, log_limit = 5*1024^2, time_limit=60*45, kwargs...)

Run the unit tests for a single package `pkg` inside of a sandbox using the Julia version
`ver`. If `do_depwarns` is `true`, deprecation warnings emitted while running the package's
tests will cause the tests to fail. Test will be forcibly interrupted after `time_limit` seconds or if the log
becomes larger than `log_limit`.

A log for the tests is written to a version-specific directory in the NewPkgEval root
directory.
"""
function run_sandboxed_test(pkg::AbstractString; ver, log_limit = 5*1024^2 #= 5 MB =#,time_limit = 45*60, do_depwarns=false, kwargs...)
    mkpath(log_path(ver))
    arg = """
        using Pkg
        # TODO: Possible to remove?
        #open("/etc/hosts", "w") do f
        #   println(f, "127.0.0.1\tlocalhost")
        #end
        # Map the local registry to the sandbox registry
        mkpath("/root/.julia/registries")
        run(`ln -s /maps/registries/General /root/.julia/registries/General`)
        # Prevent Pkg from updating registy on the Pkg.add
        ENV["CI"] = true
        Pkg.UPDATED_REGISTRY_THIS_SESSION[] = true
        Pkg.add($(repr(pkg)))
        Pkg.test($(repr(pkg)))
    """
    cmd = do_depwarns ? `--depwarn=error` : ``
    cmd = `$cmd -e $arg`
    runner, cmd = runner_sandboxed_julia(cmd; ver=ver, kwargs...)

    log = joinpath(log_path(ver), "$pkg.log")
    open(log, "w") do f
        with_mounted_shards(runner) do
            p = Base.run(pipeline(cmd, stdout=f, stderr=f); wait=false)
            t = Timer(time_limit) do timer
                process_running(p) || return # exit callback
                kill_process(p)
            end
            t2 = @async while true
                process_running(p) || break
                if stat(log).size > log_limit
                    kill_process(p)
                    break
                end
                flush(f)
                sleep(2)
            end
            s = success(p)
            close(t)
            wait(t2)
            return s
        end
    end
end

function kill_process(p)
    if BinaryBuilder.runner_override == "privileged"
        pid = getpid(p)
        Base.run(`sudo kill $pid`)
        sleep(3)
        if process_running(p)
            Base.run(`sudo kill -9 $pid`)
        end
    else
        kill(p)
        sleep(3)
        if process_running(p)
            kill(p, 9)
        end
    end
end

# Skip these packages when testing packages
const skip_list = [
    "AbstractAlgebra", # Hangs forever
    "DiscretePredictors", # Hangs forever
    "LinearLeastSquares", # Hangs forever
    "SLEEF", # Hangs forever
    "OrthogonalPolynomials", # Hangs forever
    "IndexableBitVectors",
    "LatinHypercubeSampling", # Hangs forever
    "DynamicalBilliards", # Hangs forever
    "ChangePrecision", # Hangs forever
    "Rectangle", # Hangs forever
    "Parts", # Hangs forever
    "ZippedArrays", # Hangs forever
    "Chunks", # Hangs forever
    "Electron",
    "DotOverloading",
    "ValuedTuples",
    "HCubature",
    "SequentialMonteCarlo",
    "RequirementVersions",
    "NumberedLines",
    "LazyContext",
    "RecurUnroll", # deleted, hangs
    "TypedBools", # deleted, hangs
    "LazyCall", # deleted, hangs
    "MeshCatMechanisms",
    "SessionHacker",
    "Embeddings",
    "GeoStatsDevTools",
    "DataDeps", # hangs
    "DataDepsGenerators", # hangs
    "MackeyGlass", # deleted, hangs
    "Keys", #deleted, hangs
]

# Blindly assume these packages are okay
const ok_list = [
    "BinDeps", # Not really ok, but packages may list it just as a fallback
    "Homebrew",
    "WinRPM",
    "NamedTuples", # As requested by quinnj
    "Compat",
]

# Stdlibs are assumed to be ok
append!(ok_list, readdir(Sys.STDLIB))

"""
    run(depsgraph, ninstances, version[, result]; do_depwarns=false,
        time_limit=60*45)

Run all tests for all packages in the given package dependency graph using `ninstances`
workers and the specified version of Julia. An existing result `Dict` can be specified,
in which case the function will write to that.

If the keyword argument `do_depwarns` is `true`, deprecation warnings emitted in package
tests will cause the package's tests to fail, i.e. Julia is run with `--depwarn=error`.

Tests will be forcibly interrupted after `time_limit` seconds.
"""
function run(pkgs::Vector, ninstances::Integer, ver::VersionNumber, result = Dict{String, Symbol}();
             do_depwarns=false, time_limit=60*45)
    obtain_julia(ver)
    # In case we need to provide sudo password, do that before starting the actual testing
    run_sandboxed_julia(`-e '1'`; ver=ver)
    get_registry() # make sure local registry is updated

    pkgs = copy(pkgs)
    npkgs = length(pkgs)
    ninstances = min(npkgs, ninstances)
    running = Vector{Union{Nothing, Symbol}}(nothing, ninstances)
    times = DateTime[now() for i = 1:ninstances]
    all_workers = Task[]

    done = false
    did_signal_workers = false
    function stop_work()
        if !done
            done = true
            if !did_signal_workers
                for (i, task) in enumerate(all_workers)
                    task == current_task() && continue
                    Base.istaskdone(task) && continue
                    try; schedule(task, InterruptException(); error=true); catch; end
                    running[i] = nothing
                end
                did_signal_workers = true
            end
        end
    end

    try @sync begin
        # Printer
        @async begin
            try
                io = IOContext(IOBuffer(), :color=>true)
                while (!isempty(pkgs) || !all(==(nothing), running)) && !done
                    o = count(==(:ok),      values(result))
                    f = count(==(:fail),    values(result))
                    s = count(==(:skipped), values(result))
                    print(io, "Success: ")
                    printstyled(io, o; color = :green)
                    print(io, "\tFailed: ")
                    printstyled(io, f; color = Base.error_color())
                    print(io, "\tSkipped: ")
                    printstyled(io, s; color = Base.warn_color())
                    println(io, "\tRemaining: ", npkgs - (o + f + s))
                    for i = 1:ninstances
                        r = running[i]
                        if r === nothing
                            println(io, "Worker $i: -------")
                        else
                            time = Dates.canonicalize(Dates.CompoundPeriod(now() - times[i]))
                            pop!(time.periods) # get rid of milliseconds
                            println(io, "Worker $i: $(r) running for ", time)
                        end
                    end
                    print(String(take!(io.io)))
                    sleep(1)
                    CSI = "\e["
                    print(io, "$(CSI)$(ninstances+1)A$(CSI)1G$(CSI)0J")
                end
                stop_work()
            catch e
                stop_work()
                !isa(e, InterruptException) && rethrow(e)
            end
        end

        # Workers
        for i = 1:ninstances
            push!(all_workers, @async begin
                try
                    while !isempty(pkgs) && !done
                        pkg = pop!(pkgs)
                        if pkg.name in skip_list
                            result[pkg.name] = :skip
                        else
                            running[i] = Symbol(pkg.name)
                            times[i] = now()
                            result[pkg.name] = NewPkgEval.run_sandboxed_test(pkg.name; ver=ver, time_limit=time_limit) ? :ok : :fail
                            running[i] = nothing
                        end
                    end
                catch e
                    stop_work()
                    isa(e, InterruptException) || rethrow(e)
                end
            end)
        end
    end
    catch e
        e isa InterruptException || rethrow(e)
    end
    return result
end

"""
    read_pkgs([pkgs::Vector{String}])

Read packages from the default registry and return them as a vector of tuples containing
the package name, its UUID, and a path to it. If `pkgs` is given, only collect packages
matching the names in `pkgs`
"""
function read_pkgs(pkgs::Union{Nothing, Vector{String}}=nothing)
    pkg_data = []
    for registry in (registry_path(),)
        open(joinpath(registry, "Registry.toml")) do io
            for (_uuid, pkgdata) in Pkg.Types.read_registry(joinpath(registry, "Registry.toml"))["packages"]
                uuid = UUID(_uuid)
                name = pkgdata["name"]
                if pkgs !== nothing
                    idx = findfirst(==(name), pkgs)
                    idx === nothing && continue
                    deleteat!(pkgs, idx)
                end
                path = abspath(registry, pkgdata["path"])
                push!(pkg_data, (name=name, uuid=uuid, path=path))
            end
        end
    end
    if pkgs !== nothing && !isempty(pkgs)
        @warn """did not find the following packages in the registry:\n $("  - " .* join(pkgs, '\n'))"""
    end
    pkg_data
end

end # module