text:
using Compat
using ProtoBuf
import ProtoBuf.meta
import Base: hash, isequal, ==

type __enum_AclEntryProto_AclEntryScopeProto <: ProtoEnum
    ACCESS::Int32
    DEFAULT::Int32
    __enum_AclEntryProto_AclEntryScopeProto() = new(0,1)
end #type __enum_AclEntryProto_AclEntryScopeProto
const AclEntryProto_AclEntryScopeProto = __enum_AclEntryProto_AclEntryScopeProto()

type __enum_AclEntryProto_AclEntryTypeProto <: ProtoEnum
    USER::Int32
    GROUP::Int32
    MASK::Int32
    OTHER::Int32
    __enum_AclEntryProto_AclEntryTypeProto() = new(0,1,2,3)
end #type __enum_AclEntryProto_AclEntryTypeProto
const AclEntryProto_AclEntryTypeProto = __enum_AclEntryProto_AclEntryTypeProto()

type __enum_AclEntryProto_FsActionProto <: ProtoEnum
    NONE::Int32
    EXECUTE::Int32
    WRITE::Int32
    WRITE_EXECUTE::Int32
    READ::Int32
    READ_EXECUTE::Int32
    READ_WRITE::Int32
    PERM_ALL::Int32
    __enum_AclEntryProto_FsActionProto() = new(0,1,2,3,4,5,6,7)
end #type __enum_AclEntryProto_FsActionProto
const AclEntryProto_FsActionProto = __enum_AclEntryProto_FsActionProto()

type AclEntryProto
    _type::Int32
    scope::Int32
    permissions::Int32
    name::AbstractString
    AclEntryProto() = (o=new(); fillunset(o); o)
end #type AclEntryProto
const __req_AclEntryProto = Symbol[:_type,:scope,:permissions]
meta(t::Type{AclEntryProto}) = meta(t, __req_AclEntryProto, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES)
hash(v::AclEntryProto) = ProtoBuf.protohash(v)
isequal(v1::AclEntryProto, v2::AclEntryProto) = ProtoBuf.protoisequal(v1, v2)
==(v1::AclEntryProto, v2::AclEntryProto) = ProtoBuf.protoeq(v1, v2)

type AclStatusProto
    owner::AbstractString
    group::AbstractString
    sticky::Bool
    entries::Array{AclEntryProto,1}
    permission::FsPermissionProto
    AclStatusProto() = (o=new(); fillunset(o); o)
end #type AclStatusProto
const __req_AclStatusProto = Symbol[:owner,:group,:sticky]
meta(t::Type{AclStatusProto}) = meta(t, __req_AclStatusProto, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES)
hash(v::AclStatusProto) = ProtoBuf.protohash(v)
isequal(v1::AclStatusProto, v2::AclStatusProto) = ProtoBuf.protoisequal(v1, v2)
==(v1::AclStatusProto, v2::AclStatusProto) = ProtoBuf.protoeq(v1, v2)

type AclEditLogProto
    src::AbstractString
    entries::Array{AclEntryProto,1}
    AclEditLogProto() = (o=new(); fillunset(o); o)
end #type AclEditLogProto
const __req_AclEditLogProto = Symbol[:src]
meta(t::Type{AclEditLogProto}) = meta(t, __req_AclEditLogProto, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES)
hash(v::AclEditLogProto) = ProtoBuf.protohash(v)
isequal(v1::AclEditLogProto, v2::AclEditLogProto) = ProtoBuf.protoisequal(v1, v2)
==(v1::AclEditLogProto, v2::AclEditLogProto) = ProtoBuf.protoeq(v1, v2)

type ModifyAclEntriesRequestProto
    src::AbstractString
    aclSpec::Array{AclEntryProto,1}
    ModifyAclEntriesRequestProto() = (o=new(); fillunset(o); o)
end #type ModifyAclEntriesRequestProto
const __req_ModifyAclEntriesRequestProto = Symbol[:src]
meta(t::Type{ModifyAclEntriesRequestProto}) = meta(t, __req_ModifyAclEntriesRequestProto, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES)
hash(v::ModifyAclEntriesRequestProto) = ProtoBuf.protohash(v)
isequal(v1::ModifyAclEntriesRequestProto, v2::ModifyAclEntriesRequestProto) = ProtoBuf.protoisequal(v1, v2)
==(v1::ModifyAclEntriesRequestProto, v2::ModifyAclEntriesRequestProto) = ProtoBuf.protoeq(v1, v2)

type ModifyAclEntriesResponseProto
    ModifyAclEntriesResponseProto() = (o=new(); fillunset(o); o)
end #type ModifyAclEntriesResponseProto
hash(v::ModifyAclEntriesResponseProto) = ProtoBuf.protohash(v)
isequal(v1::ModifyAclEntriesResponseProto, v2::ModifyAclEntriesResponseProto) = ProtoBuf.protoisequal(v1, v2)
==(v1::ModifyAclEntriesResponseProto, v2::ModifyAclEntriesResponseProto) = ProtoBuf.protoeq(v1, v2)

type RemoveAclRequestProto
    src::AbstractString
    RemoveAclRequestProto() = (o=new(); fillunset(o); o)
end #type RemoveAclRequestProto
const __req_RemoveAclRequestProto = Symbol[:src]
meta(t::Type{RemoveAclRequestProto}) = meta(t, __req_RemoveAclRequestProto, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES)
hash(v::RemoveAclRequestProto) = ProtoBuf.protohash(v)
isequal(v1::RemoveAclRequestProto, v2::RemoveAclRequestProto) = ProtoBuf.protoisequal(v1, v2)
==(v1::RemoveAclRequestProto, v2::RemoveAclRequestProto) = ProtoBuf.protoeq(v1, v2)

type RemoveAclResponseProto
    RemoveAclResponseProto() = (o=new(); fillunset(o); o)
end #type RemoveAclResponseProto
hash(v::RemoveAclResponseProto) = ProtoBuf.protohash(v)
isequal(v1::RemoveAclResponseProto, v2::RemoveAclResponseProto) = ProtoBuf.protoisequal(v1, v2)
==(v1::RemoveAclResponseProto, v2::RemoveAclResponseProto) = ProtoBuf.protoeq(v1, v2)

type RemoveAclEntriesRequestProto
    src::AbstractString
    aclSpec::Array{AclEntryProto,1}
    RemoveAclEntriesRequestProto() = (o=new(); fillunset(o); o)
end #type RemoveAclEntriesRequestProto
const __req_RemoveAclEntriesRequestProto = Symbol[:src]
meta(t::Type{RemoveAclEntriesRequestProto}) = meta(t, __req_RemoveAclEntriesRequestProto, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES)
hash(v::RemoveAclEntriesRequestProto) = ProtoBuf.protohash(v)
isequal(v1::RemoveAclEntriesRequestProto, v2::RemoveAclEntriesRequestProto) = ProtoBuf.protoisequal(v1, v2)
==(v1::RemoveAclEntriesRequestProto, v2::RemoveAclEntriesRequestProto) = ProtoBuf.protoeq(v1, v2)

type RemoveAclEntriesResponseProto
    RemoveAclEntriesResponseProto() = (o=new(); fillunset(o); o)
end #type RemoveAclEntriesResponseProto
hash(v::RemoveAclEntriesResponseProto) = ProtoBuf.protohash(v)
isequal(v1::RemoveAclEntriesResponseProto, v2::RemoveAclEntriesResponseProto) = ProtoBuf.protoisequal(v1, v2)
==(v1::RemoveAclEntriesResponseProto, v2::RemoveAclEntriesResponseProto) = ProtoBuf.protoeq(v1, v2)

type RemoveDefaultAclRequestProto
    src::AbstractString
    RemoveDefaultAclRequestProto() = (o=new(); fillunset(o); o)
end #type RemoveDefaultAclRequestProto
const __req_RemoveDefaultAclRequestProto = Symbol[:src]
meta(t::Type{RemoveDefaultAclRequestProto}) = meta(t, __req_RemoveDefaultAclRequestProto, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES)
hash(v::RemoveDefaultAclRequestProto) = ProtoBuf.protohash(v)
isequal(v1::RemoveDefaultAclRequestProto, v2::RemoveDefaultAclRequestProto) = ProtoBuf.protoisequal(v1, v2)
==(v1::RemoveDefaultAclRequestProto, v2::RemoveDefaultAclRequestProto) = ProtoBuf.protoeq(v1, v2)

type RemoveDefaultAclResponseProto
    RemoveDefaultAclResponseProto() = (o=new(); fillunset(o); o)
end #type RemoveDefaultAclResponseProto
hash(v::RemoveDefaultAclResponseProto) = ProtoBuf.protohash(v)
isequal(v1::RemoveDefaultAclResponseProto, v2::RemoveDefaultAclResponseProto) = ProtoBuf.protoisequal(v1, v2)
==(v1::RemoveDefaultAclResponseProto, v2::RemoveDefaultAclResponseProto) = ProtoBuf.protoeq(v1, v2)

type SetAclRequestProto
    src::AbstractString
    aclSpec::Array{AclEntryProto,1}
    SetAclRequestProto() = (o=new(); fillunset(o); o)
end #type SetAclRequestProto
const __req_SetAclRequestProto = Symbol[:src]
meta(t::Type{SetAclRequestProto}) = meta(t, __req_SetAclRequestProto, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES)
hash(v::SetAclRequestProto) = ProtoBuf.protohash(v)
isequal(v1::SetAclRequestProto, v2::SetAclRequestProto) = ProtoBuf.protoisequal(v1, v2)
==(v1::SetAclRequestProto, v2::SetAclRequestProto) = ProtoBuf.protoeq(v1, v2)

type SetAclResponseProto
    SetAclResponseProto() = (o=new(); fillunset(o); o)
end #type SetAclResponseProto
hash(v::SetAclResponseProto) = ProtoBuf.protohash(v)
isequal(v1::SetAclResponseProto, v2::SetAclResponseProto) = ProtoBuf.protoisequal(v1, v2)
==(v1::SetAclResponseProto, v2::SetAclResponseProto) = ProtoBuf.protoeq(v1, v2)

type GetAclStatusRequestProto
    src::AbstractString
    GetAclStatusRequestProto() = (o=new(); fillunset(o); o)
end #type GetAclStatusRequestProto
const __req_GetAclStatusRequestProto = Symbol[:src]
meta(t::Type{GetAclStatusRequestProto}) = meta(t, __req_GetAclStatusRequestProto, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES)
hash(v::GetAclStatusRequestProto) = ProtoBuf.protohash(v)
isequal(v1::GetAclStatusRequestProto, v2::GetAclStatusRequestProto) = ProtoBuf.protoisequal(v1, v2)
==(v1::GetAclStatusRequestProto, v2::GetAclStatusRequestProto) = ProtoBuf.protoeq(v1, v2)

type GetAclStatusResponseProto
    result::AclStatusProto
    GetAclStatusResponseProto() = (o=new(); fillunset(o); o)
end #type GetAclStatusResponseProto
const __req_GetAclStatusResponseProto = Symbol[:result]
meta(t::Type{GetAclStatusResponseProto}) = meta(t, __req_GetAclStatusResponseProto, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES)
hash(v::GetAclStatusResponseProto) = ProtoBuf.protohash(v)
isequal(v1::GetAclStatusResponseProto, v2::GetAclStatusResponseProto) = ProtoBuf.protoisequal(v1, v2)
==(v1::GetAclStatusResponseProto, v2::GetAclStatusResponseProto) = ProtoBuf.protoeq(v1, v2)

export AclEntryProto_AclEntryScopeProto, AclEntryProto_AclEntryTypeProto, AclEntryProto_FsActionProto, AclEntryProto, AclStatusProto, AclEditLogProto, ModifyAclEntriesRequestProto, ModifyAclEntriesResponseProto, RemoveAclRequestProto, RemoveAclResponseProto, RemoveAclEntriesRequestProto, RemoveAclEntriesResponseProto, RemoveDefaultAclRequestProto, RemoveDefaultAclResponseProto, SetAclRequestProto, SetAclResponseProto, GetAclStatusRequestProto, GetAclStatusResponseProto


text:
using Test

variable_function_type(::Type{<:MOI.AbstractScalarSet}) = MOI.SingleVariable
variable_function_type(::Type{<:MOI.AbstractVectorSet}) = MOI.VectorOfVariables

"""
    eval_variables(varval::Function, f::AbstractFunction)

Returns the value of function `f` if each variable index `vi` is evaluated as
`varval(vi)`. Note that `varval` should return a number, see
[`substitute_variables`](@ref) for a similar function where `varval` returns a
function.
"""
function eval_variables end
eval_variables(varval::Function, f::SVF) = varval(f.variable)
eval_variables(varval::Function, f::VVF) = varval.(f.variables)
function eval_variables(varval::Function, f::SAF)
    return mapreduce(t->eval_term(varval, t), +, f.terms, init=f.constant)
end
function eval_variables(varval::Function, f::VAF)
    out = copy(f.constants)
    for t in f.terms
        out[t.output_index] += eval_term(varval, t.scalar_term)
    end
    out
end
function eval_variables(varval::Function, f::SQF)
    init = zero(f.constant)
    lin = mapreduce(t->eval_term(varval, t), +, f.affine_terms, init=init)
    quad = mapreduce(t->eval_term(varval, t), +, f.quadratic_terms, init=init)
    return lin + quad + f.constant
end
function eval_variables(varval::Function, f::VQF)
    out = copy(f.constants)
    for t in f.affine_terms
        out[t.output_index] += eval_term(varval, t.scalar_term)
    end
    for t in f.quadratic_terms
        out[t.output_index] += eval_term(varval, t.scalar_term)
    end
    out
end
# Affine term
function eval_term(varval::Function, t::MOI.ScalarAffineTerm)
    return t.coefficient * varval(t.variable_index)
end
# Quadratic term
function eval_term(varval::Function, t::MOI.ScalarQuadraticTerm)
    tval = t.coefficient * varval(t.variable_index_1) * varval(t.variable_index_2)
    t.variable_index_1 == t.variable_index_2 ? tval/2 : tval
end

"""
    map_indices(index_map::Function, x)

Substitute any [`MOI.VariableIndex`](@ref) (resp. [`MOI.ConstraintIndex`](@ref))
in `x` by the [`MOI.VariableIndex`](@ref) (resp. [`MOI.ConstraintIndex`](@ref))
of the same type given by `index_map(x)`.

This function is used by implementations of [`MOI.copy_to`](@ref) on constraint
functions, attribute values and submittable values hence it needs to be
implemented for custom types that are meant to be used as attribute or
submittable value.
"""
function map_indices end

"""
    map_indices(variable_map::AbstractDict{T, T}, x) where {T <: MOI.Index}

Shortcut for `map_indices(vi -> variable_map[vi], x)`.
"""
function map_indices(variable_map::AbstractDict{T, T}, x) where {T <: MOI.Index}
    return map_indices(vi -> variable_map[vi], x)
end

const ObjectWithoutIndex = Union{Nothing, DataType, Number, Enum, AbstractString, MOI.AnyAttribute, MOI.AbstractSet, Function, MOI.ModelLike}
const ObjectOrTupleWithoutIndex = Union{ObjectWithoutIndex, Tuple{Vararg{ObjectWithoutIndex}}}
const ObjectOrTupleOrArrayWithoutIndex = Union{ObjectOrTupleWithoutIndex, AbstractArray{<:ObjectOrTupleWithoutIndex}, AbstractArray{<:AbstractArray{<:ObjectOrTupleWithoutIndex}}}

map_indices(::Function, x::ObjectOrTupleOrArrayWithoutIndex) = x

map_indices(index_map::Function, vi::MOI.VariableIndex) = index_map(vi)
map_indices(index_map::Function, ci::MOI.ConstraintIndex) = index_map(ci)
function map_indices(index_map::Function, array::AbstractArray{<:MOI.Index})
    return map(index_map, array)
end

map_indices(::Function, block::MOI.NLPBlockData) = block

# Terms
function map_indices(index_map::Function, t::MOI.ScalarAffineTerm)
    return MOI.ScalarAffineTerm(t.coefficient, index_map(t.variable_index))
end
function map_indices(index_map::Function, t::MOI.VectorAffineTerm)
    return MOI.VectorAffineTerm(t.output_index, map_indices(index_map, t.scalar_term))
end
function map_indices(index_map::Function, t::MOI.ScalarQuadraticTerm)
    inds = index_map.((t.variable_index_1, t.variable_index_2))
    return MOI.ScalarQuadraticTerm(t.coefficient, inds...)
end
function map_indices(index_map::Function, t::MOI.VectorQuadraticTerm)
    MOI.VectorQuadraticTerm(t.output_index, map_indices(index_map, t.scalar_term))
end

# Functions
function map_indices(index_map::Function, f::MOI.SingleVariable)
    return MOI.SingleVariable(index_map(f.variable))
end
function map_indices(index_map::Function, f::MOI.VectorOfVariables)
    return MOI.VectorOfVariables(index_map.(f.variables))
end
function map_indices(index_map::Function, f::Union{SAF, VAF})
    typeof(f)(map_indices.(index_map, f.terms), MOI.constant(f))
end
function map_indices(index_map::Function, f::Union{SQF, VQF})
    lin = map_indices.(index_map, f.affine_terms)
    quad = map_indices.(index_map, f.quadratic_terms)
    return typeof(f)(lin, quad, MOI.constant(f))
end

# Function changes
map_indices(index_map::Function, change::Union{MOI.ScalarConstantChange, MOI.VectorConstantChange}) = change
function map_indices(index_map::Function, change::MOI.ScalarCoefficientChange)
    return MOI.ScalarCoefficientChange(index_map(change.variable), change.new_coefficient)
end
function map_indices(index_map::Function, change::MOI.MultirowChange)
    return MOI.MultirowChange(index_map(change.variable), change.new_coefficients)
end

# For performance reason, we assume that the type of the function does not
# change in `substitute_variables`.
"""
    substitute_variables(variable_map::Function, x)

Substitute any [`MOI.VariableIndex`](@ref) in `x` by `variable_map(x)`. The
`variable_map` function returns either [`MOI.SingleVariable`](@ref) or
[`MOI.ScalarAffineFunction`](@ref), see [`eval_variables`](@ref) for a similar
function where `variable_map` returns a number.

This function is used by bridge optimizers on constraint functions, attribute
values and submittable values when at least one variable bridge is used hence it
needs to be implemented for custom types that are meant to be used as attribute
or submittable value.
"""
function substitute_variables end

substitute_variables(::Function, x::ObjectOrTupleOrArrayWithoutIndex) = x
substitute_variables(::Function, block::MOI.NLPBlockData) = block

# Used when submitting `HeuristicSolution`.
function substitute_variables(variable_map::Function, vis::Vector{MOI.VariableIndex})
    return substitute_variables.(variable_map, vis)
end
function substitute_variables(variable_map::Function, vi::MOI.VariableIndex)
    func = variable_map(vi)
    if func != MOI.SingleVariable(vi)
        error("Cannot substitute `$vi` as it is bridged into `$func`.")
    end
    return vi
end

function substitute_variables(variable_map::Function,
                              term::MOI.ScalarQuadraticTerm{T}) where T
    f1 = variable_map(term.variable_index_1)
    f2 = variable_map(term.variable_index_2)
    f12 = operate(*, T, f1, f2)::MOI.ScalarQuadraticFunction{T}
    coef = term.coefficient
    # The quadratic terms are evaluated as x'Qx/2 so a diagonal term should
    # be divided by 2 while an off-diagonal term appears twice in the matrix
    # and is divided by 2 so it stays the same.
    if term.variable_index_1 == term.variable_index_2
        coef /= 2
    end
    return operate!(*, T, f12, coef)
end
function substitute_variables(variable_map::Function,
                              term::MOI.ScalarAffineTerm{T}) where T
    return operate(*, T, term.coefficient,
                       variable_map(term.variable_index))::MOI.ScalarAffineFunction{T}
end
function substitute_variables(
    variable_map::Function,
    func::MOI.ScalarAffineFunction{T}) where T
    g = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm{T}[], MOI.constant(func))
    for term in func.terms
        operate!(+, T, g, substitute_variables(variable_map, term))::typeof(func)
    end
    return g
end
function substitute_variables(
    variable_map::Function,
    func::MOI.VectorAffineFunction{T}) where T
    g = MOI.VectorAffineFunction(MOI.VectorAffineTerm{T}[], copy(MOI.constant(func)))
    for term in func.terms
        sub = substitute_variables(variable_map, term.scalar_term)
        operate_output_index!(+, T, term.output_index, g, sub)::typeof(func)
    end
    return g
end
function substitute_variables(
    variable_map::Function,
    func::MOI.ScalarQuadraticFunction{T}) where T
    g = MOI.ScalarQuadraticFunction(
        MOI.ScalarAffineTerm{T}[], MOI.ScalarQuadraticTerm{T}[], MOI.constant(func))
    for term in func.affine_terms
        operate!(+, T, g, substitute_variables(variable_map, term))::typeof(func)
    end
    for term in func.quadratic_terms
        operate!(+, T, g, substitute_variables(variable_map, term))::typeof(func)
    end
    return g
end
function substitute_variables(
    variable_map::Function,
    func::MOI.VectorQuadraticFunction{T}) where T
    g = MOI.VectorQuadraticFunction(
        MOI.VectorAffineTerm{T}[], MOI.VectorQuadraticTerm{T}[], copy(MOI.constant(func)))
    for term in func.affine_terms
        sub = substitute_variables(variable_map, term.scalar_term)
        operate_output_index!(+, T, term.output_index, g, sub)::typeof(func)
    end
    for term in func.quadratic_terms
        sub = substitute_variables(variable_map, term.scalar_term)
        operate_output_index!(+, T, term.output_index, g, sub)::typeof(func)
    end
    return g
end

# Vector of constants
constant_vector(f::Union{SAF, SQF}) = [f.constant]
constant_vector(f::Union{VAF, VQF}) = f.constants

# Implements iterator interface
"""
    scalar_type(F::Type{<:MOI.AbstractVectorFunction})

Type of functions obtained by indexing objects obtained by calling `eachscalar`
on functions of type `F`.
"""
function scalar_type end
scalar_type(::Type{MOI.VectorOfVariables}) = MOI.SingleVariable
scalar_type(::Type{MOI.VectorAffineFunction{T}}) where T = MOI.ScalarAffineFunction{T}
scalar_type(::Type{MOI.VectorQuadraticFunction{T}}) where T = MOI.ScalarQuadraticFunction{T}

struct ScalarFunctionIterator{F<:MOI.AbstractVectorFunction}
    f::F
end
eachscalar(f::MOI.AbstractVectorFunction) = ScalarFunctionIterator(f)
eachscalar(f::AbstractVector) = f

function Base.iterate(it::ScalarFunctionIterator, state = 1)
    if state > length(it)
        return nothing
    else
        return (it[state], state + 1)
    end
end

function Base.length(it::ScalarFunctionIterator{<:MOI.AbstractVectorFunction})
    return MOI.output_dimension(it.f)
end
Base.eltype(it::ScalarFunctionIterator{VVF}) = SVF
Base.eltype(it::ScalarFunctionIterator{VAF{T}}) where T = SAF{T}
Base.eltype(it::ScalarFunctionIterator{VQF{T}}) where T = SQF{T}
Base.lastindex(it::ScalarFunctionIterator) = length(it)

# Define getindex for Vector functions

function Base.getindex(it::ScalarFunctionIterator{MOI.VectorOfVariables},
                       i::Integer)
    return MOI.SingleVariable(it.f.variables[i])
end
# Returns the scalar terms of output_index i
function scalar_terms_at_index(terms::Vector{<:Union{MOI.VectorAffineTerm, MOI.VectorQuadraticTerm}}, i::Int)
    I = findall(t -> t.output_index == i, terms)
    map(i -> terms[i].scalar_term, I)
end
function Base.getindex(it::ScalarFunctionIterator{<:VAF}, i::Integer)
    SAF(scalar_terms_at_index(it.f.terms, i), it.f.constants[i])
end
function Base.getindex(it::ScalarFunctionIterator{<:VQF}, i::Integer)
    lin = scalar_terms_at_index(it.f.affine_terms, i)
    quad = scalar_terms_at_index(it.f.quadratic_terms, i)
    return SQF(lin, quad, it.f.constants[i])
end

function Base.getindex(it::ScalarFunctionIterator{MOI.VectorOfVariables},
                       I::AbstractVector)
    return MOI.VectorOfVariables(it.f.variables[I])
end
function Base.getindex(it::ScalarFunctionIterator{VAF{T}}, I::AbstractVector) where T
    terms = MOI.VectorAffineTerm{T}[]
    constant = Vector{T}(undef, length(I))
    for (i, j) in enumerate(I)
        g = it[j]
        append!(terms, map(t -> MOI.VectorAffineTerm(i, t), g.terms))
        constant[i] = g.constant
    end
    return VAF(terms, constant)
end
function Base.getindex(it::ScalarFunctionIterator{VQF{T}}, I::AbstractVector) where T
    affine_terms = MOI.VectorAffineTerm{T}[]
    quadratic_terms = MOI.VectorQuadraticTerm{T}[]
    constant = Vector{T}(undef, length(I))
    for (i, j) in enumerate(I)
        g = it[j]
        append!(affine_terms, map(t -> MOI.VectorAffineTerm(i, t),
                                  g.affine_terms))
        append!(quadratic_terms, map(t -> MOI.VectorQuadraticTerm(i, t),
                                     g.quadratic_terms))
        constant[i] = g.constant
    end
    return VQF(affine_terms, quadratic_terms, constant)
end

function zero_with_output_dimension(::Type{Vector{T}}, n::Integer) where T
    return zeros(T, n)
end
function zero_with_output_dimension(::Type{MOI.VectorAffineFunction{T}}, n::Integer) where T
    return MOI.VectorAffineFunction{T}(
        MOI.VectorAffineTerm{T}[],
        zeros(T, n))
end
function zero_with_output_dimension(::Type{MOI.VectorQuadraticFunction{T}}, n::Integer) where T
    return MOI.VectorQuadraticFunction{T}(
        MOI.VectorAffineTerm{T}[],
        MOI.VectorQuadraticTerm{T}[],
        zeros(T, n))
end

"""
    unsafe_add(t1::MOI.ScalarAffineTerm, t2::MOI.ScalarAffineTerm)

Sums the coefficients of `t1` and `t2` and returns an output `MOI.ScalarAffineTerm`. It is unsafe because it uses the `variable_index` of `t1` as the `variable_index` of the output without checking that it is equal to that of `t2`.
"""
function unsafe_add(t1::MOI.ScalarAffineTerm, t2::MOI.ScalarAffineTerm)
    return MOI.ScalarAffineTerm(t1.coefficient + t2.coefficient, t1.variable_index)
end

"""
    unsafe_add(t1::MOI.ScalarQuadraticTerm, t2::MOI.ScalarQuadraticTerm)

Sums the coefficients of `t1` and `t2` and returns an output
`MOI.ScalarQuadraticTerm`. It is unsafe because it uses the `variable_index`'s
of `t1` as the `variable_index`'s of the output without checking that they are
the same (up to permutation) to those of `t2`.
"""
function unsafe_add(t1::MOI.ScalarQuadraticTerm, t2::MOI.ScalarQuadraticTerm)
    return MOI.ScalarQuadraticTerm(t1.coefficient + t2.coefficient,
                                   t1.variable_index_1,
                                   t1.variable_index_2)
end

"""
    unsafe_add(t1::MOI.VectorAffineTerm, t2::MOI.VectorAffineTerm)

Sums the coefficients of `t1` and `t2` and returns an output `MOI.VectorAffineTerm`. It is unsafe because it uses the `output_index` and `variable_index` of `t1` as the `output_index` and `variable_index` of the output term without checking that they are equal to those of `t2`.
"""
function unsafe_add(t1::VT, t2::VT) where VT <: Union{MOI.VectorAffineTerm,
                                                      MOI.VectorQuadraticTerm}
    scalar_term = unsafe_add(t1.scalar_term, t2.scalar_term)
    return VT(t1.output_index, scalar_term)
end

is_canonical(::Union{MOI.SingleVariable, MOI.VectorOfVariables}) = true

"""
    is_canonical(f::Union{ScalarAffineFunction, VectorAffineFunction})

Returns a Bool indicating whether the function is in canonical form.
See [`canonical`](@ref).
"""
function is_canonical(f::Union{SAF, VAF})
    is_strictly_sorted(f.terms, MOI.term_indices,
                       t -> !iszero(MOI.coefficient(t)))
end

"""
    is_canonical(f::Union{ScalarQuadraticFunction, VectorQuadraticFunction})

Returns a Bool indicating whether the function is in canonical form.
See [`canonical`](@ref).
"""
function is_canonical(f::Union{SQF, VQF})
    v = is_strictly_sorted(f.affine_terms, MOI.term_indices,
                           t -> !iszero(MOI.coefficient(t)))
    v &= is_strictly_sorted(f.quadratic_terms, MOI.term_indices,
                            t -> !iszero(MOI.coefficient(t)))
end

"""
    is_strictly_sorted(x::AbstractVector, by, filter)

Returns `true` if `by(x[i]) < by(x[i + 1])` and `filter(x[i]) == true` for
all indices i.
"""
function is_strictly_sorted(x::AbstractVector, by, filter)
    if isempty(x)
        return true
    end
    if !filter(first(x))
        return false
    end
    for i in eachindex(x)[2:end]
        if by(x[i]) <= by(x[i - 1])
            return false
        end
        if !filter(x[i])
            return false
        end
    end
    return true
end

"""
    canonical(f::Union{ScalarAffineFunction, VectorAffineFunction,
                       ScalarQuadraticFunction, VectorQuadraticFunction})

Returns the function in a canonical form, i.e.
* A term appear only once.
* The coefficients are nonzero.
* The terms appear in increasing order of variable where there the order of the variables is the order of their value.
* For a `AbstractVectorFunction`, the terms are sorted in ascending order of output index.

The output of `canonical` can be assumed to be a copy of `f`, even for `VectorOfVariables`.

### Examples
If `x` (resp. `y`, `z`) is `VariableIndex(1)` (resp. 2, 3).
The canonical representation of `ScalarAffineFunction([y, x, z, x, z], [2, 1, 3, -2, -3], 5)` is `ScalarAffineFunction([x, y], [-1, 2], 5)`.

"""
canonical(f::MOI.AbstractFunction) = canonicalize!(copy(f))

canonicalize!(f::Union{MOI.VectorOfVariables, MOI.SingleVariable}) = f

"""
    canonicalize!(f::Union{ScalarAffineFunction, VectorAffineFunction})

Convert a function to canonical form in-place, without allocating a copy to hold
the result. See [`canonical`](@ref).
"""
function canonicalize!(f::Union{SAF, VAF})
    sort_and_compress!(f.terms, MOI.term_indices,
                       t -> !iszero(MOI.coefficient(t)), unsafe_add)
    return f
end

"""
    canonicalize!(f::Union{ScalarQuadraticFunction, VectorQuadraticFunction})

Convert a function to canonical form in-place, without allocating a copy to hold
the result. See [`canonical`](@ref).
"""
function canonicalize!(f::Union{SQF, VQF})
    sort_and_compress!(f.affine_terms, MOI.term_indices,
                       t -> !iszero(MOI.coefficient(t)), unsafe_add)
    sort_and_compress!(f.quadratic_terms, MOI.term_indices,
                       t -> !iszero(MOI.coefficient(t)), unsafe_add)
    return f
end

"""
    sort_and_compress!(x::AbstractVector, by::Function, keep::Function, combine::Function)

Sort the vector `x` in-place using `by` as the function from elements to comparable keys, then
combine all entries for which `by(x[i]) == by(x[j])` using the function `x[i] = combine(x[i], x[j])`,
and remove any entries for which `keep(x[i]) == false`. This may result in `x` being resized to
a shorter length.
"""
function sort_and_compress!(x::AbstractVector, by, keep, combine)
    if length(x) > 0
        sort!(x, QuickSort, Base.Order.ord(isless, by, false, Base.Sort.Forward))
        i1 = firstindex(x)
        for i2 in eachindex(x)[2:end]
            if by(x[i1]) == by(x[i2])
                x[i1] = combine(x[i1], x[i2])
            else
                if !keep(x[i1])
                    x[i1] = x[i2]
                else
                    x[i1 + 1] = x[i2]
                    i1 += 1
                end
            end
        end
        if !keep(x[i1])
            i1 -= 1
        end
        resize!(x, i1)
    end
    return x
end

function test_variablenames_equal(model, variablenames)
    seen_name = Dict(name => false for name in variablenames)
    for index in MOI.get(model, MOI.ListOfVariableIndices())
        vname = MOI.get(model, MOI.VariableName(), index)
        if !haskey(seen_name, vname)
            error("Variable with name $vname present in model but not expected list of variable names.")
        end
        if seen_name[vname]
            error("Variable with name $vname present twice in model (shouldn't happen!)")
        end
        seen_name[vname] = true
    end
    for (vname,seen) in seen_name
        if !seen
            error("Did not find variable with name $vname in instance.")
        end
    end
end
function test_constraintnames_equal(model, constraintnames)
    seen_name = Dict(name => false for name in constraintnames)
    for (F,S) in MOI.get(model, MOI.ListOfConstraints())
        for index in MOI.get(model, MOI.ListOfConstraintIndices{F,S}())
            cname = MOI.get(model, MOI.ConstraintName(), index)
            if !haskey(seen_name, cname)
                error("Constraint with name $cname present in model but not expected list of constraint names.")
            end
            if seen_name[cname]
                error("Constraint with name $cname present twice in model (shouldn't happen!)")
            end
            seen_name[cname] = true
        end
    end
    for (cname,seen) in seen_name
        if !seen
            error("Did not find constraint with name $cname in instance.")
        end
    end
end

"""
    all_coefficients(p::Function, f::MOI.AbstractFunction)

Determine whether predicate `p` returns `true` for all coefficients of `f`,
returning `false` as soon as the first coefficient of `f` for which `p`
returns `false` is encountered (short-circuiting). Similar to `all`.
"""
function all_coefficients end

function all_coefficients(p::Function, f::MOI.ScalarAffineFunction)
    return p(f.constant) && all(t -> p(MOI.coefficient(t)), f.terms)
end
function all_coefficients(p::Function, f::MOI.ScalarQuadraticFunction)
    return p(f.constant) &&
        all(t -> p(MOI.coefficient(t)), f.affine_terms) &&
        all(t -> p(MOI.coefficient(t)), f.quadratic_terms)
end

"""
    isapprox_zero(f::MOI.AbstractFunction, tol)

Return a `Bool` indicating whether the function `f` is approximately zero using
`tol` as a tolerance.

## Important note

This function assumes that `f` does not contain any duplicate terms, you might
want to first call [`canonical`](@ref) if that is not guaranteed.
For instance, given
```julia
f = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1, -1], [x, x]), 0)`.
```
then `isapprox_zero(f)` is `false` but `isapprox_zero(MOIU.canonical(f))` is
`true`.
"""
function isapprox_zero end

isapprox_zero(α::AbstractFloat, tol) = -tol <= α <= tol
isapprox_zero(α::Union{Integer, Rational}, tol) = iszero(α)
function isapprox_zero(f::MOI.AbstractFunction, tol)
    return all_coefficients(α -> isapprox_zero(α, tol), f)
end

_is_constant(f::MOI.ScalarAffineFunction) = isempty(f.terms)
_is_constant(f::MOI.ScalarQuadraticFunction) = isempty(f.affine_terms) && isempty(f.quadratic_terms)

Base.iszero(::MOI.SingleVariable) = false
function Base.iszero(f::Union{MOI.ScalarAffineFunction,
                              MOI.ScalarQuadraticFunction})
    return iszero(MOI.constant(f)) && _is_constant(canonical(f))
end

Base.isone(::MOI.SingleVariable) = false
function Base.isone(f::Union{MOI.ScalarAffineFunction,
                             MOI.ScalarQuadraticFunction})
    return isone(MOI.constant(f)) && _is_constant(canonical(f))
end

"""
    test_models_equal(model1::ModelLike, model2::ModelLike, variablenames::Vector{String}, constraintnames::Vector{String})

Test that `model1` and `model2` are identical using `variablenames` as as keys for the variable names and `constraintnames` as keys for the constraint names. Uses `Base.Test` macros.
"""
function test_models_equal(model1::MOI.ModelLike, model2::MOI.ModelLike, variablenames::Vector{String}, constraintnames::Vector{String})
    # TODO: give test-friendly feedback instead of errors?
    test_variablenames_equal(model1, variablenames)
    test_variablenames_equal(model2, variablenames)
    test_constraintnames_equal(model1, constraintnames)
    test_constraintnames_equal(model2, constraintnames)

    variablemap_2to1 = Dict{VI,VI}()
    for vname in variablenames
        index1 = MOI.get(model1, VI, vname)
        index2 = MOI.get(model2, VI, vname)
        variablemap_2to1[index2] = index1
    end

    for cname in constraintnames
        index1 = MOI.get(model1, CI, cname)
        index2 = MOI.get(model2, CI, cname)
        f1 = MOI.get(model1, MOI.ConstraintFunction(), index1)
        f2 = MOI.get(model2, MOI.ConstraintFunction(), index2)
        s1 = MOI.get(model1, MOI.ConstraintSet(), index1)
        s2 = MOI.get(model2, MOI.ConstraintSet(), index2)
        @test isapprox(f1, map_indices(variablemap_2to1, f2))
        @test s1 == s2
    end
    attrs1 = MOI.get(model1, MOI.ListOfModelAttributesSet())
    attrs2 = MOI.get(model2, MOI.ListOfModelAttributesSet())
    attr_list = attrs1 ∪ attrs2
    for attr in attr_list
        value1 = MOI.get(model1, attr)
        value2 = MOI.get(model2, attr)
        if value1 isa MOI.AbstractFunction
            @test value2 isa MOI.AbstractFunction
            @test isapprox(value1, map_indices(variablemap_2to1, value2))
        else
            @test !(value2 isa MOI.AbstractFunction)
            @test value1 == value2
        end
    end
end

_keep_all(keep::Function, v::MOI.VariableIndex) = keep(v)
_keep_all(keep::Function, t::MOI.ScalarAffineTerm) = keep(t.variable_index)
_keep_all(keep::Function, t::MOI.ScalarQuadraticTerm) = keep(t.variable_index_1) && keep(t.variable_index_2)
_keep_all(keep::Function, t::Union{MOI.VectorAffineTerm, MOI.VectorQuadraticTerm}) = _keep_all(keep, t.scalar_term)
# Removes terms or variables in `vis_or_terms` that contains the variable of index `vi`
function _filter_variables(keep::Function, variables_or_terms::Vector)
    return filter(el -> _keep_all(keep, el), variables_or_terms)
end

"""
    filter_variables(keep::Function, f::AbstractFunction)

Return a new function `f` with the variable `vi` such that `!keep(vi)` removed.
"""
function filter_variables end
function filter_variables(keep::Function, f::MOI.SingleVariable)
    if !keep(f.variable)
        error("Cannot remove variable from a `SingleVariable` function of the",
              " same variable.")
    end
    return f
end
function filter_variables(keep::Function, f::MOI.VectorOfVariables)
    return MOI.VectorOfVariables(_filter_variables(keep, f.variables))
end
function filter_variables(
    keep::Function, f::Union{MOI.ScalarAffineFunction, MOI.VectorAffineFunction})
    return typeof(f)(_filter_variables(keep, f.terms), MOI.constant(f))
end
function filter_variables(
    keep, f::Union{MOI.ScalarQuadraticFunction, MOI.VectorQuadraticFunction})
    return typeof(f)(
        _filter_variables(keep, f.affine_terms),
        _filter_variables(keep, f.quadratic_terms),
        MOI.constant(f))
end

"""
    remove_variable(f::AbstractFunction, vi::VariableIndex)

Return a new function `f` with the variable vi removed.
"""
function remove_variable(f::MOI.AbstractFunction, vi::MOI.VariableIndex)
    return filter_variables(v -> v != vi, f)
end
function remove_variable(f::MOI.AbstractFunction, vis::Vector{MOI.VariableIndex})
    # Create a `Set` to test membership in `vis` in O(1).
    set = Set(vis)
    return filter_variables(vi -> !(vi in set), f)
end

"""
    modify_function(f::AbstractFunction, change::AbstractFunctionModification)

Return a new function `f` modified according to `change`.
"""
function modify_function(f::SAF, change::MOI.ScalarConstantChange)
    return SAF(f.terms, change.new_constant)
end
function modify_function(f::VAF, change::MOI.VectorConstantChange)
    return VAF(f.terms, change.new_constant)
end
function modify_function(f::SQF, change::MOI.ScalarConstantChange)
    return SQF(f.affine_terms, f.quadratic_terms, change.new_constant)
end
function modify_function(f::VQF, change::MOI.VectorConstantChange)
    return VQF(f.affine_terms, f.quadratic_terms, change.new_constant)
end
function _modifycoefficient(terms::Vector{<:MOI.ScalarAffineTerm}, variable::VI, new_coefficient)
    terms = copy(terms)
    i = something(findfirst(t -> t.variable_index == variable, terms), 0)
    if iszero(i)
        # The variable was not already in the function
        if !iszero(new_coefficient)
            push!(terms, MOI.ScalarAffineTerm(new_coefficient, variable))
        end
    else
        # The variable was already in the function
        if iszero(new_coefficient)
            deleteat!(terms, i)
        else
            terms[i] = MOI.ScalarAffineTerm(new_coefficient, variable)
        end
    end
    terms
end
function modify_function(f::SAF, change::MOI.ScalarCoefficientChange)
    lin = _modifycoefficient(f.terms, change.variable, change.new_coefficient)
    return SAF(lin, f.constant)
end
function modify_function(f::SQF, change::MOI.ScalarCoefficientChange)
    lin = _modifycoefficient(f.affine_terms, change.variable, change.new_coefficient)
    return SQF(lin, f.quadratic_terms, f.constant)
end
function _modifycoefficients(n, terms::Vector{<:MOI.VectorAffineTerm}, variable::VI, new_coefficients)
    terms = copy(terms)
    # Maps between rows in the `MOI.VectorAffineTerm`s and indices in new_coefficients
    rowmap = Dict(c[1]=>i for (i,c) in enumerate(new_coefficients))
    del = Int[]
    for i in 1:length(terms)
        if terms[i].scalar_term.variable_index == variable
            row = terms[i].output_index
            j = Base.get(rowmap, row, 0)
            if !iszero(j) # If it is zero, it means that the row should not be changed
                if iszero(new_coefficients[j][2])
                    push!(del, i)
                else
                    terms[i] = MOI.VectorAffineTerm(row, MOI.ScalarAffineTerm(new_coefficients[j][2], variable))
                end
                rowmap[row] = 0 # We only change the first term of a row
            end
        end
    end
    deleteat!(terms, del)
    for (row, j) in rowmap
        if !iszero(j)
            push!(terms, MOI.VectorAffineTerm(row, MOI.ScalarAffineTerm(new_coefficients[j][2], variable)))
        end
    end
    terms
end
function modify_function(f::VAF, change::MOI.MultirowChange)
    dim = MOI.output_dimension(f)
    coefficients = change.new_coefficients
    lin = _modifycoefficients(dim, f.terms, change.variable, coefficients)
    VAF(lin, f.constants)
end
function modify_function(f::VQF, change::MOI.MultirowChange)
    dim = MOI.output_dimension(f)
    coefficients = change.new_coefficients
    lin = _modifycoefficients(dim, f.affine_terms, change.variable, coefficients)
    return VQF(lin, f.quadratic_terms, f.constants)
end

# Arithmetic

"""
    operate(op::Function, ::Type{T},
            args::Union{T, MOI.AbstractFunction}...)::MOI.AbstractFunction where T

Returns an `MOI.AbstractFunction` representing the function resulting from the
operation `op(args...)` on functions of coefficient type `T`. No argument can be
modified.
"""
function operate end

# Without `<:Number`, Julia v1.1.1 fails at precompilation with a StackOverflowError.
operate(op::Function, ::Type{T}, α::Union{T, AbstractVector{T}}...) where {T<:Number} = op(α...)

"""
    operate!(op::Function, ::Type{T},
             args::Union{T, MOI.AbstractFunction}...)::MOI.AbstractFunction where T

Returns an `MOI.AbstractFunction` representing the function resulting from the
operation `op(args...)` on functions of coefficient type `T`. The first argument
can be modified. The return type is the same than the method
`operate(op, T, args...)` without `!`.
"""
function operate! end

operate!(op::Function, ::Type{T}, α::Union{T, AbstractVector{T}}...) where {T} = op(α...)

"""
    operate_output_index!(
        op::Function, ::Type{T}, output_index::Integer,
        func::MOI.AbstractVectorFunction
        args::Union{T, MOI.AbstractScalarFunction}...)::MOI.AbstractFunction where T

Returns an `MOI.AbstractVectorFunction` where the function at `output_index`
is the result of the operation `op` applied to the function at `output_index`
of `func` and `args`. The functions at output index different to `output_index`
are the same as the functions at the same output index in `func`. The first
argument can be modified.
"""
function operate_output_index! end

function operate_output_index!(op::Function, ::Type{T}, i::Integer, x::Vector{T}, args...) where T
    x[i] = operate!(op, T, x[i], args...)
end

"""
    promote_operation(op::Function, ::Type{T},
                      ArgsTypes::Type{<:Union{T, MOI.AbstractFunction}}...) where T

Returns the type of the `MOI.AbstractFunction` returned to the call
`operate(op, T, args...)` where the types of the arguments `args` are
`ArgsTypes`.
"""
function promote_operation end

# Helpers

function operate_term(::typeof(+), term::Union{
    MOI.ScalarAffineTerm, MOI.ScalarQuadraticTerm,
    MOI.VectorAffineTerm, MOI.VectorQuadraticTerm})
    return term
end
function operate_term(::typeof(-), term::MOI.ScalarAffineTerm)
    return MOI.ScalarAffineTerm(-term.coefficient, term.variable_index)
end
function operate_term(::typeof(-), term::MOI.ScalarQuadraticTerm)
    return MOI.ScalarQuadraticTerm(-term.coefficient,
                                   term.variable_index_1,
                                   term.variable_index_2)
end
function operate_term(::typeof(-), term::MOI.VectorAffineTerm)
    return MOI.VectorAffineTerm(term.output_index, operate_term(-, term.scalar_term))
end
function operate_term(::typeof(-), term::MOI.VectorQuadraticTerm)
    return MOI.VectorQuadraticTerm(term.output_index, operate_term(-, term.scalar_term))
end

function operate_term(::typeof(*), α::T, t::MOI.ScalarAffineTerm{T}) where T
    MOI.ScalarAffineTerm(α * t.coefficient, t.variable_index)
end
# `<:Number` is a workaround for https://github.com/JuliaOpt/MathOptInterface.jl/issues/980
function operate_term(::typeof(*), t::MOI.ScalarAffineTerm{T}, β::T) where T<:Number
    MOI.ScalarAffineTerm(t.coefficient * β, t.variable_index)
end
function operate_term(::typeof(*), α::T, t::MOI.ScalarAffineTerm{T}, β::T) where T
    MOI.ScalarAffineTerm(α * t.coefficient * β, t.variable_index)
end

function operate_term(::typeof(*), α::T, t::MOI.ScalarQuadraticTerm{T}) where T
    MOI.ScalarQuadraticTerm(α * t.coefficient, t.variable_index_1,
                            t.variable_index_2)
end
function operate_term(::typeof(*), t::MOI.ScalarQuadraticTerm{T}, β::T) where T
    MOI.ScalarQuadraticTerm(t.coefficient * β, t.variable_index_1,
                            t.variable_index_2)
end
function operate_term(::typeof(*), α::T, t::MOI.ScalarQuadraticTerm{T}, β::T) where T
    MOI.ScalarQuadraticTerm(α * t.coefficient * β, t.variable_index_1,
                            t.variable_index_2)
end

function operate_term(::typeof(*), t1::MOI.ScalarAffineTerm,
                      t2::MOI.ScalarAffineTerm)
    coef = t1.coefficient * t2.coefficient
    if t1.variable_index == t2.variable_index
        coef *= 2
    end
    MOI.ScalarQuadraticTerm(coef, t1.variable_index, t2.variable_index)
end

function operate_term(::typeof(*), α::T, t::MOI.VectorAffineTerm{T}) where T
    MOI.VectorAffineTerm(t.output_index, operate_term(*, α, t.scalar_term))
end
function operate_term(::typeof(*), t::MOI.VectorAffineTerm{T}, α::T) where T
    MOI.VectorAffineTerm(t.output_index, operate_term(*, t.scalar_term, α))
end
function operate_term(::typeof(*), α::T, t::MOI.VectorQuadraticTerm{T}) where T
    MOI.VectorQuadraticTerm(t.output_index, operate_term(*, α, t.scalar_term))
end
function operate_term(::typeof(*), t::MOI.VectorQuadraticTerm{T}, α::T) where T
    MOI.VectorQuadraticTerm(t.output_index, operate_term(*, t.scalar_term, α))
end
function operate_term(::typeof(*), t1::MOI.VectorAffineTerm,
                      t2::MOI.VectorAffineTerm)
    @assert t1.output_index == t2.output_index
    MOI.VectorQuadraticTerm(t1.output_index, operate_term(*, t1.scalar_term, t2.scalar_term))
end

function operate_term(::typeof(/), t::MOI.ScalarAffineTerm{T}, α::T) where T
    MOI.ScalarAffineTerm(t.coefficient / α, t.variable_index)
end
function operate_term(::typeof(/), t::MOI.ScalarQuadraticTerm{T}, α::T) where T
    MOI.ScalarQuadraticTerm(t.coefficient / α, t.variable_index_1,
                            t.variable_index_2)
end
function operate_term(::typeof(/), t::MOI.VectorAffineTerm{T}, α::T) where T
    MOI.VectorAffineTerm(t.output_index, operate_term(/, t.scalar_term, α))
end
function operate_term(::typeof(/), t::MOI.VectorQuadraticTerm{T}, α::T) where T
    MOI.VectorQuadraticTerm(t.output_index, operate_term(/, t.scalar_term, α))
end

# Avoid a copy in the case of +
function operate_terms(
    ::typeof(+),
    terms::Vector{<:Union{MOI.ScalarAffineTerm,
                          MOI.ScalarQuadraticTerm,
                          MOI.VectorAffineTerm,
                          MOI.VectorQuadraticTerm}})
    return terms
end
function operate_terms!(::typeof(-),
                        terms::Vector{<:Union{MOI.ScalarAffineTerm,
                                              MOI.ScalarQuadraticTerm}})
    return map!(term -> operate_term(-, term), terms, terms)
end
function operate_terms(
    ::typeof(-),
    terms::Vector{<:Union{MOI.ScalarAffineTerm,
                          MOI.ScalarQuadraticTerm,
                          MOI.VectorAffineTerm,
                          MOI.VectorQuadraticTerm}})
    return map(term -> operate_term(-, term), terms)
end
function operate_terms(::typeof(-),
        terms::Vector{<:Union{MOI.VectorAffineTerm,
                              MOI.VectorQuadraticTerm}})
    return map(term -> operate_term(-, term), terms)
end

function map_terms!(op, func::Union{MOI.ScalarAffineFunction, MOI.VectorAffineFunction})
    map!(op, func.terms, func.terms)
end
function map_terms!(op, func::Union{MOI.ScalarQuadraticFunction, MOI.VectorQuadraticFunction})
    map!(op, func.affine_terms, func.affine_terms)
    map!(op, func.quadratic_terms, func.quadratic_terms)
end

# Functions convertible to a ScalarAffineFunction
const ScalarAffineLike{T} = Union{T, MOI.SingleVariable, MOI.ScalarAffineFunction{T}}
# Functions convertible to a ScalarQuadraticFunction
const ScalarQuadraticLike{T} = Union{ScalarAffineLike{T}, MOI.ScalarQuadraticFunction{T}}

# `ScalarLike` for which `T` is defined to avoid defining, e.g.,
# `+(::SingleVariable, ::Any)` which should rather be
# `+(::SingleVariable, ::Number)`.
const TypedScalarLike{T} = Union{MOI.ScalarAffineFunction{T},
                                 MOI.ScalarQuadraticFunction{T}}
# Used for overloading Base operator functions so `T` is not in the union to
# avoid overloading e.g. `+(::Float64, ::Float64)`
const ScalarLike{T} = Union{MOI.SingleVariable, TypedScalarLike{T}}

# Functions convertible to a VectorAffineFunction
const VectorAffineLike{T} = Union{Vector{T}, MOI.VectorOfVariables, MOI.VectorAffineFunction{T}}
# Functions convertible to a VectorQuadraticFunction
const VectorQuadraticLike{T} = Union{VectorAffineLike{T}, MOI.VectorQuadraticFunction{T}}

# `VectorLike` for which `T` is defined to avoid defining, e.g.,
# `+(::VectorOfVariables, ::Any)` which should rather be
# `+(::VectorOfVariables, ::Number)`.
const TypedVectorLike{T} = Union{MOI.VectorAffineFunction{T},
                                 MOI.VectorQuadraticFunction{T}}
# Used for overloading Base operator functions so `T` is not in the union to
# avoid overloading e.g. `+(::Float64, ::Float64)`
const VectorLike{T} = Union{MOI.VectorOfVariables, TypedVectorLike{T}}

const TypedLike{T} = Union{TypedScalarLike{T}, TypedVectorLike{T}}

###################################### +/- #####################################
## promote_operation

function promote_operation(::typeof(-), ::Type{T},
                           ::Type{<:ScalarAffineLike{T}}) where T
    return MOI.ScalarAffineFunction{T}
end
function promote_operation(::Union{typeof(+), typeof(-)}, ::Type{T},
                           ::Type{<:ScalarAffineLike{T}},
                           ::Type{<:ScalarAffineLike{T}}) where T
    return MOI.ScalarAffineFunction{T}
end
function promote_operation(::typeof(-), ::Type{T},
                           ::Type{<:ScalarQuadraticLike{T}}) where T
    return MOI.ScalarQuadraticFunction{T}
end
function promote_operation(::Union{typeof(+), typeof(-)}, ::Type{T},
                           ::Type{<:ScalarQuadraticLike{T}},
                           ::Type{<:ScalarQuadraticLike{T}}) where T
    return MOI.ScalarQuadraticFunction{T}
end

## operate!
# + with at least 3 arguments
function operate!(op::typeof(+), ::Type{T}, f, g, h, args...) where T
    return operate!(+, T, operate!(op, T, f, g), h, args...)
end

# Unary -
function operate!(op::typeof(-), ::Type{T},
                  f::Union{MOI.ScalarAffineFunction{T},
                           MOI.ScalarQuadraticFunction{T}}) where T
    return MA.mutable_operate!(-, f)
end

# Scalar Affine +/-! ...
function operate!(op::Union{typeof(+), typeof(-)}, ::Type{T},
                  f::MOI.ScalarAffineFunction{T},
                  g::Union{T, MOI.SingleVariable, MOI.ScalarAffineFunction{T}}) where T
    return MA.mutable_operate!(op, f, g)
end
function operate!(op::Union{typeof(+), typeof(-)}, ::Type{T},
                  f::MOI.SingleVariable,
                  g::ScalarQuadraticLike) where T
    return operate(op, T, f, g)
end
function operate!(op::Union{typeof(+), typeof(-)}, ::Type{T},
                  f::MOI.ScalarAffineFunction{T},
                  g::MOI.ScalarQuadraticFunction{T}) where T
    return operate(op, T, f, g)
end
# Scalar Quadratic +/-! ...
function operate!(op::Union{typeof(+), typeof(-)}, ::Type{T},
                  f::MOI.ScalarQuadraticFunction{T},
                  g::Union{T, MOI.SingleVariable, MOI.ScalarAffineFunction{T}, MOI.ScalarQuadraticFunction{T}}) where T
    return MA.mutable_operate!(op, f, g)
end

## operate
# + with at least 3 arguments, can use in-place as the user cannot use
# intermediate results
function operate(op::typeof(+), ::Type{T}, f, g, h, args...) where T
    return operate!(+, T, operate(+, T, f, g), h, args...)
end

# Unary +
function operate(op::typeof(+), ::Type{T}, f::MOI.AbstractFunction) where T
    return f
end

# Scalar number +/- ...
function operate(op::typeof(+), ::Type{T}, α::T, f::ScalarLike{T}) where T
    return operate(op, T, f, α)
end
function operate(op::typeof(-), ::Type{T}, α::T, f::ScalarLike{T}) where T
    return operate!(+, T, operate(-, T, f), α)
end

# Scalar Variable +/- ...
function operate(::typeof(-), ::Type{T}, f::MOI.SingleVariable) where T
    return MOI.ScalarAffineFunction{T}(
        [MOI.ScalarAffineTerm(-one(T), f.variable)], zero(T))
end
function operate(op::Union{typeof(+), typeof(-)}, ::Type{T},
                 f::MOI.SingleVariable, α::T) where T
    return MOI.ScalarAffineFunction{T}(
        [MOI.ScalarAffineTerm(one(T), f.variable)], op(α))
end
function operate(op::Union{typeof(+), typeof(-)}, ::Type{T},
                 f::MOI.SingleVariable,
                 g::MOI.SingleVariable) where T
    return MOI.ScalarAffineFunction{T}(
        [MOI.ScalarAffineTerm(one(T), f.variable),
         MOI.ScalarAffineTerm(op(one(T)), g.variable)],
        zero(T))
end
function operate(op::typeof(+), ::Type{T},
                 f::MOI.SingleVariable,
                 g::Union{MOI.ScalarAffineFunction{T},
                          MOI.ScalarQuadraticFunction{T}}) where T
    return operate(op, T, g, f)
end
function operate(op::typeof(-), ::Type{T},
                 f::MOI.SingleVariable,
                 g::Union{MOI.ScalarAffineFunction{T},
                          MOI.ScalarQuadraticFunction{T}}) where T
    return operate!(+, T, operate(-, T, g), f)
end
# Scalar Affine +/- ...
function operate(op::Union{typeof(-)}, ::Type{T},
                 f::MOI.ScalarAffineFunction{T}) where T
    return MOI.ScalarAffineFunction(operate_terms(op, f.terms),
                                    op(f.constant))
end
function operate(op::Union{typeof(+), typeof(-)}, ::Type{T},
                 f::MOI.ScalarAffineFunction{T},
                 g::ScalarAffineLike{T}) where T
    return operate!(op, T, copy(f), g)
end
function operate(op::Union{typeof(+), typeof(-)}, ::Type{T},
                 f::MOI.ScalarAffineFunction{T},
                 g::MOI.ScalarQuadraticFunction{T}) where T
    MOI.ScalarQuadraticFunction([f.terms; operate_terms(op, g.affine_terms)],
                                operate_terms(op, g.quadratic_terms),
                                op(f.constant, g.constant))
end
# Scalar Quadratic +/- ...
function operate(op::Union{typeof(-)}, ::Type{T},
                 f::MOI.ScalarQuadraticFunction{T}) where T
    return MOI.ScalarQuadraticFunction(
        operate_terms(op, f.affine_terms),
        operate_terms(op, f.quadratic_terms),
        op(f.constant))
end
function operate(op::Union{typeof(+), typeof(-)}, ::Type{T},
                 f::MOI.ScalarQuadraticFunction{T},
                 g::ScalarQuadraticLike{T}) where T
    operate!(op, T, copy(f), g)
end

# To avoid type piracy, we add at least one `ScalarLike` outside of the `...`.
function Base.:+(arg::ScalarLike{T}, args::ScalarLike{T}...) where T
    return operate(+, T, arg, args...)
end
function Base.:+(α::T, arg::TypedScalarLike{T}, args::ScalarLike{T}...) where T
    return operate(+, T, α, arg, args...)
end
function Base.:+(α::Number, f::MOI.SingleVariable)
    return operate(+, typeof(α), α, f)
end
function Base.:+(f::TypedScalarLike{T}, α::T) where T
    return operate(+, T, f, α)
end
function Base.:+(f::MOI.SingleVariable, α::Number)
    return operate(+, typeof(α), f, α)
end
function Base.:-(arg::ScalarLike{T}, args::ScalarLike{T}...) where T
    return operate(-, T, arg, args...)
end
function Base.:-(f::TypedScalarLike{T}, α::T) where T
    return operate(-, T, f, α)
end
function Base.:-(f::MOI.SingleVariable, α::Number)
    return operate(-, typeof(α), f, α)
end
function Base.:-(α::T, f::TypedScalarLike{T}) where T
    return operate(-, T, α, f)
end
function Base.:-(α::Number, f::MOI.SingleVariable)
    return operate(-, typeof(α), α, f)
end

# Vector +/-
###############################################################################
function promote_operation(::typeof(-), ::Type{T},
    ::Type{<:VectorAffineLike{T}}) where T
    return MOI.VectorAffineFunction{T}
end
function promote_operation(::typeof(-), ::Type{T},
    ::Type{<:VectorQuadraticLike{T}}) where T
    return MOI.VectorQuadraticFunction{T}
end
function promote_operation(::Union{typeof(+), typeof(-)}, ::Type{T},
                           ::Type{<:VectorAffineLike{T}},
                           ::Type{<:VectorAffineLike{T}}) where T
    return MOI.VectorAffineFunction{T}
end
function promote_operation(::Union{typeof(+), typeof(-)}, ::Type{T},
                           ::Type{<:VectorQuadraticLike{T}},
                           ::Type{<:VectorQuadraticLike{T}}) where T
    return MOI.VectorQuadraticFunction{T}
end

# Vector Variable +/- ...
function operate!(op::Union{typeof(+), typeof(-)}, ::Type{T},
                  f::MOI.VectorOfVariables,
                  g::VectorQuadraticLike) where T
    return operate(op, T, f, g)
end
# Vector Affine +/-! ...
function operate_output_index!(
    op::Union{typeof(+), typeof(-)}, ::Type{T},
    output_index::Integer,
    f::MOI.VectorAffineFunction{T}, α::T) where T
    f.constants[output_index] = op(f.constants[output_index], α)
    return f
end
function operate!(op::Union{typeof(+), typeof(-)}, ::Type{T},
                  f::MOI.VectorAffineFunction{T},
                  g::Vector{T}) where T
    @assert MOI.output_dimension(f) == length(g)
    f.constants .= op.(f.constants, g)
    return f
end
function operate_output_index!(
    op::Union{typeof(+), typeof(-)}, ::Type{T},
    output_index::Integer,
    f::MOI.VectorAffineFunction{T},
    g::MOI.SingleVariable) where T
    push!(f.terms, MOI.VectorAffineTerm(
        output_index, MOI.ScalarAffineTerm(op(one(T)), g.variable)))
    return f
end
function operate!(op::Union{typeof(+), typeof(-)}, ::Type{T},
                  f::MOI.VectorAffineFunction{T},
                  g::MOI.VectorOfVariables) where T
    d = MOI.output_dimension(g)
    @assert MOI.output_dimension(f) == d
    append!(f.terms, MOI.VectorAffineTerm.(
                         collect(1:d),
                         MOI.ScalarAffineTerm.(op(one(T)), g.variables)))
    return f
end
function operate_output_index!(
    op::Union{typeof(+), typeof(-)}, ::Type{T},
    output_index::Integer,
    f::MOI.VectorAffineFunction{T},
    g::MOI.ScalarAffineFunction{T}) where T
    append!(f.terms, MOI.VectorAffineTerm.(
        output_index, operate_terms(op, g.terms)))
    return operate_output_index!(op, T, output_index, f, MOI.constant(g))
end
function operate!(op::Union{typeof(+), typeof(-)}, ::Type{T},
                  f::MOI.VectorAffineFunction{T},
                  g::MOI.VectorAffineFunction{T}) where T
    append!(f.terms, operate_terms(op, g.terms))
    f.constants .= op.(f.constants, g.constants)
    return f
end
function operate!(op::Union{typeof(+), typeof(-)}, ::Type{T},
                  f::MOI.VectorAffineFunction{T},
                  g::MOI.VectorQuadraticFunction{T}) where T
    return operate(op, T, f, g)
end
# Vector Quadratic +/-! ...
function operate_output_index!(
    op::Union{typeof(+), typeof(-)}, ::Type{T},
    output_index::Integer,
    f::MOI.VectorQuadraticFunction{T}, α::T) where T
    f.constants[output_index] = op(f.constants[output_index], α)
    return f
end
function operate!(op::Union{typeof(+), typeof(-)}, ::Type{T},
                  f::MOI.VectorQuadraticFunction{T},
                  g::Vector{T}) where T
    @assert MOI.output_dimension(f) == length(g)
    f.constants .= op.(f.constants, g)
    return f
end
function operate!(op::Union{typeof(+), typeof(-)}, ::Type{T},
                  f::MOI.VectorQuadraticFunction{T},
                  g::MOI.VectorOfVariables) where T
    d = MOI.output_dimension(g)
    @assert MOI.output_dimension(f) == d
    append!(f.affine_terms, MOI.VectorAffineTerm.(collect(1:d), MOI.ScalarAffineTerm.(op(one(T)), g.variables)))
    return f
end
function operate_output_index!(
    op::Union{typeof(+), typeof(-)}, ::Type{T},
    output_index::Integer,
    f::MOI.VectorQuadraticFunction{T},
    g::MOI.ScalarAffineFunction{T}) where T
    append!(f.affine_terms, MOI.VectorAffineTerm.(
        output_index, operate_terms(op, g.terms)))
    return operate_output_index!(op, T, output_index, f, MOI.constant(g))
end
function operate!(op::Union{typeof(+), typeof(-)}, ::Type{T},
                  f::MOI.VectorQuadraticFunction{T},
                  g::MOI.VectorAffineFunction{T}) where T
    append!(f.affine_terms, operate_terms(op, g.terms))
    f.constants .= op.(f.constants, g.constants)
    return f
end
function operate_output_index!(
    op::Union{typeof(+), typeof(-)}, ::Type{T},
    output_index::Integer,
    f::MOI.VectorQuadraticFunction{T},
    g::MOI.ScalarQuadraticFunction{T}) where T
    append!(f.affine_terms, MOI.VectorAffineTerm.(
        output_index, operate_terms(op, g.affine_terms)))
    append!(f.quadratic_terms, MOI.VectorQuadraticTerm.(
        output_index, operate_terms(op, g.quadratic_terms)))
    return operate_output_index!(op, T, output_index, f, MOI.constant(g))
end
function operate!(op::Union{typeof(+), typeof(-)}, ::Type{T},
                  f::MOI.VectorQuadraticFunction{T},
                  g::MOI.VectorQuadraticFunction{T}) where T
    append!(f.affine_terms, operate_terms(op, g.affine_terms))
    append!(f.quadratic_terms, operate_terms(op, g.quadratic_terms))
    f.constants .= op.(f.constants, g.constants)
    return f
end

## operate
# + with at least 3 arguments, can use in-place as the user cannot use
# intermediate results
# overload
# function operate(op::typeof(+), ::Type{T}, f, g, h, args...) where T
#     return operate!(+, T, operate(+, T, f, g), h, args...)
# end

# function operate(op::typeof(+), ::Type{T}, f::VectorOfVariables) where T
#     return f
# end

function operate(op::typeof(-), ::Type{T}, f::MOI.VectorOfVariables) where T
    d = MOI.output_dimension(f)
    return MOI.VectorAffineFunction{T}(
               MOI.VectorAffineTerm.(
                   collect(1:d),
                   MOI.ScalarAffineTerm.(-one(T), f.variables)),
               fill(zero(T),d))
end

# Vector number +/- ...
function operate(op::typeof(+), ::Type{T}, α::Vector{T}, f::VectorLike{T}) where T
    return operate(op, T, f, α)
end
function operate(op::typeof(-), ::Type{T}, α::Vector{T}, f::VectorLike{T}) where T
    return operate!(+, T, operate(-, T, f), α)
end

# Vector Variable +/- ...
function operate(op::Union{typeof(+), typeof(-)}, ::Type{T},
    f::MOI.VectorOfVariables, α::Vector{T}) where T
    d = MOI.output_dimension(f)
    @assert length(α) == d
    return MOI.VectorAffineFunction{T}(
               MOI.VectorAffineTerm.(
                   collect(1:d),
                   MOI.ScalarAffineTerm.(one(T), f.variables)),
               op.(α))
end
function operate(op::Union{typeof(+), typeof(-)}, ::Type{T},
                 f::MOI.VectorOfVariables,
                 g::MOI.VectorOfVariables) where T
    d = MOI.output_dimension(f)
    @assert MOI.output_dimension(g) == d
    return MOI.VectorAffineFunction{T}(
               vcat(
                   MOI.VectorAffineTerm.(
                       collect(1:d),
                       MOI.ScalarAffineTerm.(one(T), f.variables)),
                   MOI.VectorAffineTerm.(
                       collect(1:d),
                       MOI.ScalarAffineTerm.(op(one(T)), g.variables))),
               fill(zero(T),d))
end
function operate(op::typeof(+), ::Type{T},
    f::MOI.VectorOfVariables,
    g::Union{MOI.VectorAffineFunction{T},
                MOI.VectorQuadraticFunction{T}}) where T
    return operate(op, T, g, f)
end
function operate(op::typeof(-), ::Type{T},
    f::MOI.VectorOfVariables,
    g::Union{MOI.VectorAffineFunction{T},
                MOI.VectorQuadraticFunction{T}}) where T
    return operate!(+, T, operate(-, T, g), f)
end
# Vector Affine +/- ...
function operate(op::Union{typeof(-)}, ::Type{T},
    f::MOI.VectorAffineFunction{T}) where T
    return MOI.VectorAffineFunction(operate_terms(op, f.terms),
                      op.(f.constants))
end
function operate(op::Union{typeof(-)}, ::Type{T},
    f::MOI.VectorQuadraticFunction{T}) where T
    return MOI.VectorQuadraticFunction(
        operate_terms(op, f.affine_terms),
        operate_terms(op, f.quadratic_terms),
                      op.(f.constants))
end
function operate(op::Union{typeof(+), typeof(-)}, ::Type{T},
    f::MOI.VectorAffineFunction{T},
    g::VectorAffineLike{T}) where T
    return operate!(op, T, copy(f), g)
end
function operate(op::Union{typeof(+), typeof(-)}, ::Type{T},
    f::MOI.VectorAffineFunction{T},
    g::MOI.VectorQuadraticFunction{T}) where T
    MOI.VectorQuadraticFunction([f.terms; operate_terms(op, g.affine_terms)],
                  operate_terms(op, g.quadratic_terms),
                  op.(f.constants, g.constants))
end

# Vector Quadratic +/- ...
function operate(op::Union{typeof(+), typeof(-)}, ::Type{T},
    f::MOI.VectorQuadraticFunction{T},
    g::VectorQuadraticLike{T}) where T
    operate!(op, T, copy(f), g)
end

function Base.:+(args::VectorLike{T}...) where T
    return operate(+, T, args...)
end
# Base.:+(α::Vector{T}, f::VectorLike{T}...) is too general as it also covers
# Base.:+(α::Vector) which is type piracy
function Base.:+(α::Vector{T}, f::VectorLike{T}, g::VectorLike{T}...) where T
    return operate(+, T, α, f, g...)
end
function Base.:+(f::VectorLike{T}, α::Vector{T}) where T
    return operate(+, T, f, α)
end
function Base.:-(args::VectorLike{T}...) where T
    return operate(-, T, args...)
end
function Base.:-(f::VectorLike{T}, α::Vector{T}) where T
    return operate(-, T, f, α)
end
function Base.:-(α::Vector{T}, f::VectorLike{T}) where T
    return operate(-, T, α, f)
end

####################################### * ######################################
function promote_operation(::typeof(*), ::Type{T}, ::Type{T},
                           ::Type{<:Union{MOI.SingleVariable,
                                          MOI.ScalarAffineFunction{T}}}) where T
    return MOI.ScalarAffineFunction{T}
end
function promote_operation(::typeof(*), ::Type{T},
                           ::Type{<:Union{MOI.SingleVariable,
                                          MOI.ScalarAffineFunction{T}}},
                           ::Type{T}) where T
    return MOI.ScalarAffineFunction{T}
end
function promote_operation(::typeof(*), ::Type{T}, ::Type{T},
                           ::Type{MOI.ScalarQuadraticFunction{T}}) where T
    return MOI.ScalarQuadraticFunction{T}
end
function promote_operation(::typeof(*), ::Type{T},
                           ::Type{MOI.ScalarQuadraticFunction{T}},
                           ::Type{T}) where T
    return MOI.ScalarQuadraticFunction{T}
end
function promote_operation(::typeof(*), ::Type{T},
                           ::Type{<:Union{MOI.SingleVariable,
                                          MOI.ScalarAffineFunction{T}}},
                           ::Type{<:Union{MOI.SingleVariable,
                                          MOI.ScalarAffineFunction{T}}}) where T
    return MOI.ScalarQuadraticFunction{T}
end

function operate!(::typeof(*), ::Type{T}, f::MOI.SingleVariable, α::T) where T
    return operate(*, T, α, f)
end
function operate(::typeof(*), ::Type{T}, α::T, f::MOI.SingleVariable) where T
    return MOI.ScalarAffineFunction{T}([MOI.ScalarAffineTerm(α, f.variable)], zero(T))
end
function operate(::typeof(*), ::Type{T}, α::T, f::MOI.VectorOfVariables) where T
    return MOI.VectorAffineFunction{T}(
        [MOI.VectorAffineTerm(i, MOI.ScalarAffineTerm(α, f.variables[i]))
         for i in eachindex(f.variables)], zeros(T, MOI.output_dimension(f)))
end
function operate(::typeof(*), ::Type{T},
                 f::Union{MOI.SingleVariable, MOI.VectorOfVariables},
                 α::T) where T
    return operate(*, T, α, f)
end

# `<:Number` is a workaround for https://github.com/JuliaOpt/MathOptInterface.jl/issues/980
function operate!(::typeof(*), ::Type{T},
                  f::Union{MOI.ScalarAffineFunction{T},
                           MOI.ScalarQuadraticFunction{T}}, α::T) where T<:Number
    map_terms!(term -> operate_term(*, term, α), f)
    f.constant *= α
    return f
end
function operate!(::typeof(*), ::Type{T},
                  f::Union{MOI.VectorAffineFunction{T},
                           MOI.VectorQuadraticFunction{T}}, α::T) where T
    map_terms!(term -> operate_term(*, term, α), f)
    rmul!(f.constants, α)
    return f
end
function operate(::typeof(*), ::Type{T}, α::T, f::TypedLike{T}) where T
    return operate!(*, T, copy(f), α)
end

function operate(::typeof(*), ::Type{T}, f::MOI.SingleVariable,
                 g::MOI.SingleVariable) where T
    return MOI.ScalarQuadraticFunction(
        MOI.ScalarAffineTerm{T}[],
        [MOI.ScalarQuadraticTerm(f.variable == g.variable ? 2one(T) : one(T),
                                 f.variable, g.variable)],
        zero(T))
end

function operate(::typeof(*), ::Type{T}, f::MOI.ScalarAffineFunction{T},
                 g::MOI.SingleVariable) where T
    if iszero(f.constant)
        aff_terms = MOI.ScalarAffineTerm{T}[]
    else
        aff_terms = [MOI.ScalarAffineTerm(f.constant, g.variable)]
    end
    quad_terms = map(t -> MOI.ScalarQuadraticTerm(
        t.variable_index == g.variable ? 2t.coefficient : t.coefficient,
        t.variable_index, g.variable), f.terms)
    return MOI.ScalarQuadraticFunction(aff_terms, quad_terms, zero(T))
end
function operate(::typeof(*), ::Type{T}, f::MOI.SingleVariable,
                 g::MOI.ScalarAffineFunction{T}) where T
    return operate(*, T, g, f)
end

function operate(::typeof(*), ::Type{T}, f::MOI.ScalarAffineFunction{T},
                 g::MOI.ScalarAffineFunction{T}) where T
    nfterms = length(f.terms)
    ngterms = length(g.terms)
    quad_terms = Vector{MOI.ScalarQuadraticTerm{T}}(undef, nfterms * ngterms)
    k = 0
    for t1 in f.terms
        for t2 in g.terms
            k += 1
            quad_terms[k] = operate_term(*, t1, t2)
        end
    end
    @assert k == length(quad_terms)
    if iszero(f.constant)
        if iszero(g.constant)
            aff_terms = MOI.ScalarAffineTerm{T}[]
        else
            aff_terms = operate_term.(*, g.constant, f.terms)
        end
    else
        if iszero(g.constant)
            aff_terms = operate_term.(*, f.constant, g.terms)
        else
            aff_terms = Vector{MOI.ScalarAffineTerm{T}}(undef,
                                                        nfterms + ngterms)
            map!(t -> operate_term(*, g.constant, t), aff_terms, f.terms)
            for i in 1:ngterms
                aff_terms[nfterms + i] = operate_term(*, f.constant, g.terms[i])
            end
        end
    end
    constant = f.constant * g.constant
    return MOI.ScalarQuadraticFunction(aff_terms, quad_terms, constant)
end

Base.:*(f::MOI.AbstractFunction) = f

# To avoid type piracy, we add at least one `ScalarLike` outside of the `...`.
function Base.:*(f::ScalarLike{T}, g::ScalarLike{T}, args::ScalarLike{T}...) where T
    return operate(*, T, f, g, args...)
end
function Base.:*(f::T, g::TypedLike{T}) where T
    return operate(*, T, f, g)
end
function Base.:*(f::Number, g::Union{MOI.SingleVariable, MOI.VectorOfVariables})
    return operate(*, typeof(f), f, g)
end
function Base.:*(f::TypedLike{T}, g::T) where T
    return operate(*, T, g, f)
end
function Base.:*(f::Union{MOI.SingleVariable, MOI.VectorOfVariables}, g::Number)
    return operate(*, typeof(g), f, g)
end
# Used by sparse matrix multiplication after
# https://github.com/JuliaLang/julia/pull/33205
function Base.:*(f::TypedLike, g::Bool)
    if g
        return MA.copy_if_mutable(f)
    else
        return zero(typeof(f))
    end
end
Base.:*(f::Bool, g::TypedLike) = g * f

function Base.:^(func::MOI.ScalarAffineFunction{T}, p::Integer) where T
    if iszero(p)
        return one(MOI.ScalarQuadraticFunction{T})
    elseif isone(p)
        return convert(MOI.ScalarQuadraticFunction{T}, func)
    elseif p == 2
        return func * func
    else
        throw(ArgumentError("Cannot take $(typeof(func)) to the power $p."))
    end
end
function Base.:^(func::MOI.ScalarQuadraticFunction{T}, p::Integer) where T
    if iszero(p)
        return one(MOI.ScalarQuadraticFunction{T})
    elseif isone(p)
        return MA.mutable_copy(func)
    else
        throw(ArgumentError("Cannot take $(typeof(func)) to the power $p."))
    end
end

function LinearAlgebra.dot(f::ScalarLike, g::ScalarLike)
    return f * g
end
function LinearAlgebra.dot(α::T, func::TypedLike{T}) where T
    return α * func
end
function LinearAlgebra.dot(func::TypedLike{T}, α::T) where T
    return func * α
end

LinearAlgebra.adjoint(f::ScalarLike) = f
LinearAlgebra.transpose(f::ScalarLike) = f
LinearAlgebra.symmetric_type(::Type{F}) where {F <: ScalarLike} = F
LinearAlgebra.symmetric(f::ScalarLike, ::Symbol) = f

####################################### / ######################################
function promote_operation(::typeof(/), ::Type{T},
                           ::Type{<:Union{MOI.SingleVariable,
                                          MOI.ScalarAffineFunction{T}}},
                           ::Type{T}) where T
    MOI.ScalarAffineFunction{T}
end
function promote_operation(::typeof(/), ::Type{T},
                           ::Type{MOI.ScalarQuadraticFunction{T}},
                           ::Type{T}) where T
    MOI.ScalarQuadraticFunction{T}
end

function operate!(::typeof(/), ::Type{T}, f::MOI.SingleVariable,
                  α::T) where T
    return operate(/, T, f, α)
end
function operate(::typeof(/), ::Type{T}, f::MOI.SingleVariable,
                 α::T) where T
    return MOI.ScalarAffineFunction{T}([MOI.ScalarAffineTerm(inv(α),
                                                             f.variable)],
                                       zero(T))
end

function operate!(::typeof(/), ::Type{T},
                  f::Union{MOI.ScalarAffineFunction{T},
                           MOI.ScalarQuadraticFunction{T}},
                  α::T) where T
    map_terms!(term -> operate_term(/, term, α), f)
    f.constant /= α
    return f
end

function operate!(::typeof(/), ::Type{T},
                  f::Union{MOI.VectorAffineFunction{T},
                           MOI.VectorQuadraticFunction{T}},
                  α::T) where T
    map_terms!(term -> operate_term(/, term, α), f)
    rmul!(f.constants, inv(α))
    return f
end

function operate!(::typeof(/), ::Type{T}, f::MOI.ScalarQuadraticFunction{T},
                  α::T) where T
    f.affine_terms .= operate_term.(/, f.affine_terms, α)
    f.quadratic_terms .= operate_term.(/, f.quadratic_terms, α)
    f.constant /= α
    return f
end

function operate(::typeof(/), ::Type{T}, f::TypedLike{T}, α::T) where T
    return operate!(/, T, copy(f), α)
end
function operate(::typeof(/), ::Type{T},
                 f::Union{MOI.SingleVariable, MOI.VectorOfVariables},
                 α::T) where T
    return operate(*, T, inv(α), f)
end

function Base.:/(f::TypedLike{T}, g::T) where T
    return operate(/, T, f, g)
end
function Base.:/(f::Union{MOI.SingleVariable, MOI.VectorOfVariables}, g::Number)
    return operate(/, typeof(g), f, g)
end

## sum
function operate(::typeof(sum), ::Type{T}, vis::Vector{MOI.VariableIndex}) where T
    return MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(one(T), vis), zero(T))
end

#################### Concatenation of MOI functions: `vcat` ####################
"""
    fill_vector(vector::Vector, ::Type{T}, fill_func::Function,
                dim_func::Function, funcs) where T

Fill the vector `vector` with
`fill_func(vector, vector_offset, output_offset, func)` for each function `func`
in `funcs` where `vector_offset` (resp. `output_offset`) is the sum of
`dim_func(T, func)` (resp. `output_dim(T, func)`) of previous functions of
`func`.

    fill_vector(vector::Vector, ::Type{T}, vector_offset::Int,
                     output_offset::Int, fill_func::Function,
                     dim_func::Function, funcs...) where T

Same than previous method but starting with possible nonzero `vector_offset` and
`output_offset`.
"""
function fill_vector end

function fill_vector(vector::Vector, ::Type{T}, fill_func::Function,
                     dim_func::Function, funcs) where T
    vector_offset = 0
    output_offset = 0
    for func in funcs
        fill_func(vector, vector_offset, output_offset, func)
        vector_offset += dim_func(T, func)
        output_offset += output_dim(T, func)
    end
    @assert length(vector) == vector_offset
end
function fill_vector(vector::Vector, ::Type, vector_offset::Int,
                     output_offset::Int, fill_func::Function,
                     dim_func::Function)
    @assert length(vector) == vector_offset
end
function fill_vector(vector::Vector, ::Type{T}, vector_offset::Int,
                     output_offset::Int, fill_func::Function,
                     dim_func::Function, func, funcs...) where T
    fill_func(vector, vector_offset, output_offset, func)
    fill_vector(vector, T, vector_offset + dim_func(T, func),
                output_offset + output_dim(T, func), fill_func, dim_func,
                funcs...)
end

function fill_variables(variables::Vector{MOI.VariableIndex}, offset::Int,
                        output_offset::Int, func::MOI.SingleVariable)
    variables[offset + 1] = func.variable
end

function fill_variables(variables::Vector{MOI.VariableIndex}, offset::Int,
                        output_offset::Int, func::MOI.VectorOfVariables)
    variables[offset .+ (1:length(func.variables))] .= func.variables
end

function promote_operation(::typeof(vcat), ::Type{T},
                           ::Type{<:Union{MOI.SingleVariable,
                                          MOI.VectorOfVariables}}...) where T
    return MOI.VectorOfVariables
end
function operate(::typeof(vcat), ::Type{T},
                 funcs::Union{MOI.SingleVariable,
                              MOI.VectorOfVariables}...) where T
    out_dim = sum(func -> output_dim(T, func), funcs)
    variables = Vector{MOI.VariableIndex}(undef, out_dim)
    fill_vector(variables, T, 0, 0, fill_variables, output_dim, funcs...)
    return MOI.VectorOfVariables(variables)
end

number_of_affine_terms(::Type{T}, ::T) where T = 0
number_of_affine_terms(::Type, ::SVF) = 1
number_of_affine_terms(::Type, f::VVF) = length(f.variables)
function number_of_affine_terms(
    ::Type{T}, f::Union{SAF{T}, VAF{T}}) where T
    return length(f.terms)
end
function number_of_affine_terms(
    ::Type{T}, f::Union{SQF{T}, VQF{T}}) where T
    return length(f.affine_terms)
end

number_of_quadratic_terms(::Type{T}, ::Union{T, SVF, VVF, SAF{T}, VAF{T}}) where {T} = 0
function number_of_quadratic_terms(
    ::Type{T}, f::Union{SQF{T}, VQF{T}}) where T
    return length(f.quadratic_terms)
end

function offset_term(t::MOI.ScalarAffineTerm, offset::Int)
    return MOI.VectorAffineTerm(offset + 1, t)
end
function offset_term(t::MOI.VectorAffineTerm, offset::Int)
    return MOI.VectorAffineTerm(offset + t.output_index, t.scalar_term)
end
function offset_term(t::MOI.ScalarQuadraticTerm, offset::Int)
    return MOI.VectorQuadraticTerm(offset + 1, t)
end
function offset_term(t::MOI.VectorQuadraticTerm, offset::Int)
    return MOI.VectorQuadraticTerm(offset + t.output_index, t.scalar_term)
end

function fill_terms(::Vector{MOI.VectorAffineTerm{T}}, ::Int, ::Int, ::T) where T
end
function fill_terms(terms::Vector{MOI.VectorAffineTerm{T}}, offset::Int,
                    output_offset::Int, func::SVF) where T
    terms[offset + 1] = offset_term(
        MOI.ScalarAffineTerm(one(T), func.variable), output_offset)
end
function fill_terms(terms::Vector{MOI.VectorAffineTerm{T}}, offset::Int,
                    output_offset::Int, func::VVF) where T
    n = number_of_affine_terms(T, func)
    terms[offset .+ (1:n)] .= MOI.VectorAffineTerm.(
        output_offset .+ (1:n), MOI.ScalarAffineTerm.(one(T), func.variables))
end
function fill_terms(terms::Vector{MOI.VectorAffineTerm{T}}, offset::Int,
                    output_offset::Int, func::Union{SAF{T}, VAF{T}}) where T
    n = number_of_affine_terms(T, func)
    terms[offset .+ (1:n)] .= offset_term.(func.terms, output_offset)
end
function fill_terms(terms::Vector{MOI.VectorAffineTerm{T}}, offset::Int,
                    output_offset::Int, func::Union{SQF{T}, VQF{T}}) where T
    n = number_of_affine_terms(T, func)
    terms[offset .+ (1:n)] .= offset_term.(func.affine_terms, output_offset)
end

function fill_terms(::Vector{MOI.VectorQuadraticTerm{T}}, ::Int, ::Int,
                    ::Union{T, SVF, VVF, SAF{T}, VAF{T}}) where T
end
function fill_terms(terms::Vector{MOI.VectorQuadraticTerm{T}}, offset::Int,
                    output_offset::Int, func::Union{SQF{T}, VQF{T}}) where T
    n = number_of_quadratic_terms(T, func)
    terms[offset .+ (1:n)] .= offset_term.(func.quadratic_terms, output_offset)
end

output_dim(::Type{T}, ::T) where T = 1
output_dim(::Type, func::MOI.AbstractFunction) = MOI.output_dimension(func)
function fill_constant(constant::Vector{T}, offset::Int,
                       output_offset::Int, func::T) where T
    constant[offset + 1] = func
end
function fill_constant(constant::Vector{T}, offset::Int,
                       output_offset::Int, func::Union{SVF, VVF}) where T
end
function fill_constant(constant::Vector{T}, offset::Int,
                       output_offset::Int, func::Union{SAF{T}, SQF{T}}) where T
    constant[offset + 1] = func.constant
end
function fill_constant(constant::Vector{T}, offset::Int,
                       output_offset::Int, func::Union{VAF{T}, VQF{T}}) where T
    n = MOI.output_dimension(func)
    constant[offset .+ (1:n)] .= func.constants
end

"""
    vectorize(funcs::AbstractVector{MOI.SingleVariable})

Returns the vector of scalar affine functions in the form of a
`MOI.VectorAffineFunction{T}`.
"""
function vectorize(funcs::AbstractVector{MOI.SingleVariable})
    vars = MOI.VariableIndex[func.variable for func in funcs]
    return MOI.VectorOfVariables(vars)
end

"""
    vectorize(funcs::AbstractVector{MOI.ScalarAffineFunction{T}}) where T

Returns the vector of scalar affine functions in the form of a
`MOI.VectorAffineFunction{T}`.
"""
function vectorize(funcs::AbstractVector{MOI.ScalarAffineFunction{T}}) where T
    nterms = mapreduce(func -> number_of_affine_terms(T, func), +, funcs, init=0)
    out_dim = mapreduce(func -> output_dim(T, func), +, funcs, init=0)
    terms = Vector{MOI.VectorAffineTerm{T}}(undef, nterms)
    constant = zeros(T, out_dim)
    fill_vector(terms, T, fill_terms, number_of_affine_terms, funcs)
    fill_vector(constant, T, fill_constant, output_dim, funcs)
    return MOI.VectorAffineFunction(terms, constant)
end

"""
    vectorize(funcs::AbstractVector{MOI.ScalarQuadraticFunction{T}}) where T

Returns the vector of scalar quadratic functions in the form of a
`MOI.VectorQuadraticFunction{T}`.
"""
function vectorize(funcs::AbstractVector{MOI.ScalarQuadraticFunction{T}}) where T
    num_affine_terms = mapreduce(func -> number_of_affine_terms(T, func), +, funcs, init=0)
    num_quadratic_terms = mapreduce(func -> number_of_quadratic_terms(T, func), +, funcs, init=0)
    out_dim = mapreduce(func -> output_dim(T, func), +, funcs, init=0)
    affine_terms = Vector{MOI.VectorAffineTerm{T}}(undef, num_affine_terms)
    quadratic_terms = Vector{MOI.VectorQuadraticTerm{T}}(undef, num_quadratic_terms)
    constant = zeros(T, out_dim)
    fill_vector(affine_terms, T, fill_terms, number_of_affine_terms, funcs)
    fill_vector(quadratic_terms, T, fill_terms, number_of_quadratic_terms, funcs)
    fill_vector(constant, T, fill_constant, output_dim, funcs)
    return MOI.VectorQuadraticFunction(affine_terms, quadratic_terms, constant)
end

function promote_operation(::typeof(vcat), ::Type{T}, ::Type{T}...) where T
    return Vector{T}
end
function promote_operation(::typeof(vcat), ::Type{T},
                           ::Type{<:Union{ScalarAffineLike{T}, VVF, VAF{T}}}...) where T
    return VAF{T}
end
function promote_operation(
    ::typeof(vcat), ::Type{T},
    ::Type{<:Union{ScalarQuadraticLike{T}, VVF, VAF{T}, VQF{T}}}...) where T
    return VQF{T}
end

function operate(::typeof(vcat), ::Type{T},
                 funcs::Union{ScalarAffineLike{T}, VVF, VAF{T}}...) where T
    nterms = sum(func -> number_of_affine_terms(T, func), funcs)
    out_dim = sum(func -> output_dim(T, func), funcs)
    terms = Vector{MOI.VectorAffineTerm{T}}(undef, nterms)
    constant = zeros(T, out_dim)
    fill_vector(terms, T, 0, 0, fill_terms, number_of_affine_terms, funcs...)
    fill_vector(constant, T, 0, 0, fill_constant, output_dim, funcs...)
    return VAF(terms, constant)
end
function operate(::typeof(vcat), ::Type{T},
                 funcs::Union{ScalarQuadraticLike{T}, VVF, VAF{T}, VQF{T}}...) where T
    num_affine_terms = sum(func -> number_of_affine_terms(T, func), funcs)
    num_quadratic_terms = sum(func -> number_of_quadratic_terms(T, func), funcs)
    out_dim = sum(func -> output_dim(T, func), funcs)
    affine_terms = Vector{MOI.VectorAffineTerm{T}}(undef, num_affine_terms)
    quadratic_terms = Vector{MOI.VectorQuadraticTerm{T}}(undef, num_quadratic_terms)
    constant = zeros(T, out_dim)
    fill_vector(affine_terms, T, 0, 0, fill_terms, number_of_affine_terms, funcs...)
    fill_vector(quadratic_terms, T, 0, 0, fill_terms, number_of_quadratic_terms, funcs...)
    fill_vector(constant, T, 0, 0, fill_constant, output_dim, funcs...)
    return VQF(affine_terms, quadratic_terms, constant)
end

# Similar to `eachscalar` but faster, see
# https://github.com/JuliaOpt/MathOptInterface.jl/issues/418
function scalarize(f::MOI.VectorOfVariables, ignore_constants::Bool = false)
    MOI.SingleVariable.(f.variables)
end
function scalarize(f::MOI.VectorAffineFunction{T}, ignore_constants::Bool = false) where T
    dimension = MOI.output_dimension(f)
    constants = ignore_constants ? zeros(T, dimension) : MOI.constant(f)
    counting = count_terms(dimension, f.terms)
    functions = MOI.ScalarAffineFunction{T}[
        MOI.ScalarAffineFunction{T}(MOI.ScalarAffineTerm{T}[], constants[i]) for i in 1:dimension]
    for i in 1:dimension
        sizehint!(functions[i].terms, counting[i])
    end
    for term in f.terms
        push!(functions[term.output_index].terms, term.scalar_term)
    end
    return functions
end
function scalarize(f::MOI.VectorQuadraticFunction{T}, ignore_constants::Bool = false) where T
    dimension = MOI.output_dimension(f)
    constants = ignore_constants ? zeros(T, dimension) : MOI.constant(f)
    counting_scalars = count_terms(dimension, f.affine_terms)
    counting_quadratics = count_terms(dimension, f.quadratic_terms)
    functions = MOI.ScalarQuadraticFunction{T}[
        MOI.ScalarQuadraticFunction{T}(MOI.ScalarAffineTerm{T}[],
                                       MOI.ScalarQuadraticTerm{T}[],
                                       constants[i]) for i in 1:dimension]
    for i in 1:dimension
        sizehint!(functions[i].affine_terms, counting_scalars[i])
        sizehint!(functions[i].quadratic_terms, counting_quadratics[i])
    end
    for term in f.affine_terms
        push!(functions[term.output_index].affine_terms, term.scalar_term)
    end
    for term in f.quadratic_terms
        push!(functions[term.output_index].quadratic_terms, term.scalar_term)
    end
    return functions
end

function count_terms(counting::Vector{<:Integer}, terms::Vector{T}) where T
    for term in terms
        counting[term.output_index] += 1
    end
    return nothing
end
function count_terms(dimension::I, terms::Vector{T}) where {I,T}
    counting = zeros(I, dimension)
    count_terms(counting, terms)
    return counting
end

tol_default(T::Type{<:Union{Integer, Rational}}) = zero(T)
tol_default(T::Type{<:AbstractFloat}) = sqrt(eps(T))
convert_approx(::Type{T}, func::T; kws...) where {T} = func
function convert_approx(::Type{MOI.SingleVariable}, func::MOI.ScalarAffineFunction{T};
                        tol=tol_default(T)) where {T}
    f = canonical(func)
    i = findfirst(t -> isapprox(t.coefficient, one(T), atol=tol), f.terms)
    if abs(f.constant) > tol || i === nothing ||
        any(j -> j != i && abs(f.terms[i]) > tol, eachindex(f.terms))
        throw(InexactError(:convert_approx, MOI.SingleVariable, func))
    end
    return MOI.SingleVariable(f.terms[i].variable_index)
end
function convert_approx(::Type{MOI.VectorOfVariables}, func::MOI.VectorAffineFunction{T};
                        tol=tol_default(T)) where {T}
    return MOI.VectorOfVariables([convert_approx(MOI.SingleVariable, f, tol=tol).variable
                                  for f in scalarize(func)])
end

function Base.zero(F::Type{<:TypedScalarLike{T}}) where T
    return convert(F, zero(T))
end
function Base.one(F::Type{<:TypedScalarLike{T}}) where T
    return convert(F, one(T))
end

Base.promote_rule(::Type{F}, ::Type{T}) where {T, F<:TypedScalarLike{T}} = F

ext: jl


text:
struct ExpectedMoreColumnsError <: Exception
    msg::String
end

struct TooManyColumnsError <: Exception
    msg::String
end

"""
`CSV.validate(fullpath::Union{AbstractString,IO}, sink::Type{T}=DataFrame, args...; kwargs...)` => `typeof(sink)`

`CSV.validate(fullpath::Union{AbstractString,IO}, sink::Data.Sink; kwargs...)` => `Data.Sink`

Takes the same positional & keyword arguments as [`CSV.read`](@ref), but provides detailed information as to why reading a csv file failed. Useful for cases where reading fails and it's not clear whether it's due to a row having too many columns, or wrong types, or what have you.

"""
function validate end

function validate(s::File)
    sch = Tables.schema(s) # size, header, types
    rows, cols = tablesize(s)
    types = sch.types
    for row = 1:rows
        rowstr = ""
        for col = 1:cols
            r = parsefield(s, parsingtype(types[col]), row, col)
            rowstr *= "$(col == 1 ? "" : ", ")$(r.result)"
            Parsers.ok(r.code) || throw(Error(Parsers.Error(s.io, r), row, col))
            if col < cols
                if eof(s.io)
                    throw(ExpectedMoreColumnsError("row=$row, col=$col: expected $cols columns, parsed $col, but parsing encountered unexpected end-of-file; parsed row: '$rowstr'"))
                elseif (r.code & Parsers.NEWLINE) > 0
                    throw(ExpectedMoreColumnsError("row=$row, col=$col: expected $cols columns, parsed $col, but parsing encountered unexpected newline character; parsed row: '$rowstr'"))
                end
            else
                if !eof(s.io) && !(r.code & Parsers.NEWLINE > 0)
                    throw(TooManyColumnsError("row=$row, col=$col: expected $cols columns then a newline or EOF; parsed row: '$rowstr'"))
                end
            end
        end
    end
end

function validate(fullpath::Union{AbstractString,IO}; append::Bool=false, transforms::Dict=Dict{Int,Function}(), kwargs...)
    validate(File(fullpath; kwargs...))
    return
end


text:
function echotime end

"""
    echotime(f::TimeProfile)

Get echo time ``TE`` of the time profile `f`, which is the end of the last characteristic
interval.
"""
echotime(f::TimeProfile) = intervals(f)[end]

"""
    echotime(grad::AbstractGradient)

Get echo time ``TE`` of `gradient`.
"""
echotime(gradient::GeneralGradient) = gradient.TE
echotime(gradient::ScalarGradient) = echotime(gradient.profile)


text:
using Roots, Dierckx, NLopt, CoolProp, PyPlot
import JSON
include("../functions/rel_perms.jl")
include("../functions/frac_flow_funcs.jl")

# read the data files
lowsal=JSON.parsefile("../juliaphreeqc/data/Lowsal_data.json")
paper_title = "fathi2010"
exp_type = "coreFlood"
exp_name = ["SF20-a", "SF20-b"] # formation water and sea water
recovery_data = "coreFloodData" # -> time(day), recovery
core_flood_condition = lowsal["fathi2010"]["coreFlood"]["SF20-a"]
L_core = core_flood_condition["L(mm)"]/1000
D_core = core_flood_condition["D(mm)"]/1000
poros  = core_flood_condition["porosity"]/100
swi     = core_flood_condition["swi"]/100
soi     = core_flood_condition["soi"]/100 # 1-swi
T_c     = core_flood_condition["T(C)"]
oil_type = core_flood_condition["oil"]
mu_oil_25 = lowsal["fathi2010"]["oil"][oil_type]["viscosity(cp)"]/100 # Pa.s
t_day = lowsal["fathi2010"]["coreFloodData"]["SF20-a"]["time(day)"]
R_oil = lowsal["fathi2010"]["coreFloodData"]["SF20-a"]["recovery"]/100
injection_rate_pv_day = core_flood_condition["injectionRate(PV/day)"]
recovery_final = core_flood_condition["recovery(%OOIP)"]/100
dp_exp = 0.0

so_final = (1-recovery_final)*soi # final oil saturation
pv_inj = t_day[end]*injection_rate_pv_day
t_sec = t_day*24*3600 # s
# conversion of the input data to the right unit
A_core  = π*D_core^2/4
PV_core = L_core*A_core*poros
u_inj   = injection_rate_pv_day/(24*3600)*PV_core/A_core # Darcy velocity
T_K = T_c + 273.15
perm_ave = 0.002e-12 # m^2 permeability
mu_water = PropsSI("viscosity", "T", T_K, "P", 5e5, "H2O")
# My formulation for oil viscosity at higher Temperature
# source: P.Daučík et al., Temperature Dependence of the Viscosity of Hydrocarbon Fractions
# Acta Chimica Slovaca, Vol.1, No. 1, 2008, 43 – 57
# log(mu/mu0) = B(1/T - 1/T0)
# find B for n-Dodecane, and use it for extrapolation
mu_1   = PropsSI("viscosity", "T", 300, "P", 5e5, "n-Dodecane")
mu_2   = PropsSI("viscosity", "T", T_K, "P", 5e5, "n-Dodecane")
B_visc = log(mu_1/mu_2)/(1/300-1/T_K)
mu_oil = exp(B_visc*(1/T_K-1/(25+273.15)))*mu_oil_25 # Pa.s

# Corrections to make the data consistent;
# The recovery line must be a straight line at the begining of the experiment.
# I shift the time axis a bit to make it a straight line:
# plot(t_sec, R_oil, "-o")
(a,b) = linreg(t_sec[2:6], R_oil[2:6])
# 0 = a + b*t_shift
t_shift = -a/b
t_sec_cor = [0; t_sec[2:end]-t_shift]
# plot(t_sec_cor, R_oil, "-o")

# define the objective functions
# (xt_shock_ww, sw_shock_ww, xt_prf_ww, sw_prf_ww, t_ww, p_inj_ww, R_oil_ww)=
# frac_flow_wf(muw=mu_water, muo=mu_oil, ut=inj_vel, phi=poros_res,
#   k=perm_res, swc=swc_ww, sor=sor_ww, kro0=kro0_ww, no=no_ww, krw0=krw0_ww,
#   nw=nw_ww, sw0=swc_ww, sw_inj=1.0, L=L_res, pv_inj=inj_pv)

# define the objective functions
# param = [sor, swc, kro0, krw0, no, nw]
function error_calc(param, muw, muo, ut, phi, k, sw0, sw_inj, L, pv_inj, t_exp, R_exp, dp_exp=0.0)
  (sor, swc, kro0, krw0, no, nw) = param
  (xt_shock, sw_shock, xt_prf, sw_prf, t, p_inj, R)=
  frac_flow_wf(muw=muw, muo=muo, ut=ut, phi=phi,
    k=k, swc=swc, sor=sor, kro0=kro0, no=no, krw0=krw0,
    nw=nw, sw0=sw0, sw_inj=sw_inj, L=L, pv_inj=pv_inj)
  R_int = Spline1D(t, R, k=1)
  dp_int = Spline1D(t, p_inj, k=1)
  error_vals_R = R_int(t_exp) - R_exp
  error_vals = error_vals_R
  if dp_exp != 0.0
      error_vals_dp = (dp_int(t_exp) - dp_exp)/maximum(dp_exp) # normalize
      error_vals = [error_vals_R; error_vals_dp;]
  end
  return error_vals
end

function error_calculation(param, param_ind, muw, muo, ut, phi, k, sw0, sw_inj, L, pv_inj,
  t_exp, R_exp, dp_exp; grad_calc=true)
  eps1 = 1e-8
  error_vals = error_calc(param, muw, muo, ut, phi, k, sw0, sw_inj, L, pv_inj, t_exp, R_exp, dp_exp)
  if grad_calc
    jac  = zeros(length(error_vals), length(param_ind))
  else
    jac = zeros(0)
  end
  if grad_calc
    for j in eachindex(param_ind)
      param2 = copy(param)
      param2[param_ind[j]]+=eps1
      error_val2 = error_calc(param2, muw, muo, ut, phi, k, sw0, sw_inj, L, pv_inj, t_exp, R_exp, dp_exp)
      jac[:, j] = (error_val2 - error_vals) / eps1
    end
  end
  return error_vals, jac
end

function objective_function(param::Vector{Float64}, grad::Vector{Float64}, param_all,
  param_ind, muw, muo, ut, phi, k, sw0, sw_inj, L, pv_inj, t_exp, R_exp, dp_exp; w=1.0)
  param2 = copy(param_all)
  param2[param_ind] = param
  grad_calc = false
  if length(grad)>0
    grad_calc = true
  end
  error_val, jac = error_calculation(param2, param_ind, muw, muo, ut, phi, k, sw0,
    sw_inj, L, pv_inj, t_exp, R_exp, dp_exp, grad_calc = grad_calc)
  obj_func_value = sum(abs2, w.*error_val)
  if length(grad)>0
      grad[:] = 2.0*sum(w.*error_val.*jac, 1)
  end
  return obj_func_value
end

function plot_results(param, muw, muo, ut, phi, k, sw0, sw_inj, L, pv_inj, t_exp, R_exp, dp_exp)
  (sor, swc, kro0, krw0, no, nw) = param
  (xt_shock, sw_shock, xt_prf, sw_prf, t, p_inj, R)=
  frac_flow_wf(muw=muw, muo=muo, ut=ut, phi=phi,
    k=k, swc=swc, sor=sor, kro0=kro0, no=no, krw0=krw0,
    nw=nw, sw0=sw0, sw_inj=sw_inj, L=L, pv_inj=pv_inj)
  R_int = Spline1D(t, R, k=1)
  dp_int = Spline1D(t, p_inj, k=1)
  error_vals = R_int(t_exp) - R_exp
  t_plot = linspace(minimum(t_exp), maximum(t_exp), 100)
  figure(1)
  plot(t_plot, R_int(t_plot), t_exp, R_exp, "o")
  xlabel("time [s]")
  ylabel("Recovery factor")
  figure(2)
  plot(t_plot, dp_int(t_plot))
  if dp_exp !=0.0
      plot(t_exp, dp_exp, "o")
  end
  xlabel("time [s]")
  ylabel("pressure drop [Pa]")
end

plot_quick = param -> plot_results(param, mu_water, mu_oil, u_inj, poros, perm_ave, swi, 1.0, L_core, pv_inj, t_sec_cor, R_oil, dp_exp)
# Optimization problem:

# parameter values
# [sor, swc, kro0, krw0, no, nw]
param_all = [so_final, 0.2, 0.5, 0.5, 2, 2.4]
param_ind = [1, 2,3,4,5,6]
w = ones(length(R_oil))
w[end]=2
w[end-1]=2
w[7] = 3
obj_fun = (param, grad) -> objective_function(param::Vector{Float64}, grad::Vector{Float64}, param_all,
  param_ind, mu_water, mu_oil, u_inj, poros, perm_ave, swi, 1.0, L_core, pv_inj, t_sec_cor, R_oil, dp_exp, w=w)

# x_init = [0.2, 0.8, 0.4, 2, 2]
# x_lb = [0.1, 0.1, 0.1, 0.5, 0.5]
# x_ub = [0.5, 1.0, 1.0, 5.0, 5.0]

x_init = copy(param_all)
x_lb = [0.08, 0.08, 0.05, 0.05, 1.0, 1.0]
x_ub = [0.5, 0.5, 1.0, 1.0, 4.0, 4.0]

# algorithms
# :LD_MMA
# :LN_COBYLA
# :LD_LBFGS
# :GN_DIRECT
# :GN_DIRECT_L
opt_alg=:LD_MMA

opt1 = Opt(opt_alg, length(x_init)) # choose the algorithm
lower_bounds!(opt1, x_lb)
upper_bounds!(opt1, x_ub)
ftol_rel!(opt1, 1e-10)
ftol_abs!(opt1, 1e-10)

min_objective!(opt1, obj_fun)
(fObjOpt, paramOpt, flag) = optimize(opt1, x_init)

param_plot = copy(param_all)
param_plot[param_ind] = paramOpt
plot_quick(param_plot)

max_stars_repo_name: simulkade/peteng

max_issues_repo_name: simulkade/peteng

max_forks_repo_name: simulkade/peteng


text:
#! /usr/bin/julia

# Rosetta Code, Set consolidation

function consolidate{T}(a::Array{Set{T},1})
    1 < length(a) || return a
    b = copy(a)
    c = Set{T}[]
    while 1 < length(b)
        x = shift!(b)
        cme = true
        for (i, y) in enumerate(b)
            !isempty(intersect(x, y)) || continue
            cme = false
            b[i] = union(x, y)
            break
        end
        !cme || push!(c, x)
    end
    push!(c, b[1])
    return c
end

p = Set(["A", "B"])
q = Set(["C", "D"])
r = Set(["B", "D"])
s = Set(["H", "I", "K"])
t = Set(["F", "G", "H"])

println("p = ", p)
println("q = ", q)
println("r = ", r)
println("s = ", s)
println("t = ", t)

println("consolidate([p, q]) =\n    ", consolidate([p, q]))
println("consolidate([p, r]) =\n    ", consolidate([p, r]))
println("consolidate([p, q, r]) =\n    ", consolidate([p, q, r]))
println("consolidate([p, q, r, s, t]) =\n    ",
        consolidate([p, q, r, s, t]))

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code