text:
module TestSubset

using Test
using DataFrames
using DataFramesMeta
using Statistics

const ≅ = isequal

@testset "subset" begin
    df = DataFrame(A = [1, 2, 3, missing], B = [2, 1, 2, 1])

    x = [2, 1, 0, 0]

    @test @subset(df, :A .> 1) == df[(df.A .> 1) .=== true,:]
    @test @subset(df, :B .> 1) == df[df.B .> 1,:]
    @test @subset(df, :A .> x) == df[(df.A .> x) .=== true,:]
    @test @subset(df, :B .> x) ≅ df[df.B .> x,:]
    @test @subset(df, :A .> :B, :B .> mean(:B)) == DataFrame(A = 3, B = 2)
    @test @subset(df, :A .> 1, :B .> 1) == df[map(&, df.A .> 1, df.B .> 1),:]
    @test @subset(df, :A .> 1, :A .< 4, :B .> 1) == df[map(&, df.A .> 1, df.A .< 4, df.B .> 1),:]

    @test @subset(df, :A .> 1).A isa Vector{Union{Missing, Int}}

    @test @subset(df, cols(:A) .> 1) == df[(df.A .> 1) .=== true,:]
    @test @subset(df, cols(:B) .> 1) == df[df.B .> 1,:]
    @test @subset(df, cols(:A) .> x) == df[(df.A .> x) .=== true,:]
    @test @subset(df, cols(:B) .> x) ≅ df[df.B .> x,:]
    @test @subset(df, cols(:A) .> :B, cols(:B) .> mean(:B)) == DataFrame(A = 3, B = 2)
    @test @subset(df, cols(:A) .> 1, :B .> 1) == df[map(&, df.A .> 1, df.B .> 1),:]
    @test @subset(df, cols(:A) .> 1, :A .< 4, :B .> 1) == df[map(&, df.A .> 1, df.A .< 4, df.B .> 1),:]

    @test @subset(df, :A .> 1, :A .<= 2) == DataFrame(A = 2, B = 1)

    subdf = @view df[df.B .== 2, :]

    @test @subset(subdf, :A .== 3) == DataFrame(A = 3, B = 2)
end

@testset "subset with :block" begin
    df = DataFrame(A = [1, 2, 3, missing], B = [2, 1, 2, 1])

    d = @subset df begin
        :A .> 1
        :B .> 1
    end
    @test d ≅ @subset(df, :A .> 1, :B .> 1)

    d = @subset df begin
        cols(:A) .> 1
        :B .> 1
    end
    @test d ≅ @subset(df, :A .> 1, :B .> 1)

    d = @subset df begin
        :A .> 1
        cols(:B) .> 1
    end
    @test d ≅ @subset(df, :A .> 1, :B .> 1)

    d = @subset df begin
        begin
            :A .> 1
        end
        :B .> 1
    end
    @test d ≅ @subset(df, :A .> 1, :B .> 1)

    d = @subset df begin
        :A .> 1
        @. :B > 1
    end
    @test d ≅ @subset(df, :A .> 1, :B .> 1)
end

@testset "subset!" begin
    df = DataFrame(A = [1, 2, 3, missing], B = [2, 1, 2, 1])

    x = [2, 1, 0, 0]

    df2 = copy(df)
    @test @subset!(df2, :A .> 1) === df2
    @test df2 == df[(df.A .> 1) .=== true,:]

    @test @subset!(copy(df), :B .> 1) == df[df.B .> 1,:]
    @test @subset!(copy(df), :A .> x) == df[(df.A .> x) .=== true,:]
    @test @subset!(copy(df), :B .> x) ≅ df[df.B .> x,:]
    @test @subset!(copy(df), :A .> :B, :B .> mean(:B)) == DataFrame(A = 3, B = 2)
    @test @subset!(copy(df), :A .> 1, :B .> 1) == df[map(&, df.A .> 1, df.B .> 1),:]
    @test @subset!(copy(df), :A .> 1, :A .< 4, :B .> 1) == df[map(&, df.A .> 1, df.A .< 4, df.B .> 1),:]

    @test @subset!(copy(df), :A .> 1).A isa Vector{Union{Missing, Int}}

    @test @subset!(copy(df), cols(:A) .> 1) == df[(df.A .> 1) .=== true,:]
    @test @subset!(copy(df), cols(:B) .> 1) == df[df.B .> 1,:]
    @test @subset!(copy(df), cols(:A) .> x) == df[(df.A .> x) .=== true,:]
    @test @subset!(copy(df), cols(:B) .> x) ≅ df[df.B .> x,:]
    @test @subset!(copy(df), cols(:A) .> :B, cols(:B) .> mean(:B)) == DataFrame(A = 3, B = 2)
    @test @subset!(copy(df), cols(:A) .> 1, :B .> 1) == df[map(&, df.A .> 1, df.B .> 1),:]
    @test @subset!(copy(df), cols(:A) .> 1, :A .< 4, :B .> 1) == df[map(&, df.A .> 1, df.A .< 4, df.B .> 1),:]

    @test @subset!(copy(df), :A .> 1, :A .<= 2) == DataFrame(A = 2, B = 1)

    subdf = @view df[df.B .== 2, :]

    @test @subset!(copy(subdf), :A .== 3) == DataFrame(A = 3, B = 2)
end

@testset "subset! with :block" begin
    df = DataFrame(A = [1, 2, 3, missing], B = [2, 1, 2, 1])

    d = @subset! copy(df) begin
        :A .> 1
        :B .> 1
    end
    @test d ≅ @subset!(copy(df), :A .> 1, :B .> 1)

    d = @subset! copy(df) begin
        cols(:A) .> 1
        :B .> 1
    end
    @test d ≅ @subset!(copy(df), :A .> 1, :B .> 1)

    d = @subset! copy(df) begin
        :A .> 1
        cols(:B) .> 1
    end
    @test d ≅ @subset!(copy(df), :A .> 1, :B .> 1)

    d = @subset! copy(df) begin
        begin
            :A .> 1
        end
        :B .> 1
    end
    @test d ≅ @subset!(copy(df), :A .> 1, :B .> 1)

    d = @subset! copy(df) begin
        :A .> 1
        @. :B > 1
    end
    @test d ≅ @subset!(copy(df), :A .> 1, :B .> 1)
end

@testset "@subset with a grouped data frame" begin
    df = DataFrame(
        g = [1, 1, 1, 2, 2],
        i = 1:5,
        t = ["a", "b", "c", "c", "e"],
        y = [:v, :w, :x, :y, :z],
        c = [:g, :quote, :body, :transform, missing]
    )

    gd = groupby(df, :g)

    @test @subset(gd, :i .== first(:i)) ≅ df[[1, 4], :]
    @test @subset(gd, cols(:i) .> mean(cols(:i)), :t .== "c") ≅ df[[3], :]
    @test @subset(gd, :c .== :g) ≅ df[[], :]
end

@testset "@subset! with a grouped data frame" begin
    df = DataFrame(
        g = [1, 1, 1, 2, 2],
        i = 1:5,
        t = ["a", "b", "c", "c", "e"],
        y = [:v, :w, :x, :y, :z],
        c = [:g, :quote, :body, :transform, missing]
    )

    @test @subset!(groupby(copy(df), :g), :i .== first(:i)) ≅ df[[1, 4], :]
    @test @subset!(groupby(copy(df), :g), cols(:i) .> mean(cols(:i)), :t .== "c") ≅ df[[3], :]
    @test @subset!(groupby(copy(df), :g), :c .== :g) ≅ df[[], :]
end

end # module


text:
using Documenter, PenPlots

makedocs(sitename = "PenPlots.jl")

deploydocs(repo = "github.com/paulgb/PenPlots.jl.git")

ext: jl


text:
using GenericSchur;
using BivMatFun;
import BivMatFun.mp;

const data_directory = string(@__DIR__(), "/data")

function is_ci_test()
    return haskey(ENV, "BIVMATFUN_CI")
end

function kahan(::Type{T}, m::Integer, n::Integer, theta, pert) where T
    theta = convert(T, theta)
    pert = convert(T, pert)
    s = sin(theta)
    c = cos(theta)
    dim = min(m,n)
    A = zeros(T, m, n)
    for i = 1:m, j = 1:n
        i > dim ? A[i,j] = zero(T) :
        i > j   ? A[i,j] = zero(T) :
        i==j    ? A[i,j] = s^(i-1)+pert*eps(T)*(m-i+1) : A[i,j] = -c*s^(i-1)
    end
    return A
end
kahan(::Type{T}, n::Integer, theta, pert) where T = kahan(T, n, n, theta, pert)
kahan(::Type{T}, n::Integer) where T = kahan(T, n, n, 1.2, 25.)
kahan(args...) = kahan(Float64, args...)
kahan(::Type, args...) = throw(MethodError(kahan, Tuple(args)))

function sampling(::Type{T}, x::Vector) where T
    n = length(x)
    A = zeros(T, n, n)
    for j = 1:n, i = 1:n
        if i != j
            A[i,j] = x[i] / (x[i] - x[j])
        end
    end
    d = sum(A, dims=2)
    A = A + diagm(0 => d[:])
    return A
end
function sampling(::Type{T}, n::Integer) where T
    p = T[1:n;] / n
    return sampling(T, p)
end
sampling(args...) = sampling(Float64, args...)
sampling(::Type, args...) = throw(MethodError(sampling, Tuple(args)))

function grcar(::Type{T}, n::Integer, k::Integer = 3) where T
    # Compute grcar matrix
    G = tril(triu(ones(T, n,n)), min(k, n-1)) - diagm(-1 => ones(T, n-1))
    return G
end
grcar(args...) = grcar(Float64, args...)
grcar(::Type, args...) = throw(MethodError(grcar, Tuple(args)))

function construct_random_benchmark(f, n)
    d = 64;

    VA = mp(randn(n, n).+ 1im * randn(n, n), d);
    DA = mp(1.0 .+ rand(n) .+ 1im * randn(n), d)
    VB = mp(randn(n, n).+ 1im * randn(n, n), d);
    DB = mp(1.0 .+ rand(n).+ 1im * randn(n), d)
    C  = mp(randn(n, n).+ 1im * randn(n, n), d)

    A = convert(Matrix{ComplexF64}, VA * diagm(DA) / VA);
    B = convert(Matrix{ComplexF64}, VB * diagm(DB) / VB);

    FF = convert(Matrix{ComplexF64}, VA * BivMatFun.diag_fun(f, DA, DB, C) / VB)
    C = convert(Matrix{ComplexF64}, VA * C / VB)

    return A, B, C, FF
end

function evaluate_reference_solution(f, A, B, C)
    dgts = 128

    return BivMatFun.with_digits(dgts) do
        Y2 = similar(C)

        A = mp(A, dgts);
        B = mp(B, dgts);
        C = mp(C, dgts);

        FA = schur(A);
        FB = schur(B);

        C = FA.Z' * C * FB.Z;

        # size of the perturbation, make sure it's sufficiently far given the
        # current precision level.
        ep = 10.0^(-dgts / 2.0);

        TA = FA.T + diagm(randn(size(FA.T, 1))) * ep * norm(FA.T, Inf);
        TB = FB.T + diagm(randn(size(FB.T, 1))) * ep * norm(FB.T, Inf);
        Y = BivMatFun.diag_fun(f, TA, TB, C);

        return convert(Matrix{Complex{Float64}}, FA.Z * Y * FB.Z');
    end
end

max_stars_repo_name: numpi/BivMatFun

max_issues_repo_name: numpi/BivMatFun

max_forks_repo_name: numpi/BivMatFun


text:
# *********************************************************************************
# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification,
# are permitted provided that the following conditions are met:
#
# Redistributions of source code must retain the above copyright notice, this list
# of conditions and the following disclaimer.
#
# Redistributions in binary form must reproduce the above copyright notice, this
# list of conditions and the following disclaimer in the documentation and/or other
# materials provided with the distribution.
#
# Neither the name of the copyright holder nor the names of its contributors may be
# used to endorse or promote products derived from this software without specific
# prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
# OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
# OF THE POSSIBILITY OF SUCH DAMAGE.
# *********************************************************************************
"""
    MPCElectricLoad

    Base.@kwdef struct MPCElectricLoad
        loads_kw::Array{Real,1}
        critical_loads_kw::Union{Missing, Array{Real,1}} = missing
    end
"""
Base.@kwdef struct MPCElectricLoad
    loads_kw::Array{Real,1}
    critical_loads_kw::Union{Missing, Array{Real,1}} = missing
end

"""
    MPCFinancial

    Base.@kwdef struct MPCFinancial
        VoLL::Union{Array{R,1}, R} where R<:Real = 1.00
    end
"""
Base.@kwdef struct MPCFinancial
    VoLL::Union{Array{R,1}, R} where R<:Real = 1.00
end

"""
    MPCPV
```julia
Base.@kwdef struct MPCPV
    name::String="PV"
    size_kw::Real = 0
    prod_factor_series_kw::Union{Missing, Array{Real,1}} = missing
end
```
"""
Base.@kwdef struct MPCPV
    name::String="PV"
    size_kw::Real = 0
    prod_factor_series_kw::Union{Missing, Array{Real,1}} = missing
end

struct MPCElectricTariff
    energy_rates::AbstractVector{Float64}

    monthly_demand_rates::AbstractVector{Float64}
    time_steps_monthly::Array{Array{Int64,1},1}  # length = 0 or 12
    monthly_previous_peak_demands::AbstractVector{Float64}

    tou_demand_rates::AbstractVector{Float64}
    tou_demand_ratchet_timesteps::Array{Array{Int64,1},1}  # length = n_tou_demand_ratchets
    tou_previous_peak_demands::AbstractVector{Float64}

    fixed_monthly_charge::Float64
    annual_min_charge::Float64
    min_monthly_charge::Float64

    export_rates::DenseAxisArray{Array{Float64,1}}
    export_bins::Array{Symbol,1}
end

"""
    MPCElectricTariff(d::Dict)

function for parsing user inputs into
```julia
    struct MPCElectricTariff
        monthly_previous_peak_demands::Array{Float64,1}
        energy_rates::Array{Float64,1}

        monthly_demand_rates::Array{Float64,1}
        time_steps_monthly::Array{Array{Int64,1},1}  # length = 0 or 12

        tou_demand_rates::Array{Float64,1}
        tou_demand_ratchet_timesteps::Array{Array{Int64,1},1}  # length = n_tou_demand_ratchets
        tou_previous_peak_demands::Array{Float64,1}

        fixed_monthly_charge::Float64
        annual_min_charge::Float64
        min_monthly_charge::Float64

        export_rates::DenseAxisArray{Array{Float64,1}}
        export_bins::Array{Symbol,1}
    end
```

Keys for `d` include:

  - `energy_rates`

    - REQUIRED
    - must have length equal to `ElectricLoad.loads_kw`

  - `monthly_demand_rates`

    - default = [0]

  - `time_steps_monthly`

    - array of arrays for integer time steps that the `monthly_demand_rates` apply to
    - default = [collect(1:length(energy_rates))]

  - `monthly_previous_peak_demands`

    - default = [0]

  - `tou_demand_rates`

    - an array of time-of-use demand rates
    - must have length equal to `tou_demand_timesteps`
    - default = []

  - `tou_demand_timesteps`

    - an array of arrays for the integer time steps that apply to the `tou_demand_rates`
    - default = []

  - `tou_previous_peak_demands`

    - an array of the previous peak demands set in each time-of-use demand period
    - must have length equal to `tou_demand_timesteps`
    - default = []

  - `net_metering`

    - boolean, if `true` then customer DER export is compensated at the `energy_rates`

  - `export_rates`

    - can be a <:Real or Array{<:Real, 1}, or not provided
    - if provided, customer DER export is compensated at the `export_rates`

NOTE: if both `net_metering=true` and `export_rates` are provided then the model can choose from either option.
"""
function MPCElectricTariff(d::Dict)

    energy_rates = d["energy_rates"]

    monthly_demand_rates = get(d, "monthly_demand_rates", [0.0])
    time_steps_monthly = get(d, "time_steps_monthly", [collect(1:length(energy_rates))])
    monthly_previous_peak_demands = get(d, "monthly_previous_peak_demands", [0.0])

    tou_demand_rates = get(d, "tou_demand_rates", Float64[])
    tou_demand_timesteps = get(d, "tou_demand_timesteps", [])
    tou_previous_peak_demands = get(d, "tou_previous_peak_demands", Float64[])
    @assert length(tou_demand_rates) == length(tou_demand_timesteps) == length(tou_previous_peak_demands)

    # TODO can remove these inputs?
    fixed_monthly_charge = 0.0
    annual_min_charge = 0.0
    min_monthly_charge = 0.0

    # TODO handle tiered rates
    export_bins = [:NEM, :WHL]
    nem_rate = []
    NEM = get(d, "net_metering", false)
    if NEM
        nem_rate = [-0.999 * x for x in energy_rates]
    end
    # export_rates can be a <:Real or Array{<:Real, 1}, or not provided
    export_rates = get(d, "export_rates", nothing)
    if !isnothing(export_rates)
        export_rates = convert(Vector{Real}, export_rates)
    end
    whl_rate = create_export_rate(export_rates, length(energy_rates), 1)

    if !NEM & (sum(whl_rate) >= 0)
        export_rates = DenseAxisArray{Array{Float64,1}}(undef, [])
        export_bins = Symbol[]
    elseif !NEM
        export_bins = [:WHL]
        export_rates = DenseAxisArray([whl_rate], export_bins)
    elseif (sum(whl_rate) >= 0)
        export_bins = [:NEM]
        export_rates = DenseAxisArray([nem_rate], export_bins)
    else
        export_bins = [:NEM, :WHL]  # TODO add :EXC to align with REopt Lite API
        export_rates = DenseAxisArray([nem_rate, whl_rate], export_bins)
    end

    MPCElectricTariff(
        energy_rates,
        monthly_demand_rates,
        time_steps_monthly,
        monthly_previous_peak_demands,
        tou_demand_rates,
        tou_demand_timesteps,
        tou_previous_peak_demands,
        fixed_monthly_charge,
        annual_min_charge,
        min_monthly_charge,
        export_rates,
        export_bins
    )
end

"""
    MPCElecStorage
```julia
Base.@kwdef struct MPCElecStorage
    size_kw::Float64
    size_kwh::Float64
    charge_efficiency::Float64
    discharge_efficiency::Float64
    soc_min_pct::Float64
    soc_init_pct::Float64
    can_grid_charge::Bool = true
    grid_charge_efficiency::Float64
end
```
"""
Base.@kwdef struct MPCElecStorage
    size_kw::Float64
    size_kwh::Float64
    charge_efficiency::Float64
    discharge_efficiency::Float64
    soc_min_pct::Float64
    soc_init_pct::Float64
    can_grid_charge::Bool = true
    grid_charge_efficiency::Float64
end

Base.@kwdef struct MPCStorage
    types::Array{Symbol,1} = [:elec]
    size_kw::Dict{Symbol, Float64}
    size_kwh::Dict{Symbol, Float64}
    charge_efficiency::Dict{Symbol, Float64} = Dict(:elec => 0.96 * 0.975^2)
    discharge_efficiency::Dict{Symbol, Float64} = Dict(:elec => 0.96 * 0.975^2)
    soc_min_pct::Dict{Symbol, Float64} = Dict(:elec => 0.2)
    soc_init_pct::Dict{Symbol, Float64} = Dict(:elec => 0.5)
    can_grid_charge::Array{Symbol,1} = [:elec]
    grid_charge_efficiency::Dict{Symbol, Float64} = Dict(:elec => 0.96 * 0.975^2)
end

"""

NOTE: d must have symbolic keys
"""
function MPCStorage(d::Dict)
    d2 = Dict()
    d2[:can_grid_charge] = get(d, :can_grid_charge, false) ? [:elec] : Symbol[]
    if haskey(d, :can_grid_charge)
        pop!(d, :can_grid_charge)
    end
    # have to convert to all d values to DenseAxisArray's with storage type as Axis
    # (only modeling elec storage in MPC for now)
    for (k,v) in d
        d2[k] = Dict(:elec => convert(Float64, v))
    end

    return MPCStorage(; d2...)
end

"""
    MPCGenerator

struct with inner constructor:
```julia
function MPCGenerator(;
    size_kw::Real,
    fuel_cost_per_gallon::Float64 = 3.0,
    fuel_slope_gal_per_kwh::Float64 = 0.076,
    fuel_intercept_gal_per_hr::Float64 = 0.0,
    fuel_avail_gal::Float64 = 660.0,
    min_turn_down_pct::Float64 = 0.0,  # TODO change this to non-zero value
    only_runs_during_grid_outage::Bool = true,
    sells_energy_back_to_grid::Bool = false,
    om_cost_per_kwh::Float64=0.0,
    )
```
"""
struct MPCGenerator <: AbstractGenerator
    size_kw
    max_kw
    fuel_cost_per_gallon
    fuel_slope_gal_per_kwh
    fuel_intercept_gal_per_hr
    fuel_avail_gal
    min_turn_down_pct
    only_runs_during_grid_outage
    sells_energy_back_to_grid
    om_cost_per_kwh

    function MPCGenerator(;
        size_kw::Real,
        fuel_cost_per_gallon::Float64 = 3.0,
        fuel_slope_gal_per_kwh::Float64 = 0.076,
        fuel_intercept_gal_per_hr::Float64 = 0.0,
        fuel_avail_gal::Float64 = 660.0,
        min_turn_down_pct::Float64 = 0.0,  # TODO change this to non-zero value
        only_runs_during_grid_outage::Bool = true,
        sells_energy_back_to_grid::Bool = false,
        om_cost_per_kwh::Float64=0.0,
        )

        max_kw = size_kw

        new(
            size_kw,
            max_kw,
            fuel_cost_per_gallon,
            fuel_slope_gal_per_kwh,
            fuel_intercept_gal_per_hr,
            fuel_avail_gal,
            min_turn_down_pct,
            only_runs_during_grid_outage,
            sells_energy_back_to_grid,
            om_cost_per_kwh,
        )
    end
end

max_stars_repo_name: NLaws/REoptLite

max_issues_repo_name: NLaws/REoptLite

max_forks_repo_name: NLaws/REoptLite


text:
"""
    simulate
    - solves 1-step feasibility problem for T time steps
    - initial configurations: q1, q2 (note this can encode initial velocity)
    - time step: h
"""
function simulate(q1, q2, T, h)
    println("simulation")

    # initialize histories
    q = [q1, q2]
    n = [0.0]
    b = [zeros(2)]
    Δq1 = []
    Δq2 = []
    Δu1 = []

    u1 = zeros(3)

    # step
    for t = 1:T
        println("   t = $t")
        q3, n1, b1, _Δq1, _Δq2, _Δu1, status = step(q[end-1], q[end], u1, h)

        if !status
            @error "failed step (t = $t)"
            return q, n, b
        else
            push!(q, q3)
            push!(n, n1)
            push!(b, b1)
            push!(Δq1, _Δq1)
            push!(Δq2, _Δq2)
            push!(Δu1, _Δu1)
        end
    end

    return q, n, b, Δq1, Δq2, Δu1
end


text:
module TestMOIWrapper

using Gurobi
using Random
using Test

const MOI  = Gurobi.MOI
const MOIT = MOI.Test

const GRB_ENV = isdefined(Main, :GRB_ENV) ? Main.GRB_ENV : Gurobi.Env()

const OPTIMIZER = MOI.Bridges.full_bridge_optimizer(
    begin
        model = Gurobi.Optimizer(GRB_ENV)
        MOI.set(model, MOI.Silent(), true)
        # We set `DualReductions = 0` so that we never return
        # `INFEASIBLE_OR_UNBOUNDED`.
        MOI.set(model, MOI.RawParameter("DualReductions"), 0)
        MOI.set(model, MOI.RawParameter("QCPDual"), 1)
        MOI.set(model, MOI.RawParameter("InfUnbdInfo"), 1)
        MOI.set(model, MOI.RawParameter("NonConvex"), 2)
        model
    end,
    Float64
)
const CONFIG = MOIT.TestConfig()

function test_basic_constraint_tests()
    MOIT.basic_constraint_tests(OPTIMIZER, CONFIG; exclude = [
        (MOI.VectorOfVariables, MOI.SecondOrderCone),
        (MOI.VectorOfVariables, MOI.RotatedSecondOrderCone),
        (MOI.VectorOfVariables, MOI.GeometricMeanCone),
        (MOI.VectorAffineFunction{Float64}, MOI.SecondOrderCone),
        (MOI.VectorAffineFunction{Float64}, MOI.RotatedSecondOrderCone),
        (MOI.VectorAffineFunction{Float64}, MOI.GeometricMeanCone),
        (MOI.VectorQuadraticFunction{Float64}, MOI.SecondOrderCone),
        (MOI.VectorQuadraticFunction{Float64}, MOI.RotatedSecondOrderCone),
        (MOI.VectorQuadraticFunction{Float64}, MOI.GeometricMeanCone),
    ])
    # TODO(odow): bugs deleting SOC variables. See also the
    # `delete_soc_variables` test.
    MOIT.basic_constraint_tests(
        OPTIMIZER,
        CONFIG;
        include = [
            (MOI.VectorOfVariables, MOI.SecondOrderCone),
            (MOI.VectorOfVariables, MOI.RotatedSecondOrderCone),
            (MOI.VectorOfVariables, MOI.GeometricMeanCone),
            (MOI.VectorAffineFunction{Float64}, MOI.SecondOrderCone),
            (MOI.VectorAffineFunction{Float64}, MOI.RotatedSecondOrderCone),
            (MOI.VectorAffineFunction{Float64}, MOI.GeometricMeanCone),
            (MOI.VectorQuadraticFunction{Float64}, MOI.SecondOrderCone),
            (MOI.VectorQuadraticFunction{Float64}, MOI.RotatedSecondOrderCone),
            (MOI.VectorQuadraticFunction{Float64}, MOI.GeometricMeanCone),
        ],
        delete = false
    )
end

function test_unittest()
    MOIT.unittest(OPTIMIZER, MOIT.TestConfig(atol=1e-6), [
        # TODO(odow): bug! We can't delete a vector of variables if one is in
        # a second order cone.
        "delete_soc_variables",
    ])
end

function test_modificationtest()
    MOIT.modificationtest(OPTIMIZER, CONFIG)
end

function test_contlineartest()
    MOIT.contlineartest(OPTIMIZER, MOIT.TestConfig(basis = true))
end

function test_contquadratictest()
    MOIT.contquadratictest(OPTIMIZER, MOIT.TestConfig(atol=1e-3, rtol=1e-3))
end

function test_conictest()
    MOIT.lintest(OPTIMIZER, CONFIG)
    MOIT.soctest(OPTIMIZER, MOIT.TestConfig(duals = false, atol=1e-3), ["soc3"])
    MOIT.soc3test(
        OPTIMIZER,
        MOIT.TestConfig(duals = false, infeas_certificates = false, atol = 1e-3)
    )
    MOIT.rsoctest(OPTIMIZER, MOIT.TestConfig(duals = false, atol=5e-3))
    MOIT.geomeantest(OPTIMIZER, MOIT.TestConfig(duals = false, atol=1e-3))
end

function test_intlinear()
    MOIT.intlineartest(OPTIMIZER, CONFIG)
end

function test_solvername()
    @test MOI.get(OPTIMIZER, MOI.SolverName()) == "Gurobi"
end

function test_default_objective_test()
    MOIT.default_objective_test(OPTIMIZER)
end

function test_default_status_test()
    MOIT.default_status_test(OPTIMIZER)
end

function test_nametest()
    MOIT.nametest(OPTIMIZER)
end

function test_validtest()
    MOIT.validtest(OPTIMIZER)
end

function test_emptytest()
    MOIT.emptytest(OPTIMIZER)
end

function test_orderedindicestest()
    MOIT.orderedindicestest(OPTIMIZER)
end

function test_copytest()
    MOIT.copytest(
        OPTIMIZER,
        MOI.Bridges.full_bridge_optimizer(Gurobi.Optimizer(GRB_ENV), Float64)
    )
end

function test_scalar_function_constant_not_zero()
    MOIT.scalar_function_constant_not_zero(OPTIMIZER)
end

function test_start_values_test()
    model = Gurobi.Optimizer(GRB_ENV)
    MOI.set(model, MOI.Silent(), true)
    x = MOI.add_variables(model, 2)
    @test MOI.supports(model, MOI.VariablePrimalStart(), MOI.VariableIndex)
    @test MOI.get(model, MOI.VariablePrimalStart(), x[1]) === nothing
    @test MOI.get(model, MOI.VariablePrimalStart(), x[2]) === nothing
    Gurobi._update_if_necessary(model)
    p = Ref{Cdouble}()
    Gurobi.GRBgetdblattrelement(model.inner, "Start", Cint(0), p)
    @test p[] == Gurobi.GRB_UNDEFINED
    MOI.set(model, MOI.VariablePrimalStart(), x[1], 1.0)
    MOI.set(model, MOI.VariablePrimalStart(), x[2], nothing)
    @test MOI.get(model, MOI.VariablePrimalStart(), x[1]) == 1.0
    @test MOI.get(model, MOI.VariablePrimalStart(), x[2]) === nothing
    Gurobi._update_if_necessary(model)
    Gurobi.GRBgetdblattrelement(model.inner, "Start", Cint(1), p)
    @test p[] == Gurobi.GRB_UNDEFINED
    MOI.optimize!(model)
    @test MOI.get(model, MOI.ObjectiveValue()) == 0.0
    # We don't support ConstraintDualStart or ConstraintPrimalStart yet.
    # @test_broken MOIT.start_values_test(Gurobi.Optimizer(GRB_ENV), OPTIMIZER)
end

function test_supports_constrainttest()
    # supports_constrainttest needs VectorOfVariables-in-Zeros,
    # MOIT.supports_constrainttest(Gurobi.Optimizer(GRB_ENV), Float64, Float32)
    # but supports_constrainttest is broken via bridges:
    MOI.empty!(OPTIMIZER)
    MOI.add_variable(OPTIMIZER)
    @test  MOI.supports_constraint(OPTIMIZER, MOI.SingleVariable, MOI.EqualTo{Float64})
    @test  MOI.supports_constraint(OPTIMIZER, MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64})
    # This test is broken for some reason:
    @test_broken !MOI.supports_constraint(OPTIMIZER, MOI.ScalarAffineFunction{Int}, MOI.EqualTo{Float64})
    @test !MOI.supports_constraint(OPTIMIZER, MOI.ScalarAffineFunction{Int}, MOI.EqualTo{Int})
    @test !MOI.supports_constraint(OPTIMIZER, MOI.SingleVariable, MOI.EqualTo{Int})
    @test  MOI.supports_constraint(OPTIMIZER, MOI.VectorOfVariables, MOI.Zeros)
    @test !MOI.supports_constraint(OPTIMIZER, MOI.VectorOfVariables, MOI.EqualTo{Float64})
    @test !MOI.supports_constraint(OPTIMIZER, MOI.SingleVariable, MOI.Zeros)
    @test !MOI.supports_constraint(OPTIMIZER, MOI.VectorOfVariables, MOIT.UnknownVectorSet)
end

function test_set_lower_bound_twice()
    MOIT.set_lower_bound_twice(OPTIMIZER, Float64)
end

function test_set_upper_bound_twice()
    MOIT.set_upper_bound_twice(OPTIMIZER, Float64)
end

function test_User_limit_handling_issue_140()
    # Verify that we return the correct status codes when a mixed-integer
    # problem has been solved to a *feasible* but not necessarily optimal
    # solution. To do that, we will set up an intentionally dumbed-down
    # Gurobi Gurobi.Optimizer (with all heuristics and pre-solve turned off) and
    # ask it to solve a classic knapsack problem. Setting SolutionLimit=1
    # forces the solver to return after its first feasible MIP solution,
    # which tests the right part of the code without relying on potentially
    # flaky or system-dependent time limits.
    m = Gurobi.Optimizer(GRB_ENV)
    MOI.set(m, MOI.Silent(), true)
    MOI.set(m, MOI.RawParameter("SolutionLimit"), 1)
    MOI.set(m, MOI.RawParameter("Heuristics"), 0)
    MOI.set(m, MOI.RawParameter("Presolve"), 0)
    N = 100
    x = MOI.add_variables(m, N)
    for xi in x
        MOI.add_constraint(m, MOI.SingleVariable(xi), MOI.ZeroOne())
        MOI.set(m, MOI.VariablePrimalStart(), xi, 0.0)
    end
    # Given a collection of items with individual weights and values,
    # maximize the total value carried subject to the constraint that
    # the total weight carried is less than 10.
    Random.seed!(1)
    item_weights = rand(N)
    item_values = rand(N)
    MOI.add_constraint(m,
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(item_weights, x), 0.0),
        MOI.LessThan(10.0))
    MOI.set(m, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(.-item_values, x), 0.0))
    MOI.optimize!(m)

    @test MOI.get(m, MOI.TerminationStatus()) == MOI.SOLUTION_LIMIT
    # We should have a primal feasible solution:
    @test MOI.get(m, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
    # But we have no dual status:
    @test MOI.get(m, MOI.DualStatus()) == MOI.NO_SOLUTION
end

function test_Constant_objective_issue_111()
    m = Gurobi.Optimizer(GRB_ENV)
    x = MOI.add_variable(m)
    MOI.set(m, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm{Float64}[], 2.0))
    @test MOI.get(m, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}()).constant == 2.0
    p = Ref{Cdouble}()
    Gurobi.GRBgetdblattr(m, "ObjCon", p)
    @test p[] == 2.0

    MOI.modify(m, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarConstantChange(3.0))
    @test MOI.get(m, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}()).constant == 3.0
    Gurobi.GRBgetdblattr(m.inner, "ObjCon", p)
    @test p[] == 3.0
end

function test_user_provided_env()
    model_1 = Gurobi.Optimizer(GRB_ENV)
    @test model_1.env === GRB_ENV
    model_2 = Gurobi.Optimizer(GRB_ENV)
    @test model_2.env === GRB_ENV
    # Check that finalizer doesn't touch GRB_ENV when manually provided.
    finalize(model_1)
    @test GRB_ENV.ptr_env != C_NULL
end

function test_MULTI_ENV()
    # Gurobi tests should pass if the function begins with test_MULTI_ENV and
    # this specific error is thrown.
    error("Gurobi Error 10009: Failed to obtain a valid license")
end

function test_MULTI_ENV_automatic_env()
    model_1 = Gurobi.Optimizer()
    model_2 = Gurobi.Optimizer()
    @test model_1.env !== model_2.env
    # Check that env is finalized with model when not supplied manually.
    finalize(model_1)
    @test model_1.env.ptr_env == C_NULL
end

function test_user_provided_env_empty()
    model = Gurobi.Optimizer(GRB_ENV)
    @test model.env === GRB_ENV
    @test GRB_ENV.ptr_env != C_NULL
    MOI.empty!(model)
    @test model.env === GRB_ENV
    @test GRB_ENV.ptr_env != C_NULL
end

function test_MULTI_ENV_automatic_env_empty()
    model = Gurobi.Optimizer()
    env = model.env
    MOI.empty!(model)
    @test model.env === env
    @test env.ptr_env != C_NULL
end

function test_MULTI_ENV_manual_finalize()
    env = Gurobi.Env()
    model = Gurobi.Optimizer(env)
    finalize(env)
    @test env.finalize_called
    finalize(model)
    @test env.ptr_env == C_NULL
end

function test_Conflict_refiner_bound_bound()
    model = Gurobi.Optimizer(GRB_ENV)
    MOI.set(model, MOI.Silent(), true)
    x = MOI.add_variable(model)
    c1 = MOI.add_constraint(model, MOI.SingleVariable(x), MOI.GreaterThan(2.0))
    c2 = MOI.add_constraint(model, MOI.SingleVariable(x), MOI.LessThan(1.0))

    # Getting the results before the conflict refiner has been called must return an error.
    @test MOI.get(model, Gurobi.ConflictStatus()) == -1
    @test MOI.get(model, MOI.ConflictStatus()) == MOI.COMPUTE_CONFLICT_NOT_CALLED
    @test_throws ErrorException MOI.get(model, MOI.ConstraintConflictStatus(), c1)

    # Once it's called, no problem.
    MOI.compute_conflict!(model)
    @test MOI.get(model, Gurobi.ConflictStatus()) == 0
    @test MOI.get(model, MOI.ConflictStatus()) == MOI.CONFLICT_FOUND
    @test MOI.get(model, MOI.ConstraintConflictStatus(), c1) == MOI.IN_CONFLICT
    @test MOI.get(model, MOI.ConstraintConflictStatus(), c2) == MOI.IN_CONFLICT
end

function test_Conflict_refiner_bound_affine()
    model = Gurobi.Optimizer(GRB_ENV)
    MOI.set(model, MOI.Silent(), true)
    x = MOI.add_variable(model)
    c1 = MOI.add_constraint(model, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0], [x]), 0.0), MOI.GreaterThan(2.0))
    c2 = MOI.add_constraint(model, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0], [x]), 0.0), MOI.LessThan(1.0))

    # Getting the results before the conflict refiner has been called must return an error.
    @test MOI.get(model, Gurobi.ConflictStatus()) == -1
    @test MOI.get(model, MOI.ConflictStatus()) == MOI.COMPUTE_CONFLICT_NOT_CALLED
    @test_throws ErrorException MOI.get(model, MOI.ConstraintConflictStatus(), c1)

    # Once it's called, no problem.
    MOI.compute_conflict!(model)
    @test MOI.get(model, Gurobi.ConflictStatus()) == 0
    @test MOI.get(model, MOI.ConflictStatus()) == MOI.CONFLICT_FOUND
    @test MOI.get(model, MOI.ConstraintConflictStatus(), c1) == MOI.IN_CONFLICT
    @test MOI.get(model, MOI.ConstraintConflictStatus(), c2) == MOI.IN_CONFLICT
end

function test_Conflict_refiner_invalid_interval()
    model = Gurobi.Optimizer(GRB_ENV)
    MOI.set(model, MOI.Silent(), true)
    x = MOI.add_variable(model)
    c1 = MOI.add_constraint(
        model, MOI.SingleVariable(x), MOI.Interval(1.0, 0.0)
    )
    # Getting the results before the conflict refiner has been called must return an error.
    @test MOI.get(model, Gurobi.ConflictStatus()) == -1
    @test MOI.get(model, MOI.ConflictStatus()) == MOI.COMPUTE_CONFLICT_NOT_CALLED
    @test_throws ErrorException MOI.get(model, MOI.ConstraintConflictStatus(), c1)

    # Once it's called, no problem.
    MOI.compute_conflict!(model)
    @test MOI.get(model, Gurobi.ConflictStatus()) == 0
    @test MOI.get(model, MOI.ConflictStatus()) == MOI.CONFLICT_FOUND
    @test MOI.get(model, MOI.ConstraintConflictStatus(), c1) == MOI.IN_CONFLICT
end

function test_Conflict_refiner_affine_affine()
    model = Gurobi.Optimizer(GRB_ENV)
    MOI.set(model, MOI.Silent(), true)
    x = MOI.add_variable(model)
    y = MOI.add_variable(model)
    b1 = MOI.add_constraint(model, MOI.SingleVariable(x), MOI.GreaterThan(0.0))
    b2 = MOI.add_constraint(model, MOI.SingleVariable(y), MOI.GreaterThan(0.0))
    cf1 = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 1.0], [x, y]), 0.0)
    c1 = MOI.add_constraint(model, cf1, MOI.LessThan(-1.0))
    cf2 = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, -1.0], [x, y]), 0.0)
    c2 = MOI.add_constraint(model, cf2, MOI.GreaterThan(1.0))

    # Getting the results before the conflict refiner has been called must return an error.
    @test MOI.get(model, Gurobi.ConflictStatus()) == -1
    @test MOI.get(model, MOI.ConflictStatus()) == MOI.COMPUTE_CONFLICT_NOT_CALLED
    @test_throws ErrorException MOI.get(model, MOI.ConstraintConflictStatus(), c1)

    # Once it's called, no problem.
    MOI.compute_conflict!(model)
    @test MOI.get(model, Gurobi.ConflictStatus()) == 0
    @test MOI.get(model, MOI.ConflictStatus()) == MOI.CONFLICT_FOUND
    @test MOI.get(model, MOI.ConstraintConflictStatus(), b1) == MOI.IN_CONFLICT
    @test MOI.get(model, MOI.ConstraintConflictStatus(), b2) == MOI.IN_CONFLICT
    @test MOI.get(model, MOI.ConstraintConflictStatus(), c1) == MOI.IN_CONFLICT
    @test MOI.get(model, MOI.ConstraintConflictStatus(), c2) == MOI.NOT_IN_CONFLICT
end

function test_Conflict_refiner_equalto()
    model = Gurobi.Optimizer(GRB_ENV)
    MOI.set(model, MOI.Silent(), true)
    x = MOI.add_variable(model)
    y = MOI.add_variable(model)
    b1 = MOI.add_constraint(model, MOI.SingleVariable(x), MOI.GreaterThan(0.0))
    b2 = MOI.add_constraint(model, MOI.SingleVariable(y), MOI.GreaterThan(0.0))
    cf1 = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 1.0], [x, y]), 0.0)
    c1 = MOI.add_constraint(model, cf1, MOI.EqualTo(-1.0))
    cf2 = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, -1.0], [x, y]), 0.0)
    c2 = MOI.add_constraint(model, cf2, MOI.GreaterThan(1.0))

    # Getting the results before the conflict refiner has been called must return an error.
    @test MOI.get(model, Gurobi.ConflictStatus()) == -1
    @test MOI.get(model, MOI.ConflictStatus()) == MOI.COMPUTE_CONFLICT_NOT_CALLED
    @test_throws ErrorException MOI.get(model, MOI.ConstraintConflictStatus(), c1)

    # Once it's called, no problem.
    MOI.compute_conflict!(model)
    @test MOI.get(model, Gurobi.ConflictStatus()) == 0
    @test MOI.get(model, MOI.ConflictStatus()) == MOI.CONFLICT_FOUND
    @test MOI.get(model, MOI.ConstraintConflictStatus(), b1) == MOI.IN_CONFLICT
    @test MOI.get(model, MOI.ConstraintConflictStatus(), b2) == MOI.IN_CONFLICT
    @test MOI.get(model, MOI.ConstraintConflictStatus(), c1) == MOI.IN_CONFLICT
    @test MOI.get(model, MOI.ConstraintConflictStatus(), c2) == MOI.NOT_IN_CONFLICT
end

function test_Conflict_refiner_outside_conflict()
    model = Gurobi.Optimizer(GRB_ENV)
    MOI.set(model, MOI.Silent(), true)
    x = MOI.add_variable(model)
    y = MOI.add_variable(model)
    z = MOI.add_variable(model)
    b1 = MOI.add_constraint(model, MOI.SingleVariable(x), MOI.GreaterThan(0.0))
    b2 = MOI.add_constraint(model, MOI.SingleVariable(y), MOI.GreaterThan(0.0))
    b3 = MOI.add_constraint(model, MOI.SingleVariable(z), MOI.GreaterThan(0.0))
    cf1 = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 1.0], [x, y]), 0.0)
    c1 = MOI.add_constraint(model, cf1, MOI.LessThan(-1.0))
    cf2 = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, -1.0, 1.0], [x, y, z]), 0.0)
    c2 = MOI.add_constraint(model, cf2, MOI.GreaterThan(1.0))

    # Getting the results before the conflict refiner has been called must return an error.
    @test MOI.get(model, Gurobi.ConflictStatus()) == -1
    @test MOI.get(model, MOI.ConflictStatus()) == MOI.COMPUTE_CONFLICT_NOT_CALLED
    @test_throws ErrorException MOI.get(model, MOI.ConstraintConflictStatus(), c1)

    # Once it's called, no problem.
    MOI.compute_conflict!(model)
    @test MOI.get(model, Gurobi.ConflictStatus()) == 0
    @test MOI.get(model, MOI.ConflictStatus()) == MOI.CONFLICT_FOUND
    @test MOI.get(model, MOI.ConstraintConflictStatus(), b1) == MOI.IN_CONFLICT
    @test MOI.get(model, MOI.ConstraintConflictStatus(), b2) == MOI.IN_CONFLICT
    @test MOI.get(model, MOI.ConstraintConflictStatus(), b3) == MOI.NOT_IN_CONFLICT
    @test MOI.get(model, MOI.ConstraintConflictStatus(), c1) == MOI.IN_CONFLICT
    @test MOI.get(model, MOI.ConstraintConflictStatus(), c2) == MOI.NOT_IN_CONFLICT
end

function test_Conflict_refiner_no_conflict()
    model = Gurobi.Optimizer(GRB_ENV)
    MOI.set(model, MOI.Silent(), true)
    x = MOI.add_variable(model)
    c1 = MOI.add_constraint(model, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0], [x]), 0.0), MOI.GreaterThan(1.0))
    c2 = MOI.add_constraint(model, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0], [x]), 0.0), MOI.LessThan(2.0))

    # Getting the results before the conflict refiner has been called must return an error.
    @test MOI.get(model, Gurobi.ConflictStatus()) == -1
    @test MOI.get(model, MOI.ConflictStatus()) == MOI.COMPUTE_CONFLICT_NOT_CALLED
    @test_throws ErrorException MOI.get(model, MOI.ConstraintConflictStatus(), c1)

    # TODO(odow): bypass Gurobi's IIS checker when underlying model is
    # feasible.
    # @test_throws Gurobi.GurobiError MOI.compute_conflict!(model)
    # @test MOI.get(model, Gurobi.ConflictStatus()) == Gurobi.IIS_NOT_INFEASIBLE
    # @test MOI.get(model, MOI.ConflictStatus()) == MOI.NO_CONFLICT_EXISTS
    # @test MOI.get(model, MOI.ConstraintConflictStatus(), c1) == MOI.NOT_IN_CONFLICT
    # @test MOI.get(model, MOI.ConstraintConflictStatus(), c2) == MOI.NOT_IN_CONFLICT
end

function test_RawParameter()
    model = Gurobi.Optimizer(GRB_ENV)
    @test MOI.get(model, MOI.RawParameter("OutputFlag")) == 1
    MOI.set(model, MOI.RawParameter("OutputFlag"), 0)
    @test MOI.get(model, MOI.RawParameter("OutputFlag")) == 0
end

function test_QCPDual_1()
    model = Gurobi.Optimizer(GRB_ENV)
    MOI.set(model, MOI.Silent(), true)
    MOI.set(model, MOI.RawParameter("QCPDual"), 1)
    MOI.Utilities.loadfromstring!(model, """
    variables: x, y, z
    minobjective: 1.0 * x + 1.0 * y + 1.0 * z
    c1: x + y == 2.0
    c2: x + y + z >= 0.0
    c3: 1.0 * x * x + -1.0 * y * y + -1.0 * z * z >= 0.0
    c4: x >= 0.0
    c5: y >= 0.0
    c6: z >= 0.0
    """)
    MOI.optimize!(model)
    @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL
    @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
    @test MOI.get(model, MOI.DualStatus()) == MOI.FEASIBLE_POINT
    c1 = MOI.get(model, MOI.ConstraintIndex, "c1")
    c2 = MOI.get(model, MOI.ConstraintIndex, "c2")
    c3 = MOI.get(model, MOI.ConstraintIndex, "c3")
    @test MOI.get(model, MOI.ConstraintDual(), c1) ≈ 1.0 atol=1e-6
    @test MOI.get(model, MOI.ConstraintDual(), c2) ≈ 0.0 atol=1e-6
    @test MOI.get(model, MOI.ConstraintDual(), c3) ≈ 0.0 atol=1e-6
end

function test_QCPDual_default()
    model = Gurobi.Optimizer(GRB_ENV)
    MOI.set(model, MOI.Silent(), true)
    MOI.Utilities.loadfromstring!(model, """
    variables: x, y, z
    minobjective: 1.0 * x + 1.0 * y + 1.0 * z
    c1: x + y == 2.0
    c2: x + y + z >= 0.0
    c3: 1.0 * x * x + -1.0 * y * y + -1.0 * z * z >= 0.0
    c4: x >= 0.0
    c5: y >= 0.0
    c6: z >= 0.0
    """)
    MOI.optimize!(model)
    @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL
    @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
    @test MOI.get(model, MOI.DualStatus()) == MOI.NO_SOLUTION
    c1 = MOI.get(model, MOI.ConstraintIndex, "c1")
    c2 = MOI.get(model, MOI.ConstraintIndex, "c2")
    c3 = MOI.get(model, MOI.ConstraintIndex, "c3")
    @test_throws ErrorException MOI.get(model, MOI.ConstraintDual(), c1)
    @test_throws ErrorException MOI.get(model, MOI.ConstraintDual(), c2)
    @test_throws ErrorException MOI.get(model, MOI.ConstraintDual(), c3)
end

function test_Add_and_delete_constraints()
    model = Gurobi.Optimizer(GRB_ENV)
    x = MOI.add_variables(model, 2)
    cs = MOI.add_constraints(
        model,
        [MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x[i])], 0.0) for i in 1:2],
        MOI.EqualTo.([0.0, 0.0])
    )
    @test MOI.get(model, MOI.NumberOfConstraints{
        MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}
    }()) == 2
    MOI.delete(model, cs)
    @test iszero(MOI.get(model, MOI.NumberOfConstraints{
        MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}
    }()))
end

function test_Buffered_deletion_test()
    # Check if the VarHintVal of variables (could be any attribute that is not
    # cached, the most common as: the bounds, name, and the start, are all
    # cached) are right if we remove some variables and add new ones without
    # updating the model (and that they stay correct after updating the model,
    # here done by `MOI.optimize!`, as the update needs to be done inside MOI
    # code, that is aware of the lazy updates).
    model = Gurobi.Optimizer(GRB_ENV)
    MOI.set(model, MOI.Silent(), true)
    vars = MOI.add_variables(model, 3)
    vars_obj = [7.0, 42.0, -0.5]
    obj_attr = Gurobi.VariableAttribute("Obj")
    MOI.set.(model, obj_attr, vars, vars_obj)
    @test all(MOI.is_valid.(model, vars))
    @test MOI.get.(model, obj_attr, vars) == vars_obj
    MOI.delete(model, vars[[1, 3]])
    @test !MOI.is_valid(model, vars[1])
    @test !MOI.is_valid(model, vars[3])
    fourth_var = MOI.add_variable(model)
    fourth_var_obj = -77.0
    MOI.set(model, obj_attr, fourth_var, fourth_var_obj)
    # Check before updating the model and after the delete+insert.
    @test !MOI.is_valid(model, vars[1])
    @test !MOI.is_valid(model, vars[3])
    @test MOI.is_valid(model, vars[2])
    @test vars_obj[2] == MOI.get(model, obj_attr, vars[2])
    @test MOI.is_valid(model, fourth_var)
    @test fourth_var_obj == MOI.get(model, obj_attr, fourth_var)
    # Then optimize to force an update.
    MOI.optimize!(model)
    # And check again.
    @test MOI.is_valid(model, vars[2])
    @test vars_obj[2] == MOI.get(model, obj_attr, vars[2])
    @test MOI.is_valid(model, fourth_var)
    @test fourth_var_obj == MOI.get(model, obj_attr, fourth_var)
end

function test_extra_name_Variables()
    model = Gurobi.Optimizer(GRB_ENV)
    x = MOI.add_variables(model, 3)
    MOI.set(model, MOI.VariableName(), x[1], "x1")
    @test MOI.get(model, MOI.VariableIndex, "x1") == x[1]
    MOI.set(model, MOI.VariableName(), x[1], "x2")
    @test MOI.get(model, MOI.VariableIndex, "x1") === nothing
    @test MOI.get(model, MOI.VariableIndex, "x2") == x[1]
    MOI.set(model, MOI.VariableName(), x[2], "x1")
    @test MOI.get(model, MOI.VariableIndex, "x1") == x[2]
    MOI.set(model, MOI.VariableName(), x[3], "xα")
    @test MOI.get(model, MOI.VariableIndex, "xα") == x[3]
    MOI.set(model, MOI.VariableName(), x[1], "x1")
    @test_throws ErrorException MOI.get(model, MOI.VariableIndex, "x1")
end

function test_extra_name_Variable_bounds()
    model = Gurobi.Optimizer(GRB_ENV)
    x = MOI.add_variable(model)
    c1 = MOI.add_constraint(model, MOI.SingleVariable(x), MOI.GreaterThan(0.0))
    c2 = MOI.add_constraint(model, MOI.SingleVariable(x), MOI.LessThan(1.0))
    MOI.set(model, MOI.ConstraintName(), c1, "c1")
    @test MOI.get(model, MOI.ConstraintIndex, "c1") == c1
    MOI.set(model, MOI.ConstraintName(), c1, "c2")
    @test MOI.get(model, MOI.ConstraintIndex, "c1") === nothing
    @test MOI.get(model, MOI.ConstraintIndex, "c2") == c1
    MOI.set(model, MOI.ConstraintName(), c2, "c1")
    @test MOI.get(model, MOI.ConstraintIndex, "c1") == c2
    MOI.set(model, MOI.ConstraintName(), c1, "c1")
    @test_throws ErrorException MOI.get(model, MOI.ConstraintIndex, "c1")
end

function test_extra_name_Affine_constraints()
    model = Gurobi.Optimizer(GRB_ENV)
    x = MOI.add_variable(model)
    f = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x)], 0.0)
    c1 = MOI.add_constraint(model, f, MOI.GreaterThan(0.0))
    c2 = MOI.add_constraint(model, f, MOI.LessThan(1.0))
    MOI.set(model, MOI.ConstraintName(), c1, "c1")
    @test MOI.get(model, MOI.ConstraintIndex, "c1") == c1
    MOI.set(model, MOI.ConstraintName(), c1, "c2")
    @test MOI.get(model, MOI.ConstraintIndex, "c1") === nothing
    @test MOI.get(model, MOI.ConstraintIndex, "c2") == c1
    MOI.set(model, MOI.ConstraintName(), c2, "c1")
    @test MOI.get(model, MOI.ConstraintIndex, "c1") == c2
    MOI.set(model, MOI.ConstraintName(), c1, "c1")
    @test_throws ErrorException MOI.get(model, MOI.ConstraintIndex, "c1")
end

function test_ConstraintAttribute()
    model = Gurobi.Optimizer(GRB_ENV)
    MOI.Utilities.loadfromstring!(model, """
variables: x
minobjective: x
c1: x >= 0.0
c2: 2x >= 1.0
c3: x in Integer()
""")
    c2 = MOI.get(model, MOI.ConstraintIndex, "c2")
    c3 = MOI.get(model, MOI.ConstraintIndex, "c3")
    # Linear constraints are supported - one test for each different type.
    # Integer attribute
    MOI.set(model, Gurobi.ConstraintAttribute("Lazy"), c2, 2)
    @test MOI.get(model, Gurobi.ConstraintAttribute("Lazy"), c2) == 2
    # Real attribute
    MOI.set(model, Gurobi.ConstraintAttribute("RHS"), c2, 2.0)
    @test MOI.get(model, Gurobi.ConstraintAttribute("RHS"), c2) == 2.0
    # Char attribute
    MOI.set(model, Gurobi.ConstraintAttribute("Sense"), c2, '<')
    @test MOI.get(model, Gurobi.ConstraintAttribute("Sense"), c2) == '<'
    # String attribute
    MOI.set(model, Gurobi.ConstraintAttribute("ConstrName"), c2, "c4")
    @test MOI.get(model, Gurobi.ConstraintAttribute("ConstrName"), c2) == "c4"
    # Things that should fail follow.
    # Non-linear constraints are not supported.
    @test_throws(
        MOI.SetAttributeNotAllowed(Gurobi.ConstraintAttribute("Lazy")),
        MOI.set(model, Gurobi.ConstraintAttribute("Lazy"), c3, 1)
    )
    # Getting/setting a non-existing attribute.
    attr = Gurobi.ConstraintAttribute("Non-existing")
    err = ErrorException("Gurobi Error 10004: Unknown attribute 'Non-existing'")
    @test_throws err MOI.set(model, attr, c2, 1)
    @test_throws err MOI.get(model, attr, c2)
    # Setting an attribute to a value of the wrong type.
    @test_throws(
        ArgumentError("Attribute Lazy requires Int64 arguments. Provided argument was of type Float64."),
        MOI.set(model, Gurobi.ConstraintAttribute("Lazy"), c2, 1.5)
    )
end

function test_VariableAttribute()
    model = Gurobi.Optimizer(GRB_ENV)
    MOI.Utilities.loadfromstring!(model, """
variables: x
minobjective: x
c1: x >= 0.0
c2: 2x >= 1.0
c3: x in Integer()
""")
    x = MOI.get(model, MOI.VariableIndex, "x")
    # Setting attributes of each type
    # Integer attribute
    MOI.set(model, Gurobi.VariableAttribute("VarHintPri"), x, 2)
    @test MOI.get(model, Gurobi.VariableAttribute("VarHintPri"), x) == 2
    # Real Attribute
    MOI.set(model, Gurobi.VariableAttribute("LB"), x, 2.0)
    @test MOI.get(model, Gurobi.VariableAttribute("LB"), x) == 2.0
    # Char Attribute
    MOI.set(model, Gurobi.VariableAttribute("VType"), x, 'B')
    @test MOI.get(model, Gurobi.VariableAttribute("VType"), x) == 'B'
    # String Attribute
    MOI.set(model, Gurobi.VariableAttribute("VarName"), x, "my_var")
    @test MOI.get(model, Gurobi.VariableAttribute("VarName"), x) == "my_var"
    # Things that should fail follow.
    # Getting/setting a non-existing attribute.
    attr = Gurobi.VariableAttribute("Non-existing")
    err = ErrorException("Gurobi Error 10004: Unknown attribute 'Non-existing'")
    @test_throws err MOI.set(model, attr, x, 1)
    @test_throws err MOI.get(model, attr, x)
    # Setting an attribute to a value of the wrong type.
    @test_throws(
        ArgumentError("Attribute BranchPriority requires Int64 arguments. Provided argument was of type Float64."),
        MOI.set(model, Gurobi.VariableAttribute("BranchPriority"), x, 1.5)
    )
end

function test_ModelAttribute()
    model = Gurobi.Optimizer(GRB_ENV)
    MOI.Utilities.loadfromstring!(model, """
variables: x
minobjective: x
c1: x >= 0.0
c2: 2x >= 1.0
c3: x in Integer()
""")
    # Setting attributes of each type
    # Integer attribute
    MOI.set(model, Gurobi.ModelAttribute("ModelSense"), -1)
    @test MOI.get(model, Gurobi.ModelAttribute("ModelSense")) == -1
    # Real Attribute
    MOI.set(model, Gurobi.ModelAttribute("ObjCon"), 3.0)
    @test MOI.get(model, Gurobi.ModelAttribute("ObjCon")) == 3.0
    # String Attribute
    MOI.set(model, Gurobi.ModelAttribute("ModelName"), "My model")
    @test MOI.get(model, Gurobi.ModelAttribute("ModelName")) == "My model"
    # Things that should fail follow.
    # Getting/setting a non-existing attribute.
    attr = Gurobi.ModelAttribute("Non-existing")
    err = ErrorException("Gurobi Error 10004: Unknown attribute 'Non-existing'")
    @test_throws err MOI.set(model, attr, 1)
    @test_throws err MOI.get(model, attr)
    # Setting an attribute to a value of the wrong type.
    @test_throws(
        ArgumentError("Attribute NumStart requires $(Int) arguments. Provided argument was of type Float64."),
        MOI.set(model, Gurobi.ModelAttribute("NumStart"), 4.0)
    )
end

function test_soc_no_initial_bound()
    model = Gurobi.Optimizer(GRB_ENV)
    MOI.set(model, MOI.Silent(), true)
    t = MOI.add_variable(model)
    x = MOI.add_variables(model, 2)
    MOI.add_constraints(
        model, MOI.SingleVariable.(x), MOI.GreaterThan.(3.0:4.0)
    )
    c_soc = MOI.add_constraint(
        model, MOI.VectorOfVariables([t; x]), MOI.SecondOrderCone(3)
    )
    MOI.set(model, MOI.ObjectiveFunction{MOI.SingleVariable}(), MOI.SingleVariable(t))
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.optimize!(model)
    @test MOI.get(model, MOI.VariablePrimal(), t) == 5.0
    MOI.delete(model, c_soc)
    MOI.optimize!(model)
    @test MOI.get(model, MOI.TerminationStatus()) == MOI.DUAL_INFEASIBLE
end

function test_soc_nonnegative_initial_bound()
    model = Gurobi.Optimizer(GRB_ENV)
    MOI.set(model, MOI.Silent(), true)
    t = MOI.add_variable(model)
    x = MOI.add_variables(model, 2)
    MOI.add_constraint(model, MOI.SingleVariable(t), MOI.GreaterThan(1.0))
    MOI.add_constraints(
        model, MOI.SingleVariable.(x), MOI.GreaterThan.(3.0:4.0)
    )
    c_soc = MOI.add_constraint(
        model, MOI.VectorOfVariables([t; x]), MOI.SecondOrderCone(3)
    )
    MOI.set(model, MOI.ObjectiveFunction{MOI.SingleVariable}(), MOI.SingleVariable(t))
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.optimize!(model)
    @test MOI.get(model, MOI.VariablePrimal(), t) == 5.0
    MOI.delete(model, c_soc)
    MOI.optimize!(model)
    @test MOI.get(model, MOI.VariablePrimal(), t) == 1.0
end

function test_soc_negative_initial_bound()
    model = Gurobi.Optimizer(GRB_ENV)
    MOI.set(model, MOI.Silent(), true)
    t = MOI.add_variable(model)
    x = MOI.add_variables(model, 2)
    MOI.add_constraint(model, MOI.SingleVariable(t), MOI.GreaterThan(-1.0))
    MOI.add_constraints(
        model, MOI.SingleVariable.(x), MOI.GreaterThan.(3.0:4.0)
    )
    c_soc = MOI.add_constraint(
        model, MOI.VectorOfVariables([t; x]), MOI.SecondOrderCone(3)
    )
    MOI.optimize!(model)
    MOI.set(model, MOI.ObjectiveFunction{MOI.SingleVariable}(), MOI.SingleVariable(t))
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    @test MOI.get(model, MOI.VariablePrimal(), t) == 5.0
    MOI.delete(model, c_soc)
    MOI.optimize!(model)
    @test MOI.get(model, MOI.VariablePrimal(), t) == -1.0
end

function test_soc_nonnegative_post_bound()
    model = Gurobi.Optimizer(GRB_ENV)
    MOI.set(model, MOI.Silent(), true)
    t = MOI.add_variable(model)
    x = MOI.add_variables(model, 2)
    MOI.add_constraints(
        model, MOI.SingleVariable.(x), MOI.GreaterThan.(3.0:4.0)
    )
    MOI.add_constraint(
        model, MOI.VectorOfVariables([t; x]), MOI.SecondOrderCone(3)
    )
    c_lb = MOI.add_constraint(model, MOI.SingleVariable(t), MOI.GreaterThan(6.0))
    MOI.set(model, MOI.ObjectiveFunction{MOI.SingleVariable}(), MOI.SingleVariable(t))
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.optimize!(model)
    @test MOI.get(model, MOI.VariablePrimal(), t) == 6.0
    MOI.delete(model, c_lb)
    MOI.optimize!(model)
    @test MOI.get(model, MOI.VariablePrimal(), t) == 5.0
end

function test_soc_negative_post_bound()
    model = Gurobi.Optimizer(GRB_ENV)
    MOI.set(model, MOI.Silent(), true)
    t = MOI.add_variable(model)
    x = MOI.add_variables(model, 2)
    MOI.add_constraints(
        model, MOI.SingleVariable.(x), MOI.GreaterThan.(3.0:4.0)
    )
    c_soc = MOI.add_constraint(
        model, MOI.VectorOfVariables([t; x]), MOI.SecondOrderCone(3)
    )
    MOI.add_constraint(model, MOI.SingleVariable(t), MOI.GreaterThan(-6.0))
    MOI.set(model, MOI.ObjectiveFunction{MOI.SingleVariable}(), MOI.SingleVariable(t))
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.optimize!(model)
    @test MOI.get(model, MOI.VariablePrimal(), t) == 5.0
    MOI.delete(model, c_soc)
    MOI.optimize!(model)
    @test MOI.get(model, MOI.VariablePrimal(), t) == -6.0
end

function test_soc_negative_post_bound_ii()
    model = Gurobi.Optimizer(GRB_ENV)
    MOI.set(model, MOI.Silent(), true)
    t = MOI.add_variable(model)
    x = MOI.add_variables(model, 2)
    MOI.add_constraints(
        model, MOI.SingleVariable.(x), MOI.GreaterThan.(3.0:4.0)
    )
    c_soc = MOI.add_constraint(
        model, MOI.VectorOfVariables([t; x]), MOI.SecondOrderCone(3)
    )
    c_lb = MOI.add_constraint(model, MOI.SingleVariable(t), MOI.GreaterThan(-6.0))
    MOI.set(model, MOI.ObjectiveFunction{MOI.SingleVariable}(), MOI.SingleVariable(t))
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.optimize!(model)
    @test MOI.get(model, MOI.VariablePrimal(), t) == 5.0
    MOI.delete(model, c_lb)
    MOI.optimize!(model)
    @test MOI.get(model, MOI.VariablePrimal(), t) == 5.0
    MOI.delete(model, c_soc)
    MOI.optimize!(model)
    @test MOI.get(model, MOI.TerminationStatus()) == MOI.DUAL_INFEASIBLE
end

function test_soc_negative_post_bound_iii()
    # This test was added because add_constraint and add_constraints had
    # different implementations and add_constraints failed where
    # add_constraint succeeded.
    model = Gurobi.Optimizer(GRB_ENV)
    MOI.set(model, MOI.Silent(), true)
    t = MOI.add_variable(model)
    x = MOI.add_variables(model, 2)
    c_soc = MOI.add_constraint(
        model, MOI.VectorOfVariables([t; x]), MOI.SecondOrderCone(3)
    )
    c_lbs = MOI.add_constraints(
        model, MOI.SingleVariable.([t; x]), MOI.GreaterThan.([-6.0, 3.0, 4.0])
    )
    c_lb = first(c_lbs)
    MOI.set(model, MOI.ObjectiveFunction{MOI.SingleVariable}(), MOI.SingleVariable(t))
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.optimize!(model)
    @test MOI.get(model, MOI.VariablePrimal(), t) == 5.0
    MOI.delete(model, c_lb)
    MOI.optimize!(model)
    @test MOI.get(model, MOI.VariablePrimal(), t) == 5.0
    MOI.delete(model, c_soc)
    MOI.optimize!(model)
    @test MOI.get(model, MOI.TerminationStatus()) == MOI.DUAL_INFEASIBLE
end

function test_2_soc()
    model = Gurobi.Optimizer(GRB_ENV)
    MOI.set(model, MOI.Silent(), true)
    t = MOI.add_variable(model)
    x = MOI.add_variables(model, 2)
    MOI.add_constraint(model, MOI.SingleVariable(t), MOI.GreaterThan(-1.0))
    MOI.add_constraints(
        model, MOI.SingleVariable.(x), MOI.GreaterThan.([4.0, 3.0])
    )
    c_soc_1 = MOI.add_constraint(
        model, MOI.VectorOfVariables([t, x[1]]), MOI.SecondOrderCone(2)
    )
    c_soc_2 = MOI.add_constraint(
        model, MOI.VectorOfVariables([t, x[2]]), MOI.SecondOrderCone(2)
    )
    MOI.optimize!(model)
    MOI.set(model, MOI.ObjectiveFunction{MOI.SingleVariable}(), MOI.SingleVariable(t))
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    @test MOI.get(model, MOI.VariablePrimal(), t) == 4.0
    MOI.delete(model, c_soc_1)
    MOI.optimize!(model)
    @test MOI.get(model, MOI.VariablePrimal(), t) == 3.0
    MOI.delete(model, c_soc_2)
    MOI.optimize!(model)
    @test MOI.get(model, MOI.VariablePrimal(), t) == -1.0
end

function test_Duplicate_names_Variables()
    model = Gurobi.Optimizer(GRB_ENV)
    (x, y, z) = MOI.add_variables(model, 3)
    MOI.set(model, MOI.VariableName(), x, "x")
    MOI.set(model, MOI.VariableName(), y, "x")
    MOI.set(model, MOI.VariableName(), z, "z")
    @test MOI.get(model, MOI.VariableIndex, "z") == z
    @test_throws ErrorException MOI.get(model, MOI.VariableIndex, "x")
    MOI.set(model, MOI.VariableName(), y, "y")
    @test MOI.get(model, MOI.VariableIndex, "x") == x
    @test MOI.get(model, MOI.VariableIndex, "y") == y
    MOI.set(model, MOI.VariableName(), z, "x")
    @test_throws ErrorException MOI.get(model, MOI.VariableIndex, "x")
    MOI.delete(model, x)
    @test MOI.get(model, MOI.VariableIndex, "x") == z
end

function test_Duplicate_names_SingleVariable()
    model = Gurobi.Optimizer(GRB_ENV)
    x = MOI.add_variables(model, 3)
    c = MOI.add_constraints(model, MOI.SingleVariable.(x), MOI.GreaterThan(0.0))
    MOI.set(model, MOI.ConstraintName(), c[1], "x")
    MOI.set(model, MOI.ConstraintName(), c[2], "x")
    MOI.set(model, MOI.ConstraintName(), c[3], "z")
    @test MOI.get(model, MOI.ConstraintIndex, "z") == c[3]
    @test_throws ErrorException MOI.get(model, MOI.ConstraintIndex, "x")
    MOI.set(model, MOI.ConstraintName(), c[2], "y")
    @test MOI.get(model, MOI.ConstraintIndex, "x") == c[1]
    @test MOI.get(model, MOI.ConstraintIndex, "y") == c[2]
    MOI.set(model, MOI.ConstraintName(), c[3], "x")
    @test_throws ErrorException MOI.get(model, MOI.ConstraintIndex, "x")
    MOI.delete(model, c[1])
    @test MOI.get(model, MOI.ConstraintIndex, "x") == c[3]
    MOI.set(model, MOI.ConstraintName(), c[2], "x")
    @test_throws ErrorException MOI.get(model, MOI.ConstraintIndex, "x")
    MOI.delete(model, x[3])
    @test MOI.get(model, MOI.ConstraintIndex, "x") == c[2]
end

function test_Duplicate_names_ScalarAffineFunction()
    model = Gurobi.Optimizer(GRB_ENV)
    x = MOI.add_variables(model, 3)
    fs = [
        MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, xi)], 0.0)
        for xi in x
    ]
    c = MOI.add_constraints(model, fs, MOI.GreaterThan(0.0))
    MOI.set(model, MOI.ConstraintName(), c[1], "x")
    MOI.set(model, MOI.ConstraintName(), c[2], "x")
    MOI.set(model, MOI.ConstraintName(), c[3], "z")
    @test MOI.get(model, MOI.ConstraintIndex, "z") == c[3]
    @test_throws ErrorException MOI.get(model, MOI.ConstraintIndex, "x")
    MOI.set(model, MOI.ConstraintName(), c[2], "y")
    @test MOI.get(model, MOI.ConstraintIndex, "x") == c[1]
    @test MOI.get(model, MOI.ConstraintIndex, "y") == c[2]
    MOI.set(model, MOI.ConstraintName(), c[3], "x")
    @test_throws ErrorException MOI.get(model, MOI.ConstraintIndex, "x")
    MOI.delete(model, c[1])
    @test MOI.get(model, MOI.ConstraintIndex, "x") == c[3]
end

function test_Duals_with_equal_bounds_250()
    model = Gurobi.Optimizer(GRB_ENV)
    MOI.set(model, MOI.Silent(), true)
    x = MOI.add_variable(model)
    xl = MOI.add_constraint(model, MOI.SingleVariable(x), MOI.GreaterThan(1.0))
    xu = MOI.add_constraint(model, MOI.SingleVariable(x), MOI.LessThan(1.0))
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.set(model, MOI.ObjectiveFunction{MOI.SingleVariable}(), MOI.SingleVariable(x))
    MOI.optimize!(model)
    @test MOI.get(model, MOI.ConstraintDual(), xl) == 1.0
    @test MOI.get(model, MOI.ConstraintDual(), xu) == 0.0
end

function test_Objective_functions()
    model = Gurobi.Optimizer(GRB_ENV)
    x = MOI.add_variable(model)
    @test MOI.get(model, MOI.ObjectiveSense()) == MOI.FEASIBILITY_SENSE
    @test MOI.get(model, MOI.ListOfModelAttributesSet()) == Any[MOI.ObjectiveSense()]
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.set(model, MOI.ObjectiveFunction{MOI.SingleVariable}(), MOI.SingleVariable(x))
    @test MOI.get(model, MOI.ListOfModelAttributesSet()) ==
        Any[MOI.ObjectiveSense(), MOI.ObjectiveFunction{MOI.SingleVariable}()]
    MOI.set(model, MOI.ObjectiveSense(), MOI.FEASIBILITY_SENSE)
    @test MOI.get(model, MOI.ListOfModelAttributesSet()) == Any[MOI.ObjectiveSense()]
end

function test_FEASIBILITY_SENSE_zeros_objective()
    model = Gurobi.Optimizer(GRB_ENV)
    MOI.set(model, MOI.Silent(), true)
    x = MOI.add_variable(model)
    MOI.add_constraint(model, MOI.SingleVariable(x), MOI.GreaterThan(1.0))
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.set(model, MOI.ObjectiveFunction{MOI.SingleVariable}(), MOI.SingleVariable(x))
    MOI.optimize!(model)
    @test MOI.get(model, MOI.ObjectiveValue()) == 1.0
    MOI.set(model, MOI.ObjectiveSense(), MOI.FEASIBILITY_SENSE)
    MOI.optimize!(model)
    @test MOI.get(model, MOI.ObjectiveValue()) == 0.0
end

function test_Attributes()
    model = Gurobi.Optimizer(GRB_ENV)
    MOI.set(model, MOI.Silent(), true)
    MOI.optimize!(model)
    @test MOI.get(model, MOI.NodeCount()) == 0
    @test MOI.get(model, MOI.BarrierIterations()) == 0
    @test MOI.get(model, MOI.SimplexIterations()) == 0
end

function test_GRBterminate()
    model = Gurobi.Optimizer(GRB_ENV)
    MOI.set(model, MOI.Silent(), true)
    x = MOI.add_variable(model)
    MOI.set(model, Gurobi.CallbackFunction(), (cb_data, cb_where) -> begin
        GRBterminate(model)
    end)
    MOI.optimize!(model)
    @test MOI.get(model, MOI.TerminationStatus()) == MOI.INTERRUPTED
end

"""
    test_InterruptException()

This test simulates an InterruptException being thrown. It is a little
complicated due to the delayed handling of GRBterminate, which _schedules_ a
request for termination, rather than terminating immediately. This means Gurobi
may continue to call the callback after the interruption.

First, we must ensure that InterruptException() is only thrown once. Double
interrupting would interrupt our handling of the first interrupt!

Second, if the model is too simplisitic, Gurobi may be able to prove optimality
after we have interrupted, but before it has decided to actually exit the solve.
"""
function test_InterruptException()
    model = Gurobi.Optimizer(GRB_ENV)
    MOI.set(model, MOI.Silent(), true)
    x = MOI.add_variable(model)
    MOI.add_constraint(model, MOI.SingleVariable(x), MOI.Integer())
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.SingleVariable}(),
        MOI.SingleVariable(x),
    )
    MOI.set(model, MOI.RawParameter("LazyConstraints"), 1)
    interrupt_thrown = false
    i = 0.0
    MOI.set(model, Gurobi.CallbackFunction(), (cb_data, cb_where) -> begin
        if cb_where != Gurobi.GRB_CB_MIPSOL
            return
        end
        MOI.submit(
            model,
            MOI.LazyConstraint(cb_data),
            MOI.ScalarAffineFunction{Float64}(
                [MOI.ScalarAffineTerm(1.0, x)], 0.0
            ),
            MOI.GreaterThan{Float64}(i)
        )
        i += 1
        if !interrupt_thrown
            interrupt_thrown = true
            throw(InterruptException())
        end
    end)
    MOI.optimize!(model)
    @test MOI.get(model, MOI.TerminationStatus()) == MOI.INTERRUPTED
end

function test_indicator_name()
    model = Gurobi.Optimizer(GRB_ENV)
    MOI.set(model, MOI.Silent(), true)
    x = MOI.add_variables(model, 2)
    MOI.add_constraint(model, MOI.SingleVariable(x[1]), MOI.ZeroOne())
    f = MOI.VectorAffineFunction(
        [
            MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, x[1])),
            MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(2.0, x[2])),
        ],
        [0.0, 0.0],
    )
    s = MOI.IndicatorSet{MOI.ACTIVATE_ON_ONE}(MOI.GreaterThan(1.0))
    c = MOI.add_constraint(model, f, s)
    MOI.set(model, MOI.ConstraintName(), c, "my_indicator")
    @test MOI.get(model, MOI.ConstraintName(), c) == "my_indicator"
end

function test_indicator_on_one()
    model = Gurobi.Optimizer(GRB_ENV)
    MOI.set(model, MOI.Silent(), true)
    x = MOI.add_variables(model, 2)
    MOI.add_constraint(model, MOI.SingleVariable(x[1]), MOI.ZeroOne())
    f = MOI.VectorAffineFunction(
        [
            MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, x[1])),
            MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(2.0, x[2])),
        ],
        [0.0, 0.0],
    )
    s = MOI.IndicatorSet{MOI.ACTIVATE_ON_ONE}(MOI.GreaterThan(1.0))
    c = MOI.add_constraint(model, f, s)
    @test MOI.get(model, MOI.ConstraintSet(), c) == s
    @test isapprox(MOI.get(model, MOI.ConstraintFunction(), c), f)
end

function test_indicator_on_zero()
    model = Gurobi.Optimizer(GRB_ENV)
    MOI.set(model, MOI.Silent(), true)
    x = MOI.add_variables(model, 2)
    MOI.add_constraint(model, MOI.SingleVariable(x[1]), MOI.ZeroOne())
    f = MOI.VectorAffineFunction(
        [
            MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, x[1])),
            MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(2.0, x[2])),
        ],
        [0.0, 0.0],
    )
    s = MOI.IndicatorSet{MOI.ACTIVATE_ON_ZERO}(MOI.GreaterThan(1.0))
    c = MOI.add_constraint(model, f, s)
    @test MOI.get(model, MOI.ConstraintSet(), c) == s
    @test isapprox(MOI.get(model, MOI.ConstraintFunction(), c), f)
end

function test_indicator_nonconstant_x()
    model = Gurobi.Optimizer(GRB_ENV)
    MOI.set(model, MOI.Silent(), true)
    x = MOI.add_variables(model, 2)
    MOI.add_constraint(model, MOI.SingleVariable(x[1]), MOI.ZeroOne())
    f = MOI.VectorAffineFunction(
        [
            MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(-1.0, x[1])),
            MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(2.0, x[2])),
        ],
        [0.0, 0.0],
    )
    s = MOI.IndicatorSet{MOI.ACTIVATE_ON_ONE}(MOI.GreaterThan(1.0))
    @test_throws ErrorException MOI.add_constraint(model, f, s)
end

function test_indicator_too_many_indicators()
    model = Gurobi.Optimizer(GRB_ENV)
    MOI.set(model, MOI.Silent(), true)
    x = MOI.add_variables(model, 2)
    MOI.add_constraint(model, MOI.SingleVariable(x[1]), MOI.ZeroOne())
    f = MOI.VectorAffineFunction(
        [
            MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(-1.0, x[1])),
            MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(2.0, x[2])),
        ],
        [0.0, 0.0],
    )
    s = MOI.IndicatorSet{MOI.ACTIVATE_ON_ONE}(MOI.GreaterThan(1.0))
    @test_throws ErrorException MOI.add_constraint(model, f, s)
end

function test_indicator_nonconstant()
    model = Gurobi.Optimizer(GRB_ENV)
    MOI.set(model, MOI.Silent(), true)
    x = MOI.add_variables(model, 2)
    MOI.add_constraint(model, MOI.SingleVariable(x[1]), MOI.ZeroOne())
    f = MOI.VectorAffineFunction(
        [
            MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, x[1])),
            MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(2.0, x[2])),
        ],
        [1.0, 0.0],
    )
    s = MOI.IndicatorSet{MOI.ACTIVATE_ON_ONE}(MOI.GreaterThan(1.0))
    @test_throws ErrorException MOI.add_constraint(model, f, s)
end

function test_farkas_dual_min()
    model = Gurobi.Optimizer(GRB_ENV)
    MOI.set(model, MOI.Silent(), true)
    MOI.set(model, MOI.RawParameter("InfUnbdInfo"), 1)
    x = MOI.add_variables(model, 2)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.SingleVariable}(),
        MOI.SingleVariable(x[1]),
    )
    clb = MOI.add_constraint.(
        model, MOI.SingleVariable.(x), MOI.GreaterThan(0.0)
    )
    c = MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([2.0, 1.0], x), 0.0),
        MOI.LessThan(-1.0),
    )
    MOI.optimize!(model)
    @test MOI.get(model, MOI.TerminationStatus()) == MOI.INFEASIBLE
    @test MOI.get(model, MOI.DualStatus()) == MOI.INFEASIBILITY_CERTIFICATE
    clb_dual = MOI.get.(model, MOI.ConstraintDual(), clb)
    c_dual = MOI.get(model, MOI.ConstraintDual(), c)
    @show clb_dual, c_dual
    @test clb_dual[1] > -1e-6
    @test clb_dual[2] > -1e-6
    @test c_dual[1] < 1e-6
    @test clb_dual[1] ≈ -2 * c_dual atol = 1e-6
    @test clb_dual[2] ≈ -c_dual atol = 1e-6
end

function test_farkas_dual_min_interval()
    model = Gurobi.Optimizer(GRB_ENV)
    MOI.set(model, MOI.Silent(), true)
    MOI.set(model, MOI.RawParameter("InfUnbdInfo"), 1)
    x = MOI.add_variables(model, 2)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.SingleVariable}(),
        MOI.SingleVariable(x[1]),
    )
    clb = MOI.add_constraint.(
        model, MOI.SingleVariable.(x), MOI.Interval(0.0, 10.0)
    )
    c = MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([2.0, 1.0], x), 0.0),
        MOI.LessThan(-1.0),
    )
    MOI.optimize!(model)
    @test MOI.get(model, MOI.TerminationStatus()) == MOI.INFEASIBLE
    @test MOI.get(model, MOI.DualStatus()) == MOI.INFEASIBILITY_CERTIFICATE
    clb_dual = MOI.get.(model, MOI.ConstraintDual(), clb)
    c_dual = MOI.get(model, MOI.ConstraintDual(), c)
    @show clb_dual, c_dual
    @test clb_dual[1] > -1e-6
    @test clb_dual[2] > -1e-6
    @test c_dual[1] < 1e-6
    @test clb_dual[1] ≈ -2 * c_dual atol = 1e-6
    @test clb_dual[2] ≈ -c_dual atol = 1e-6
end

function test_farkas_dual_min_equalto()
    model = Gurobi.Optimizer(GRB_ENV)
    MOI.set(model, MOI.Silent(), true)
    MOI.set(model, MOI.RawParameter("InfUnbdInfo"), 1)
    x = MOI.add_variables(model, 2)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.SingleVariable}(),
        MOI.SingleVariable(x[1]),
    )
    clb = MOI.add_constraint.(model, MOI.SingleVariable.(x), MOI.EqualTo(0.0))
    c = MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([2.0, 1.0], x), 0.0),
        MOI.LessThan(-1.0),
    )
    MOI.optimize!(model)
    @test MOI.get(model, MOI.TerminationStatus()) == MOI.INFEASIBLE
    @test MOI.get(model, MOI.DualStatus()) == MOI.INFEASIBILITY_CERTIFICATE
    clb_dual = MOI.get.(model, MOI.ConstraintDual(), clb)
    c_dual = MOI.get(model, MOI.ConstraintDual(), c)
    @show clb_dual, c_dual
    @test clb_dual[1] > -1e-6
    @test clb_dual[2] > -1e-6
    @test c_dual[1] < 1e-6
    @test clb_dual[1] ≈ -2 * c_dual atol = 1e-6
    @test clb_dual[2] ≈ -c_dual atol = 1e-6
end

function test_farkas_dual_min_ii()
    model = Gurobi.Optimizer(GRB_ENV)
    MOI.set(model, MOI.Silent(), true)
    MOI.set(model, MOI.RawParameter("InfUnbdInfo"), 1)
    x = MOI.add_variables(model, 2)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(-1.0, x[1])], 0.0),
    )
    clb = MOI.add_constraint.(
        model, MOI.SingleVariable.(x), MOI.LessThan(0.0)
    )
    c = MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([-2.0, -1.0], x), 0.0),
        MOI.LessThan(-1.0),
    )
    MOI.optimize!(model)
    @test MOI.get(model, MOI.TerminationStatus()) == MOI.INFEASIBLE
    @test MOI.get(model, MOI.DualStatus()) == MOI.INFEASIBILITY_CERTIFICATE
    clb_dual = MOI.get.(model, MOI.ConstraintDual(), clb)
    c_dual = MOI.get(model, MOI.ConstraintDual(), c)
    @show clb_dual, c_dual
    @test clb_dual[1] < 1e-6
    @test clb_dual[2] < 1e-6
    @test c_dual[1] < 1e-6
    @test clb_dual[1] ≈ 2 * c_dual atol = 1e-6
    @test clb_dual[2] ≈ c_dual atol = 1e-6
end

function test_farkas_dual_max()
    model = Gurobi.Optimizer(GRB_ENV)
    MOI.set(model, MOI.Silent(), true)
    MOI.set(model, MOI.RawParameter("InfUnbdInfo"), 1)
    x = MOI.add_variables(model, 2)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.SingleVariable}(),
        MOI.SingleVariable(x[1]),
    )
    clb = MOI.add_constraint.(
        model, MOI.SingleVariable.(x), MOI.GreaterThan(0.0)
    )
    c = MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([2.0, 1.0], x), 0.0),
        MOI.LessThan(-1.0),
    )
    MOI.optimize!(model)
    @test MOI.get(model, MOI.TerminationStatus()) == MOI.INFEASIBLE
    @test MOI.get(model, MOI.DualStatus()) == MOI.INFEASIBILITY_CERTIFICATE
    clb_dual = MOI.get.(model, MOI.ConstraintDual(), clb)
    c_dual = MOI.get(model, MOI.ConstraintDual(), c)
    @show clb_dual, c_dual
    @test clb_dual[1] > -1e-6
    @test clb_dual[2] > -1e-6
    @test c_dual[1] < 1e-6
    @test clb_dual[1] ≈ -2 * c_dual atol = 1e-6
    @test clb_dual[2] ≈ -c_dual atol = 1e-6
end

function test_farkas_dual_max_ii()
    model = Gurobi.Optimizer(GRB_ENV)
    MOI.set(model, MOI.Silent(), true)
    MOI.set(model, MOI.RawParameter("InfUnbdInfo"), 1)
    x = MOI.add_variables(model, 2)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(-1.0, x[1])], 0.0),
    )
    clb = MOI.add_constraint.(
        model, MOI.SingleVariable.(x), MOI.LessThan(0.0)
    )
    c = MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([-2.0, -1.0], x), 0.0),
        MOI.LessThan(-1.0),
    )
    MOI.optimize!(model)
    @test MOI.get(model, MOI.TerminationStatus()) == MOI.INFEASIBLE
    @test MOI.get(model, MOI.DualStatus()) == MOI.INFEASIBILITY_CERTIFICATE
    clb_dual = MOI.get.(model, MOI.ConstraintDual(), clb)
    c_dual = MOI.get(model, MOI.ConstraintDual(), c)
    @show clb_dual, c_dual
    @test clb_dual[1] < 1e-6
    @test clb_dual[2] < 1e-6
    # Confirmed bug in Gurobi <=9.1 and fixed in next release. When updating for new
    # release, confirm fixed and add a conditional check on `_GRB_VERSION`.
    @test_broken c_dual[1] < 1e-6
    @test_broken clb_dual[1] ≈ 2 * c_dual atol = 1e-6
    @test_broken clb_dual[2] ≈ c_dual atol = 1e-6
end

function _build_basis_model()
    T = Float64
    model = Gurobi.Optimizer(GRB_ENV)
    MOI.set(model, MOI.Silent(), true)
    # Min -x
    # s.t. x + y <= 1
    # x, y >= 0

    x = MOI.add_variable(model)
    y = MOI.add_variable(model)

    cf = MOI.ScalarAffineFunction{T}(
        MOI.ScalarAffineTerm{T}.([one(T), one(T)], [x, y]),
        zero(T),
    )
    c = MOI.add_constraint(model, cf, MOI.LessThan(one(T)))

    vc1 = MOI.add_constraint(
        model,
        MOI.SingleVariable(x),
        MOI.GreaterThan(zero(T)),
    )
    vc2 = MOI.add_constraint(
        model,
        MOI.SingleVariable(y),
        MOI.GreaterThan(zero(T)),
    )
    objf = MOI.ScalarAffineFunction{T}(
        MOI.ScalarAffineTerm{T}.([-one(T), zero(T)], [x, y]),
        zero(T),
    )
    MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}}(), objf)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    return model, vc1, vc2, c
end

function test_set_basis()
    # The following is an indirect, brittle way to test the desired behavior.
    # Gurobi appears to not allow you to query VBasis/CBasis attributes without
    # optimizing the problem
    # (https://support.gurobi.com/hc/en-us/community/posts/360075729911-Getting-VBasis-CBasis-attributes-without-optimizing).
    # The following just verifies that the problem is solved with 0 simplex iterations
    # when we seed the optimal solution and 1 simplex iteration when we seed a
    # suboptimal feasible basis. Ideally, we would remove the call to MOI.optimize!
    # and just verify that, after setting the basis, we can get the same basis
    # statuses back.
    let
        model, vc1, vc2, c = _build_basis_model()
        MOI.set(model, MOI.ConstraintBasisStatus(), vc1, MOI.BASIC)
        MOI.set(model, MOI.ConstraintBasisStatus(), vc2, MOI.NONBASIC)
        MOI.set(model, MOI.ConstraintBasisStatus(), c, MOI.NONBASIC)
        MOI.optimize!(model)
        @test MOI.get(model, MOI.SimplexIterations()) == 0
    end

    let
        model, vc1, vc2, c = _build_basis_model()
        MOI.set(model, MOI.ConstraintBasisStatus(), vc1, MOI.NONBASIC)
        MOI.set(model, MOI.ConstraintBasisStatus(), vc2, MOI.NONBASIC)
        MOI.set(model, MOI.ConstraintBasisStatus(), c, MOI.BASIC)
        MOI.optimize!(model)
        @test MOI.get(model, MOI.SimplexIterations()) == 1
    end
end

end

runtests(TestMOIWrapper)


text:
using Test
using PolynomialBasis
using ImplicitDomainQuadrature
# using Revise
using CutCell
include("useful_routines.jl")

polyorder = 1
numqp = 2

basis = TensorProductBasis(2, polyorder)
levelset = InterpolatingPolynomial(1, basis)
nf = CutCell.number_of_basis_functions(basis)
mesh = CutCell.Mesh([0.0, 0.0], [3.0, 1.0], [3, 1], nf)
nodalcoordinates = CutCell.nodal_coordinates(mesh)
nodalconnectivity = CutCell.nodal_connectivity(mesh)

normal = [1.0, 0.0]
x0 = [0.5, 0.0]
tol = perturbation = 1e-3
levelsetcoeffs = plane_distance_function(nodalcoordinates, normal, x0)

cellsign = CutCell.cell_sign!(levelset, levelsetcoeffs, nodalconnectivity, tol, perturbation)
cellmap = CutCell.cell_map(mesh,1)

interfacequads = CutCell.InterfaceQuadratures(
    cellsign,
    levelset,
    levelsetcoeffs,
    nodalconnectivity,
    numqp,
    cellmap,
)
@test size(interfacequads.quads) == (2,1)
testcelltoquad = [1, 0, 0]
@test allequal(interfacequads.celltoquad, testcelltoquad)

testnormals = repeat(normal,inner=(1,numqp))
@test allapprox(testnormals,CutCell.interface_normals(interfacequads,1))

@test interfacequads[1,1] ≈ interfacequads[-1,1]


text:
println("Test metrics")

import Distributions.SpecialFunctions.logabsbeta
import Distributions.StatsBase.harmmean
lbeta(a, b) = logabsbeta(a, b)[1]

@testset "metrics" begin
  @testset "dic" begin
    x = randn(10000)
    ll = normlogpdf.(x)
    @test isfinite(dic(ll))
  end

  @testset "logharmmean" begin
    Random.seed!(0)
    x = randn(100)
    @test log(harmmean(exp.(x))) ≈ MCMC.logharmmeanexp(x) atol=1e-6
  end

  @testset "log_bayes_factor" begin
    Random.seed!(0)
    n = 100
    y = rand(Bernoulli(.7), n)
    sumy = sum(y)

    m0 = Beta(70, 100)
    m1 = Beta(100, 70)

    compute_logc_true(m) = lbeta(m.α + sumy, m.β + n - sumy) - lbeta(m.α, m.β)
    c0_true = compute_logc_true(m0)
    c1_true = compute_logc_true(m1)

    true_log_bf = c1_true - c0_true

    post_samps_m0 = rand(Beta(m0.α + sumy, m0.β + n - sumy), 100_000)
    post_samps_m1 = rand(Beta(m1.α + sumy, m1.β + n - sumy), 100_000)
    ll0 = map(p -> sum(logpdf.(Bernoulli(p), y)), post_samps_m0)
    ll1 = map(p -> sum(logpdf.(Bernoulli(p), y)), post_samps_m1)

    approx_log_bf = log_bayes_factor(ll0, ll1)

    @test c0_true ≈ MCMC.logharmmeanexp(ll0)  atol=0.2
    @test c1_true ≈ MCMC.logharmmeanexp(ll1)  atol=0.2

    @test true_log_bf ≈ approx_log_bf atol=0.2
  end

  @testset "logmeanexp" begin
    Random.seed!(0)
    x = randn(10000)
    @test logmeanexp(x) ≈ log(mean(exp.(x)))
  end

  @testset "hellinger" begin
    Random.seed!(0)
    d1 = Normal(0, 1)
    d2 = Normal(1, 1)
    @test hellinger(d1, d2, Int(1e6)) ≈ hellinger(d2, d1, Int(1e6))  atol=0.02
    @test hellinger(d1, d2, Int(1e6)) ≈ 0.34 atol=0.01
  end
end


text:
module SurfaceRegression
import ..AbstractOperator.get_misfit
using ..AbstractOperator, ..CommonToAll
using PyPlot, LinearAlgebra, StatsBase

import ..Model, ..Options, ..OptionsStat, ..OptionsNonstat

export Surface, SurfaceWithDifferentData

mutable struct Surface <: Operator2D
    d      :: Array{Float64}
    useML  :: Bool
    σ      :: Array{Float64}
    r      :: Array{Float64}
    select :: Array{Bool}
end

function Surface(d::Array{Float64} ;useML=false, σ=1.0)
    if length(d) != length(σ)
        @assert length(σ) == 1
        σ = σ[1]*ones(length(d))
    end
    r = d[.!isnan.(d)]
    select = .!isnan.(d)
    Surface(d, useML, σ, r, select)
end

function get_misfit(m::Model, opt::Options, surface::Surface)
    chi2by2 = 0.0
    if !opt.debug
        d, σ, select, r  = surface.d, surface.σ, surface.select, surface.r
        r .= (m.fstar[vec(select)] .- d[select])./σ[select]
        chi2 = r'*r
        if surface.useML
            n = length(r)
            chi2by2 = 0.5*n*log(chi2)
        else
            chi2by2 = 0.5*chi2
        end
    end
    return chi2by2
end

mutable struct SurfaceWithDifferentData <: Operator2D
    surfaces :: Array{Surface, 1}
end

function get_misfit(m::Model, opt::Options, surface::SurfaceWithDifferentData)
    chi2by2 = 0.0
    for s in surface.surfaces
        chi2by2 += get_misfit(m, opt, s)
    end
    return chi2by2
end

function relativeuncertainty(opt::Options, surface::SurfaceWithDifferentData; burninfrac=0.5, temperaturenum=1)
    M = assembleTat1(opt, :fstar, burninfrac=burninfrac, temperaturenum=temperaturenum)
    Σ = [[s.σ[s.select] for m in M] for s in surface.surfaces]
    scaling = [zeros(length(M)) for s in surface.surfaces]
    for (im, m) in enumerate(M)
        for (is, s) in enumerate(surface.surfaces)
            d, σ, select, r  = s.d, s.σ, s.select, s.r
            r .= (vec(m[vec(select)]) .- d[select])./σ[select]
            n = length(r)
            scaling[is][im] = sqrt(r'r/n)
            Σ[is][im] .= scaling[is][im]*s.σ[s.select]
        end
    end
    finalΣ = [vcat(s...) for s in Σ]
    scaling, finalΣ
end

function slice_image_posterior( M::AbstractArray, opt::Options, roworcol::Symbol;
                                rowcolnum = 1,
                                nbins = 50,
                                rhomin=Inf,
                                rhomax=-Inf,
                                qp1=0.05,
                                qp2=0.95,
                                pdfnormalize=false,
                                temperaturenum=1)
    @assert temperaturenum == 1
    Mslices = Array{Array{Float64}, 1}(undef, length(M))
    x = unique(opt.xall[1,:])
    y = unique(opt.xall[2,:])
    for i in 1:length(M)
        m = reshape(M[i],length(y), length(x))
        if roworcol == :row
            Mslices[i] = m[rowcolnum,:]
        else
            Mslices[i] = m[:,rowcolnum]
        end
    end
    himage, edges, CI, = gethimage(Mslices, opt, burninfrac=0.0, temperaturenum=temperaturenum,
                nbins=nbins, rhomin=rhomin, rhomax=rhomax, qp1=qp1, qp2=qp2,
                pdfnormalize=pdfnormalize)
end

function get_image_quantile(M::AbstractArray, q=0.5)
    n = length(M[1])
    inrows = size(M[1], 2) == 1
    if inrows
        mall = hcat(M...)
    else
        mall = vcat(M...)
    end
    quantM = zeros(n)
    for i = 1:n
        quantM[i] = inrows ? quantile(mall[i,:], q) : quantile(mall[:,i], q)
    end
    quantM
end

function plot_surface_posterior(optns::OptionsNonstat, opts::OptionsStat;
                        rownum = 10,
                        colnum = 10,
                        temperaturenum = 1,
                        nbins = 50,
                        burninfrac=0.5,
                        qp1=0.05,
                        qp2=0.95,
                        cmappdf = "magma",
                        cmapmean = "bone",
                        figsizerows=[11,6],
                        figsizecols=[7.5,9.5],
                        pdfnormalize=false,
                        fsize=10,
                        domean = false,
                        getquant = 0.5,
                        concisefigsize=(10,5),
                        property="pixel value",
                        property_units="")
    @assert temperaturenum == 1
    @assert 0.0 < getquant < 1.0
    M = assembleTat1(optns, :fstar, burninfrac=burninfrac, temperaturenum=temperaturenum)
    if domean
        mns = mean(M)
    else
        mns = get_image_quantile(M, getquant)
    end
    x = vec(unique(optns.xall[1,:]))
    delx = diff(x)[1]
    xmesh = x[1]-delx/2:delx:x[end]+delx/2
    y = vec(unique(optns.xall[2,:]))
    dely = diff(y)[1]
    ymesh = y[1]-dely/2:dely:y[end]+dely/2
    mns = reshape(mns, length(y), length(x))
    CIimage = reshape(getCIimage(M, qp1, qp2), size(mns))
    himage_r_ns, edges_r_ns, CI_r_ns = slice_image_posterior(M, optns, :row, rowcolnum=rownum, nbins = nbins, qp1=qp1, qp2=qp2,
                                pdfnormalize=pdfnormalize, temperaturenum=temperaturenum)
    himage_c_ns, edges_c_ns, CI_c_ns = slice_image_posterior(M, optns, :col, rowcolnum=colnum, nbins = nbins, qp1=qp1, qp2=qp2,
                                pdfnormalize=pdfnormalize, temperaturenum=temperaturenum)
    M = assembleTat1(opts, :fstar, burninfrac=burninfrac, temperaturenum=temperaturenum)
    if domean
        ms = 0.5*log10.(mean(M))
    else
        a, b = [m[1,:] for m in M], [m[2,:] for m in M]
        ms = 0.5*log10.([get_image_quantile(a, getquant) get_image_quantile(b, getquant)])'
    end
    msx = reshape(ms[1,:], length(y), length(x))
    msy = reshape(ms[2,:], length(y), length(x))
    Mx, My = Array{Array{Float64}, 1}(undef, length(M)), Array{Array{Float64}, 1}(undef, length(M))
    for i = 1:length(M)
        Mx[i] = M[i][1,:]
        My[i] = M[i][2,:]
    end
    himage_x_c, edges_x_c, CI_x_c = slice_image_posterior(Mx, opts, :col, rowcolnum=colnum, nbins = nbins, qp1=qp1, qp2=qp2,
                                pdfnormalize=pdfnormalize, temperaturenum=temperaturenum)
    himage_y_c, edges_y_c, CI_y_c = slice_image_posterior(My, opts, :col, rowcolnum=colnum, nbins = nbins, qp1=qp1, qp2=qp2,
                                pdfnormalize=pdfnormalize, temperaturenum=temperaturenum)
    himage_x_r, edges_x_r, CI_x_r = slice_image_posterior(Mx, opts, :row, rowcolnum=rownum, nbins = nbins, qp1=qp1, qp2=qp2,
                                pdfnormalize=pdfnormalize, temperaturenum=temperaturenum)
    himage_y_r, edges_y_r, CI_y_r = slice_image_posterior(My, opts, :row, rowcolnum=rownum, nbins = nbins, qp1=qp1, qp2=qp2,
                                pdfnormalize=pdfnormalize, temperaturenum=temperaturenum)
    f = figure(figsize=figsizerows)
    s1 = subplot(231)
    im1 = s1.imshow(mns, extent=[xmesh[1], xmesh[end], ymesh[1], ymesh[end]], origin="lower", cmap=cmapmean)
    s1.plot([x[1], x[end]], [y[rownum], y[rownum]], "--k", alpha=0.5)
    central_tendency = domean ? ("Mean") : ("Percentile "*"$(round(Int, getquant*100))")
    s1.set_title(central_tendency*" "*property)
    s1.set_ylabel("distance y")
    s1.set_xlabel("distance x")
    cb1 = colorbar(im1, ax=s1)
    s2 = subplot(232, sharex=s1, sharey=s1)
    im2 = s2.imshow(msx, extent=[xmesh[1], xmesh[end], ymesh[1], ymesh[end]], origin="lower", cmap=cmapmean)
    s2.plot([x[1], x[end]], [y[rownum], y[rownum]], "--k", alpha=0.5)
    s2.set_title(central_tendency*" "*L"\log_{10}λ_x")
    s2.set_ylabel("distance y")
    s2.set_xlabel("distance x")
    cb2 = colorbar(im2, ax=s2)
    s3 = subplot(233, sharex=s1, sharey=s1)
    im3 = s3.imshow(msy, extent=[xmesh[1], xmesh[end], ymesh[1], ymesh[end]], origin="lower", cmap=cmapmean)
    s3.plot([x[1], x[end]], [y[rownum], y[rownum]], "--k", alpha=0.5)
    s3.set_title(central_tendency*" "*L"\log_{10}λ_y")
    s3.set_ylabel("distance y")
    s3.set_xlabel("distance x")
    cb3 = colorbar(im3, ax=s3)
    s4 = subplot(234, sharex=s1)
    im4 = s4.pcolormesh(xmesh[:], edges_r_ns[:], himage_r_ns', cmap=cmappdf)
    s4.plot(x, mns[rownum,:][:],"--w", alpha=0.5)
    s4.plot(x, CI_r_ns,"--c", alpha=0.5)
    s4.set_title(property*" PDF")
    s4.set_ylabel(property)
    s4.set_xlabel("distance x")
    cb4 = colorbar(im4, ax=s4)
    cb4.ax.set_title("pdf")
    s5 = subplot(235, sharex=s1)
    im5 = s5.pcolormesh(xmesh[:], edges_x_r[:], himage_x_r', cmap=cmappdf)
    cb5 = colorbar(im5, ax=s5)
    cb5.ax.set_title("pdf")
    s5.plot(x, msx[rownum,:][:],"--w", alpha=0.5)
    s5.plot(x, CI_x_r,"--c", alpha=0.5)
    s5.set_title(L"\log_{10}λ_x"*" PDF")
    s5.set_ylabel(L"\log_{10}λ_x")
    s5.set_xlabel("distance x")
    s6 = subplot(236, sharex=s1)
    im6 = s6.pcolormesh(xmesh[:], edges_y_r[:], himage_y_r', cmap=cmappdf)
    cb6 = colorbar(im6, ax=s6)
    cb6.ax.set_title("pdf")
    s6.plot(x, msy[rownum,:][:],"--w", alpha=0.5)
    s6.plot(x, CI_y_r,"--c", alpha=0.5)
    s6.set_title(L"\log_{10}λ_y"*" PDF")
    s6.set_ylabel(L"\log_{10}λ_y")
    s6.set_xlabel("distance x")
    nicenup(f, fsize=fsize)

    f = figure(figsize=figsizecols)
    s1 = subplot(321)
    im1 = s1.imshow(mns, extent=[xmesh[1], xmesh[end], ymesh[1], ymesh[end]], origin="lower", cmap=cmapmean)
    s1.plot([y[colnum], y[colnum]], [y[1], y[end]], "--k", alpha=0.5)
    s1.set_title(central_tendency*" "*property)
    s1.set_ylabel("distance y")
    s1.set_xlabel("distance x")
    # s1.set_aspect(aspect)
    cb1 = colorbar(im1, ax=s1)
    s2 = subplot(323, sharex=s1, sharey=s1)
    im2 = s2.imshow(msx, extent=[xmesh[1], xmesh[end], ymesh[1], ymesh[end]], origin="lower", cmap=cmapmean)
    s2.plot([y[colnum], y[colnum]], [y[1], y[end]], "--k", alpha=0.5)
    s2.set_title(central_tendency*" "*L"\log_{10}λ_x")
    s2.set_ylabel("distance y")
    s2.set_xlabel("distance x")
    # s2.set_aspect(aspect)
    s3 = subplot(325, sharex=s1, sharey=s1)
    cb2 = colorbar(im2, ax=s2)
    im3 = s3.imshow(msy, extent=[xmesh[1], xmesh[end], ymesh[1], ymesh[end]], origin="lower", cmap=cmapmean)
    s3.plot([y[colnum], y[colnum]], [y[1], y[end]], "--k", alpha=0.5)
    # s3.set_aspect(aspect)
    s3.set_title(central_tendency*" "*L"\log_{10}λ_y")
    s3.set_ylabel("distance y")
    s3.set_xlabel("distance x")
    cb3 = colorbar(im3, ax=s3)
    s4 = subplot(322, sharey=s1)
    im4 = s4.pcolormesh(edges_c_ns[:], ymesh[:], himage_c_ns, cmap=cmappdf)
    s4.plot(mns[:,colnum][:], y, "--w", alpha=0.5)
    s4.plot(CI_c_ns, y, "--c", alpha=0.5)
    s4.set_title(property*" PDF")
    s4.set_xlabel(property)
    s4.set_ylabel("distance y")
    cb4 = colorbar(im4, ax=s4)
    cb4.ax.set_title("pdf")
    s5 = subplot(324, sharey=s1)
    im5 = s5.pcolormesh(edges_x_c[:], ymesh[:], himage_x_c, cmap=cmappdf)
    cb5 = colorbar(im5, ax=s5)
    cb5.ax.set_title("pdf")
    s5.plot(msx[:,colnum][:], y, "--w", alpha=0.5)
    s5.plot(CI_x_c, y, "--c", alpha=0.5)
    s5.set_title(L"\log_{10}λ_x"*" PDF")
    s5.set_xlabel(L"\log_{10}λ_x")
    s5.set_ylabel("distance y")
    s6 = subplot(326, sharey=s1)
    im6 = s6.pcolormesh(edges_y_c[:], ymesh[:], himage_y_c, cmap=cmappdf)
    cb6 = colorbar(im6, ax=s6)
    cb6.ax.set_title("pdf")
    s6.plot(msy[:,colnum][:], y, "--w", alpha=0.5)
    s6.plot(CI_y_c, y, "--c", alpha=0.5)
    s6.set_title(L"\log_{10}λ_y"*" PDF")
    s6.set_xlabel(L"\log_{10}λ_y")
    s6.set_ylabel("distance y")
    nicenup(f, fsize=fsize)

    f, ax = plt.subplots(1, 2, sharex=true, sharey=true, figsize=concisefigsize)
    im1 = ax[1].imshow(mns, extent=[xmesh[1], xmesh[end], ymesh[1], ymesh[end]], origin="lower", cmap=cmapmean)
    central_tendency = domean ? ("Mean") : ("Percentile "*"$(round(Int, getquant*100))")
    ax[1].set_title(central_tendency*" "*property)
    ax[1].set_ylabel("distance y")
    ax[1].set_xlabel("distance x")
    cb1 = colorbar(im1, ax=ax[1])
    cb1.ax.set_xlabel(property_units)
    im2 = ax[2].imshow(CIimage, extent=[xmesh[1], xmesh[end], ymesh[1], ymesh[end]], origin="lower", cmap=cmappdf)
    ax[2].set_title("Credible interval: $(round(Int, 100*qp2)) - $(round(Int, 100*qp1))")
    ax[2].set_xlabel("distance x")
    cb2 = colorbar(im2, ax=ax[2])
    cb2.ax.set_xlabel(property_units)
end

function getCIimage(marray, qp1, qp2)
    q1 = get_image_quantile(marray, qp1)
    q2 = get_image_quantile(marray, qp2)
    return q2-q1
end

function plot_surface_posterior(opt::Options;
                        rownum = 10,
                        colnum = 10,
                        temperaturenum = 1,
                        nbins = 50,
                        burninfrac=0.5,
                        qp1=0.05,
                        qp2=0.95,
                        cmappdf = "magma",
                        cmapmean = "bone",
                        figsize=[7.5,5.9],
                        pdfnormalize=false,
                        fsize=10,
                        aspect=1.0,
                        domean = false,
                        getquant = 0.5,
                        property = "pixel value",
                        property_units="",
                        concisefigsize=(10,5))
    @assert temperaturenum == 1
    @assert 0.0 < getquant < 1.0
    if typeof(opt) == OptionsStat
        @assert opt.updatenonstat == false
        @assert opt.needλ²fromlog == false
    end
    x = vec(unique(opt.xall[1,:]))
    delx = diff(x)[1]
    xmesh = x[1]-delx/2:delx:x[end]+delx/2
    y = vec(unique(opt.xall[2,:]))
    dely = diff(y)[1]
    ymesh = y[1]-dely/2:dely:y[end]+dely/2
    M = assembleTat1(opt, :fstar, burninfrac=burninfrac, temperaturenum=temperaturenum)
    if domean
        m = mean(M)
    else
        a = [m' for m in M]
        m = get_image_quantile(a, getquant)
    end
    m = reshape(m, length(y), length(x))
    CIimage = reshape(getCIimage(M, qp1, qp2), size(m))
    himage_r, edges_r, CI_r = slice_image_posterior(M, opt, :row, rowcolnum=rownum, nbins = nbins, qp1=qp1, qp2=qp2,
                                pdfnormalize=pdfnormalize, temperaturenum=temperaturenum)
    himage_c, edges_c, CI_c = slice_image_posterior(M, opt, :col, rowcolnum=colnum, nbins = nbins, qp1=qp1, qp2=qp2,
                                pdfnormalize=pdfnormalize, temperaturenum=temperaturenum)
    f = figure(figsize=figsize)
    s1 = subplot(223)
    central_tendency = domean ? ("Mean") : ("Percentile "*"$(round(Int, getquant*100))")
    im1 = s1.imshow(m, extent=[xmesh[1], xmesh[end], ymesh[1], ymesh[end]], origin="lower", cmap=cmapmean)
    s1.plot([x[1], x[end]], [y[rownum], y[rownum]], "--k", alpha=0.5)
    s1.plot([y[colnum], y[colnum]], [y[1], y[end]], "--k", alpha=0.5)
    s1.set_title(central_tendency*" "*property)
    s1.set_ylabel("distance y")
    # s1.set_aspect(aspect)
    s1.set_xlabel("distance x")
    cb1 = colorbar(im1, ax=s1)
    s2 = subplot(221, sharex=s1)
    im2 = s2.pcolormesh(xmesh[:], edges_r[:], himage_r', cmap=cmappdf)
    s2.plot(x, m[rownum,:][:],"--w", alpha=0.5)
    s2.plot(x, CI_r, "--c", alpha=0.5)
    s2.set_title(property*" PDF")
    s2.set_ylabel(property)
    s2.set_xlabel("distance x")
    cb2 = colorbar(im2, ax=s2)
    cb2.ax.set_title("pdf")
    s3 = subplot(224, sharey=s1)
    im3 = s3.pcolormesh(edges_c[:], ymesh[:], himage_c, cmap=cmappdf)
    s3.plot(m[:,colnum][:], y, "--w", alpha=0.5)
    s3.plot(CI_c, y, "--c", alpha=0.5)
    s3.set_title(property*" PDF")
    s3.set_xlabel(property)
    s3.set_ylabel("distance y")
    cb3 = colorbar(im3, ax=s3)
    cb3.ax.set_title("pdf")
    nicenup(f, fsize=fsize)

    f, ax = plt.subplots(1, 2, sharex=true, sharey=true, figsize=concisefigsize)
    im1 = ax[1].imshow(m, extent=[xmesh[1], xmesh[end], ymesh[1], ymesh[end]], origin="lower", cmap=cmapmean)
    central_tendency = domean ? ("Mean") : ("Percentile "*"$(round(Int, getquant*100))")
    ax[1].set_title(central_tendency*" "*property)
    ax[1].set_ylabel("distance y")
    ax[1].set_xlabel("distance x")
    cb1 = colorbar(im1, ax=ax[1])
    cb1.ax.set_xlabel(property_units)
    im2 = ax[2].imshow(CIimage, extent=[xmesh[1], xmesh[end], ymesh[1], ymesh[end]], origin="lower", cmap=cmappdf)
    ax[2].set_title("Credible interval: $(round(Int, 100*qp2)) - $(round(Int, 100*qp1))")
    ax[2].set_xlabel("distance x")
    cb2 = colorbar(im2, ax=ax[2])
    cb2.ax.set_xlabel(property_units)
end

end

max_stars_repo_name: GeoscienceAustralia/HiQGA

max_issues_repo_name: GeoscienceAustralia/HiQGA


text:
# This file is a part of SimilaritySearch.jl
# License is Apache 2.0: https://www.apache.org/licenses/LICENSE-2.0.txt

export Laesa

mutable struct Laesa{T} <: Index
    db::Vector{T}
    pivots::Vector{T}
    table::Matrix{Float64} # rows: number of pivots; cols: number of objects
end

"""
    fit(::Type{Laesa}, dist::Function, db::AbstractVector{T}, pivots::Vector{T})
    fit(::Type{Laesa}, dist::Function, db::AbstractVector{T}, numPivots::Integer)

Creates a `Laesa` index with the given pivots. If the number of pivots is specified,
then they will be randomly selected from the dataset.
"""
function fit(::Type{Laesa}, dist::Function, db::AbstractVector{T}, pivots::AbstractVector{T})  where T
    @info "Creating a pivot table with $(length(pivots)) pivots and distance=$(dist)"
    table = Matrix{Float64}(undef, length(pivots), length(db))

    for j in 1:length(db)
        for i in 1:length(pivots)
            table[i, j] = dist(db[j], pivots[i])
        end
    end

    Laesa(db, pivots, table)
end

function fit(::Type{Laesa}, dist::Function, db::AbstractVector{T}, numPivots::Integer) where T
    pivots = rand(db, numPivots)
    fit(Laesa, dist, db, pivots)
end

"""
    search(index::Laesa, dist::Function, q, res::KnnResult)

Solves a query with the Laesa index.
"""
function search(index::Laesa, dist::Function, q, res::KnnResult)
    dqp = [dist(q, piv) for piv in index.pivots]
    for i in 1:length(index.db)
        dpu = @view index.table[:, i]

        evaluate = true
        for pivID in 1:length(index.pivots)
            if abs(dqp[pivID] - dpu[pivID]) > covrad(res)
                evaluate = false
                break
            end
        end

        if evaluate
            d = dist(q, index.db[i])
            push!(res, i, d)
        end
    end

    return res
end

"""
    push!(index::Laesa, dist::Function, obj)

Inserts `obj` into the index
"""
function push!(index::Laesa, dist::Function, obj)
    push!(index.db, obj)
    vec = Vector{Float64}(undef, length(index.pivots))
    for pivID in 1:length(vec)
        vec[pivID] = dist(index.pivots[pivID], obj)
    end

    index.table = hcat(index.table, vec)
    length(index.db)
end