text:
"""
    seed_g(T::Type,j::Int64,N::Int64)

Creates a `x::SVector{N,T}` object that is one at `x[j]` and zero everywhere else.
"""
function seed_g(T::Type,j::Int64,N::Int64)
    return SVector{N,T}([i == j ? 1.0 : 0.0 for i=1:N])
end

"""
    grad(x::SMCg{N,T},j,n) where {N,T}

sets convex and concave (sub)gradients of length `n` of `x` to be `1` at index `j`
"""
function grad(x::SMCg{N,T},j::Int64) where {N,T<:AbstractFloat}
  sv_grad::SVector{N,T} = seed_g(T,j,N)
  return SMCg{N,T}(x.cc,x.cv,sv_grad,sv_grad,x.Intv,x.cnst,x.IntvBox,x.xref)
end

"""
    zgrad(x::SMCg{N,T},n::Int64) where {N,T}

sets convex and concave (sub)gradients of length `n` to be zero
"""
function zgrad(x::SMCg{N,T}) where {N,T<:AbstractFloat}
  grad::SVector{N,T} = @SVector zeros(N)
  return SMCg{N,T}(x.cc,x.cv,grad,grad,x.Intv,x.cnst,x.IntvBox,x.xref)
end

function convert(::Type{SMCg{N,T}},x::S) where {S<:Integer,N,T<:AbstractFloat}
          seed::SVector{N,T} = @SVector zeros(T,N)
          SMCg{N,T}(convert(T,x),convert(T,x),seed,seed,Interval(convert(Interval{T},x)),
                    false,[emptyinterval(T)],[zero(T)])
end
function convert(::Type{SMCg{N,T}},x::S) where {S<:AbstractFloat,N,T<:AbstractFloat}
          seed::SVector{N,T} = @SVector zeros(T,N)
          SMCg{N,T}(convert(T,x),convert(T,x),seed,seed,Interval(convert(Interval{T},x)),
                    false,[emptyinterval(T)],[zero(T)])
end
function convert(::Type{SMCg{N,T}},x::S) where {S<:Interval,N,T<:AbstractFloat}
          seed::SVector{N,T} = @SVector zeros(T,N)
          SMCg{N,T}(convert(T,x.hi),convert(T,x.lo),seed,seed,convert(Interval{T},x),
                    false,[emptyinterval(T)],[zero(T)])
end

promote_rule(::Type{SMCg{N,T}}, ::Type{S}) where {S<:Integer,N,T<:AbstractFloat} = SMCg{N,T}
promote_rule(::Type{SMCg{N,T}}, ::Type{S}) where {S<:AbstractFloat,N,T<:AbstractFloat} = SMCg{N,T}
promote_rule(::Type{SMCg{N,T}}, ::Type{S}) where {S<:Interval,N,T<:AbstractFloat} = SMCg{N,T}
promote_rule(::Type{SMCg{N,T}}, ::Type{S}) where {S<:Real,N,T<:AbstractFloat} = SMCg{N,T}

"""
    mid3(x::T,y::T,z::T)

Calculates the midpoint of three numbers returning the value and the index.
"""
function mid3(x::T,y::T,z::T) where {T<:AbstractFloat}
  (((x>=y)&&(y>=z))||((z>=y)&&(y>=x))) && (return y,2)
  (((y>=x)&&(x>=z))||((z>=x)&&(x>=y))) && (return x,1)
  return z,3
end

"""
    mid_grad(cc_grad::SVector{N,T}, cv_grad::SVector{N,T}, id::Int64)

Takes the concave relaxation gradient 'cc_grad', the convex relaxation gradient
'cv_grad', and the index of the midpoint returned 'id' and outputs the appropriate
gradient according to McCormick relaxation rules.
"""
function mid_grad(cc_grad::SVector{N,T}, cv_grad::SVector{N,T}, id::Int64) where {N,T<:AbstractFloat}
  if (id == 1)
    return cc_grad
  elseif (id == 2)
    return cv_grad
  elseif (id == 3)
    return zero(cc_grad)
  else
    error("Invalid mid3 position")
  end
end

"""
    line_seg(x0::T,x1::T,y1::T,x2::T,y2::T)

Calculates the value of the line segment between `(x1,y1)` and `(x2,y2)` at `x = x0`.
"""
function line_seg(x0::T,x1::T,y1::T,x2::T,y2::T) where {T<:AbstractFloat}
   if (x2-x1) == zero(T)
     return y1
   else
     return y1*((x2-x0)/(x2-x1)) + y2*((x0-x1)/(x2-x1))
   end
end

"""
    dline_seg(x0::T,x1::T,y1::T,x2::T,y2::T)

Calculates the value of the slope line segment between `(x1,y1)` and `(x2,y2)`
defaults to evaluating the derivative of the function if the interval is tight.
"""
function dline_seg(x0::T,x1::T,y1::T,x2::T,y2::T,d::T) where {T<:AbstractFloat}
    if (x2 == x1)
      return d
    else
      return (y2-y1)/(x2-x1)
    end
end

"""
    grad_calc(cv::T,cc::T,int1::Int64,int2::Int64,dcv::SVector{N,T},dcc::SVector{N,T}) where {N,T}

(Sub)gradient calculation function. Takes the convex gradient, 'cv', the
concave gradient, 'cc', the mid index values 'int1,int2', and the derivative of
the convex and concave envelope functions 'dcv,dcc'.
"""
function grad_calc(cv::SVector{N,T},cc::SVector{N,T},int1::Int64,int2::Int64,dcv::T,dcc::T) where {N,T<:AbstractFloat}
  cv_grad::SVector{N,T} = dcv*( int1==1 ? cv :( int1==2 ? cv : zeros(cv)))
  cc_grad::SVector{N,T} = dcc*( int2==1 ? cc :( int2==2 ? cc : zeros(cv)))
  return cv_grad, cc_grad
end

"""
    tighten_subgrad(cc,cv,cc_grad,cv_grad,Xintv,Xbox,xref)

Tightens the interval bounds using subgradients. Inputs:
* `cc::T`: concave bound
* `cv::T`: convex bound
* `cc_grad::SVector{N,T}`: subgradient/gradient of concave bound
* `cv_grad::SVector{N,T}`: subgradient/gradient of convex bound
* `Xintv::Interval{T}`: Interval domain of function
* `Xbox::Vector{Interval{T}}`: Original decision variable bounds
* `xref::Vector{T}`: Reference point in Xbox
"""
function tighten_subgrad(cc::T,cv::T,cc_grad::SVector{N,T},cv_grad::SVector{N,T},
                         Xintv::Interval{T},Xbox::Vector{Interval{T}},xref::Vector{T}) where {N,T<:AbstractFloat}
  if (length(Xbox)>0 && Xbox[1]!=∅)
    upper_refine::Interval{T} = Interval{T}(cc)
    lower_refine::Interval{T} = Interval{T}(cv)
    for i=1:length(Xbox)
      upper_refine = upper_refine + cc_grad[i]*(Xbox[i]-xref[i])
      lower_refine = lower_refine + cv_grad[i]*(Xbox[i]-xref[i])
    end
    return Interval{T}(max(lower_refine.lo,Xintv.lo),min(upper_refine.hi,Xintv.hi))
  else
    return Xintv
  end
end

"""
    outer_rnd!(Intv::Interval{T})

Outer rounds the interval `Intv` by `MC_param.outer_param`.
"""
function outer_rnd(Intv::Interval{T}) where {T<:AbstractFloat}
  return Interval{T}(Intv.lo-MC_param.outer_param, Intv.hi+MC_param.outer_param)
end

"""
    isequal(x::S,y::S,atol::T,rtol::T)

Checks that `x` and `y` are equal to with absolute tolerance `atol` and relative
tolerance `rtol`.
"""
function isequal(x::S,y::S,atol::T,rtol::T) where {S,T<:AbstractFloat}
  return (abs(x-y) < (atol + 0.5*abs(x+y)*rtol))
end

"""
    Intv(x::SMCg{N,T})
"""
function Intv(x::SMCg{N,T}) where {N,T<:AbstractFloat}
  return x.Intv
end


text:
# get unique values of css_sys in graph.jl.
# It is denoted here as cs.
function cstopcs(   cs::Vector{Vector{T}},
    cs_LUT::Vector{Vector{Vector{Int}}}) where T <: Real

    N_subsystems = length(cs)
    p_cs = Vector{Vector{T}}(undef, N_subsystems)

    for m = 1:N_subsystems

    p_cs[m] = Vector{T}(undef, length(cs_LUT[m]))

    for i = 1:length(cs_LUT[m])
    inds = cs_LUT[m][i]

    # take average.
    p_cs[m][i] = sum( cs[m][k] for k in inds)/length(inds)
    end

    end

    return p_cs
end

function pcstocs(   p_cs::Vector{Vector{T}},
    cs_LUT::Vector{Vector{Vector{Int}}},
    cs_len_sys::Vector{Int}) where T <: Real

    N_subsystems = length(p_cs)

    cs = Vector{Vector{T}}(undef, N_subsystems)
    pcstocs!(cs, p_cs, cs_LUT, cs_len_sys)

    return cs
end

function pcstocs!(  cs::Vector{Vector{T}},
    p_cs::Vector{Vector{T}},
    cs_LUT::Vector{Vector{Vector{Int}}},
    cs_len_sys::Vector{Int}) where T <: Real
    #
    N_subsystems = length(p_cs)

    for m = 1:N_subsystems

    cs[m] = Vector{T}(undef, cs_len_sys[m])

    # for i = 1:length(p_cs[m])
    #     inds = cs_LUT[m][i]
    #
    #     for k in inds
    #         cs[m][k] = p_cs[m][i]
    #     end
    # end
    pcstocs!(cs[m], p_cs[m], cs_LUT[m])
end

return nothing
end

function pcstocs!(  cs_m::Vector{T},
    p_cs_m::Vector{T},
    cs_LUT_m::Vector{Vector{Int}}) where T <: Real

    @assert length(p_cs_m) == length(cs_LUT_m) <= length(cs_m)

    for i = 1:length(p_cs_m)
    inds = cs_LUT_m[i]

    for k in inds
    cs_m[k] = p_cs_m[i]
    end
    end

    return nothing
end

ext: jl


text:
push!(LOAD_PATH,"../src/")

using Documenter
using DocumenterCitations
using Plots
using HOODESolver

ENV["GKSwstype"] = "100"

bib = CitationBibliography(joinpath(@__DIR__, "references.bib"))
makedocs(
    bib,
    sitename = "HOODESolver.jl",
    authors="Yves Mocquard",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://ymocquar.github.io/HOODESolver.jl",
        assets=String[],
    ),
    modules = [HOODESolver],
    pages = ["Documentation"    => "index.md",
             "Numerical Method" => "numerical_method.md",
             "Quickstart"       => "quickstart.md",
             "Charged Particle" => "charged_particle.md",
             "Future work"      => "future_work.md",
             "Types"            => "types.md",
             "Functions"        => "functions.md"],
    repo = "https://github.com/ymocquar/HOODESolver.jl/blob/{commit}{path}#{line}"
)

deploydocs(;
branch = "gh-pages",
    devbranch = "master",
    repo="github.com/ymocquar/HOODESolver.jl",
)


text:
export CollisionChecker, DiscreteCollisionChecker, SweptCollisionChecker
export is_free_state, is_free_motion, is_free_path, inflate, addobstacle, addblocker, close, closest

abstract type CollisionChecker end
abstract type DiscreteCollisionChecker <: CollisionChecker end
abstract type SweptCollisionChecker <: CollisionChecker end

include("collisioncheckers/SAT2D.jl")
include("collisioncheckers/robots2D.jl")
include("collisioncheckers/bvls.jl")
include("collisioncheckers/boxesND.jl")


text:
module OmegaGrav
using HCubature
using DoubleExponentialFormulas
using Dierckx
using SpecialFunctions
using MatterPower
using HaloMF
export ograv_pk, ograv_halo
export otherm_upp, otherm_ks
export by
include("ograv.jl")
include("otherm.jl")
include("dndlnMh.jl")
end # module


text:
const NUM_CORES = Sys.CPU_THREADS >>> 1

const CACHE_COUNT = (
    NUM_CORES,
    NUM_CORES,
    1,
    0
)

const CACHE_LEVELS = 3

"""
L₁, L₂, L₃, L₄ cache size

Warning: this is a generic fallback
Assuming
Split: 32 KiB L₁/core, 0.5 MiB L₂ /core
Shared 1 MiB L₃/core
"""
const CACHE_SIZE = ( 32 * (1 << 10), (1 << 20) >>> 1, (1 << 20) * NUM_CORES )

const L₁CACHE = (size = CACHE_SIZE[1], depth = 1, linesize = 64, associativity = 8, type = :Data)
const L₂CACHE = (size = CACHE_SIZE[2], depth = 2, linesize = 64, associativity = 8, type = :Unified)
const L₃CACHE = (size = CACHE_SIZE[3], depth = 3, linesize = 64, associativity = 8, type = :Unified)
const L₄CACHE = (size = nothing, depth = nothing, linesize = nothing, associativity = nothing, type = nothing)


text:
abstract type AbstractQRLayout <: MemoryLayout end

"""
   QRCompactWYLayout{SLAY,TLAY}()

represents a Compact-WY QR factorization whose
factors are stored with layout SLAY and τ stored with layout TLAY
"""
struct QRCompactWYLayout{SLAY,TLAY} <: AbstractQRLayout end
"""
    QRPackedLayout{SLAY,TLAY}()

represents a Packed QR factorization whose
factors are stored with layout SLAY and τ stored with layout TLAY
"""
struct QRPackedLayout{SLAY,TLAY} <: AbstractQRLayout end

MemoryLayout(::Type{<:LinearAlgebra.QRCompactWY{<:Any,MAT}}) where MAT =
    QRCompactWYLayout{typeof(MemoryLayout(MAT)),DenseColumnMajor}()
MemoryLayout(::Type{<:LinearAlgebra.QR{<:Any,MAT}}) where MAT =
    QRPackedLayout{typeof(MemoryLayout(MAT)),DenseColumnMajor}()

function materialize!(L::Ldiv{<:QRCompactWYLayout,<:Any,<:Any,<:AbstractVector})
    A,b = L.A, L.B
    ldiv!(UpperTriangular(A.R), view(lmul!(adjoint(A.Q), b), 1:size(A, 2)))
    b
end

function materialize!(L::Ldiv{<:QRCompactWYLayout,<:Any,<:Any,<:AbstractMatrix})
    A,B = L.A, L.B
    ldiv!(UpperTriangular(A.R), view(lmul!(adjoint(A.Q), B), 1:size(A, 2), 1:size(B, 2)))
    B
end

# Julia implementation similar to xgelsy
function materialize!(Ldv::Ldiv{<:QRPackedLayout,<:Any,<:Any,<:AbstractMatrix{T}}) where T
    A,B = Ldv.A,Ldv.B
    m, n = size(A)
    minmn = min(m,n)
    mB, nB = size(B)
    lmul!(adjoint(A.Q), view(B, 1:m, :))
    mB ≥ n || throw(DimensionMismatch("Number of rows in B must exceed number of columns in A"))
    R = A.R
    @inbounds begin
        if n > m # minimum norm solution
            τ = zeros(T,m)
            for k = m:-1:1 # Trapezoid to triangular by elementary operation
                x = view(R, k, [k; m + 1:n])
                τk = reflector!(x)
                τ[k] = conj(τk)
                for i = 1:k - 1
                    vRi = R[i,k]
                    for j = m + 1:n
                        vRi += R[i,j]*x[j - m + 1]'
                    end
                    vRi *= τk
                    R[i,k] -= vRi
                    for j = m + 1:n
                        R[i,j] -= vRi*x[j - m + 1]
                    end
                end
            end
        end
        ldiv!(UpperTriangular(view(R, :, 1:minmn)), view(B, 1:minmn, :))
        if n > m # Apply elementary transformation to solution
            B[m + 1:mB,1:nB] .= zero(T)
            for j = 1:nB
                for k = 1:m
                    vBj = B[k,j]
                    for i = m + 1:n
                        vBj += B[i,j]*R[k,i]'
                    end
                    vBj *= τ[k]
                    B[k,j] -= vBj
                    for i = m + 1:n
                        B[i,j] -= R[k,i]*vBj
                    end
                end
            end
        end
    end
    return B
end
materialize!(Ldv::Ldiv{<:QRPackedLayout,<:Any,<:Any,<:AbstractVector{T}}) where T =
    ldiv!(Ldv.A, reshape(Ldv.B, length(Ldv.B), 1))[:]

abstract type AbstractQLayout <: MemoryLayout end
struct QRPackedQLayout{SLAY,TLAY} <: AbstractQLayout end
struct AdjQRPackedQLayout{SLAY,TLAY} <: AbstractQLayout end
struct QRCompactWYQLayout{SLAY,TLAY} <: AbstractQLayout end
struct AdjQRCompactWYQLayout{SLAY,TLAY} <: AbstractQLayout end

MemoryLayout(::Type{<:LinearAlgebra.QRPackedQ{<:Any,S}}) where {S,T} =
    QRPackedQLayout{typeof(MemoryLayout(S)),DenseColumnMajor}()
MemoryLayout(::Type{<:LinearAlgebra.QRCompactWYQ{<:Any,S}}) where {S,T} =
    QRCompactWYQLayout{typeof(MemoryLayout(S)),DenseColumnMajor}()

adjointlayout(::Type, ::QRPackedQLayout{SLAY,TLAY}) where {SLAY,TLAY} = AdjQRPackedQLayout{SLAY,TLAY}()
adjointlayout(::Type, ::QRCompactWYQLayout{SLAY,TLAY}) where {SLAY,TLAY} = AdjQRCompactWYQLayout{SLAY,TLAY}()

copy(M::Lmul{<:AbstractQLayout}) = copyto!(similar(M), M)

function copyto!(dest::AbstractArray{T}, M::Lmul{<:AbstractQLayout}) where T
    A,B = M.A,M.B
    if size(dest,1) == size(B,1)
        copyto!(dest, B)
    else
        copyto!(view(dest,1:size(B,1),:), B)
        zero!(@view(dest[size(B,1)+1:end,:]))
    end
    lmul!(A,dest)
end

function copyto!(dest::AbstractArray, M::Ldiv{<:AbstractQLayout})
    A,B = M.A,M.B
    copyto!(dest, B)
    ldiv!(A,dest)
end

materialize!(M::Lmul{LAY}) where LAY<:AbstractQLayout = error("Overload materialize!(::Lmul{$(LAY)})")
materialize!(M::Rmul{LAY}) where LAY<:AbstractQLayout = error("Overload materialize!(::Rmul{$(LAY)})")

materialize!(M::Ldiv{<:AbstractQLayout}) = materialize!(Lmul(M.A',M.B))

materialize!(M::Lmul{<:QRPackedQLayout{<:AbstractColumnMajor,<:AbstractColumnMajor},<:AbstractColumnMajor,<:AbstractMatrix{T},<:AbstractVecOrMat{T}}) where T<:BlasFloat =
    LAPACK.ormqr!('L','N',M.A.factors,M.A.τ,M.B)

materialize!(M::Lmul{<:QRCompactWYQLayout{<:AbstractColumnMajor,<:AbstractColumnMajor},<:AbstractColumnMajor,<:AbstractMatrix{T},<:AbstractVecOrMat{T}}) where T<:BlasFloat =
    LAPACK.gemqrt!('L','N',M.A.factors,M.A.T,M.B)

function materialize!(M::Lmul{<:QRCompactWYQLayout{<:AbstractColumnMajor,<:AbstractColumnMajor},<:AbstractRowMajor,<:AbstractMatrix{T},<:AbstractMatrix{T}}) where T<:BlasReal
    LAPACK.gemqrt!('R','T',M.A.factors,M.A.T,transpose(M.B))
    M.B
end
function materialize!(M::Lmul{<:QRCompactWYQLayout{<:AbstractColumnMajor,<:AbstractColumnMajor},<:ConjLayout{<:AbstractRowMajor},<:AbstractMatrix{T},<:AbstractMatrix{T}}) where T<:BlasComplex
    LAPACK.gemqrt!('R','C',M.A.factors,M.A.T,(M.B)')
    M.B
end

function materialize!(M::Lmul{<:QRPackedQLayout})
    A,B = M.A, M.B
    require_one_based_indexing(B)
    mA, nA = size(A.factors)
    mB, nB = size(B,1), size(B,2)
    if mA != mB
        throw(DimensionMismatch("matrix A has dimensions ($mA,$nA) but B has dimensions ($mB, $nB)"))
    end
    Afactors = A.factors
    @inbounds begin
        for k = min(mA,nA):-1:1
            for j = 1:nB
                vBj = B[k,j]
                for i = k+1:mB
                    vBj += conj(Afactors[i,k])*B[i,j]
                end
                vBj = A.τ[k]*vBj
                B[k,j] -= vBj
                for i = k+1:mB
                    B[i,j] -= Afactors[i,k]*vBj
                end
            end
        end
    end
    B
end

### QcB
materialize!(M::Lmul{<:AdjQRPackedQLayout{<:AbstractStridedLayout,<:AbstractStridedLayout},<:AbstractStridedLayout,<:AbstractMatrix{T},<:AbstractVecOrMat{T}}) where T<:BlasFloat =
    (A = M.A.parent; LAPACK.ormqr!('L','T',A.factors,A.τ,M.B))
materialize!(M::Lmul{<:AdjQRPackedQLayout{<:AbstractStridedLayout,<:AbstractStridedLayout},<:AbstractStridedLayout,<:AbstractMatrix{T},<:AbstractVecOrMat{T}}) where T<:BlasComplex =
    (A = M.A.parent; LAPACK.ormqr!('L','C',A.factors,A.τ,M.B))
materialize!(M::Lmul{<:AdjQRCompactWYQLayout{<:AbstractStridedLayout,<:AbstractStridedLayout},<:AbstractStridedLayout,<:AbstractMatrix{T},<:AbstractVecOrMat{T}}) where T<:BlasFloat =
    (A = M.A.parent; LAPACK.gemqrt!('L','T',A.factors,A.T,M.B))
materialize!(M::Lmul{<:AdjQRCompactWYQLayout{<:AbstractStridedLayout,<:AbstractStridedLayout},<:AbstractStridedLayout,<:AbstractMatrix{T},<:AbstractVecOrMat{T}}) where T<:BlasComplex =
    (A = M.A.parent; LAPACK.gemqrt!('L','C',A.factors,A.T,M.B))
function materialize!(M::Lmul{<:AdjQRPackedQLayout})
    adjA,B = M.A, M.B
    require_one_based_indexing(B)
    A = adjA.parent
    mA, nA = size(A.factors)
    mB, nB = size(B,1), size(B,2)
    if mA != mB
        throw(DimensionMismatch("matrix A has dimensions ($mA,$nA) but B has dimensions ($mB, $nB)"))
    end
    Afactors = A.factors
    @inbounds begin
        for k = 1:min(mA,nA)
            for j = 1:nB
                vBj = B[k,j]
                for i = k+1:mB
                    vBj += conj(Afactors[i,k])*B[i,j]
                end
                vBj = conj(A.τ[k])*vBj
                B[k,j] -= vBj
                for i = k+1:mB
                    B[i,j] -= Afactors[i,k]*vBj
                end
            end
        end
    end
    B
end

## AQ
materialize!(M::Rmul{<:AbstractStridedLayout,<:QRPackedQLayout{<:AbstractStridedLayout,<:AbstractStridedLayout},<:AbstractVecOrMat{T},<:AbstractMatrix{T}}) where T<:BlasFloat =
    LAPACK.ormqr!('R', 'N', M.B.factors, M.B.τ, M.A)
materialize!(M::Rmul{<:AbstractStridedLayout,<:QRCompactWYQLayout{<:AbstractStridedLayout,<:AbstractStridedLayout},<:AbstractVecOrMat{T},<:AbstractMatrix{T}}) where T<:BlasFloat =
    LAPACK.gemqrt!('R','N', M.B.factors, M.B.T, M.A)
function materialize!(M::Rmul{<:Any,<:QRPackedQLayout})
    A,Q = M.A,M.B
    mQ, nQ = size(Q.factors)
    mA, nA = size(A,1), size(A,2)
    if nA != mQ
        throw(DimensionMismatch("matrix A has dimensions ($mA,$nA) but matrix Q has dimensions ($mQ, $nQ)"))
    end
    Qfactors = Q.factors
    @inbounds begin
        for k = 1:min(mQ,nQ)
            for i = 1:mA
                vAi = A[i,k]
                for j = k+1:mQ
                    vAi += A[i,j]*Qfactors[j,k]
                end
                vAi = vAi*Q.τ[k]
                A[i,k] -= vAi
                for j = k+1:nA
                    A[i,j] -= vAi*conj(Qfactors[j,k])
                end
            end
        end
    end
    A
end

### AQc
materialize!(M::Rmul{<:AbstractStridedLayout,<:AdjQRPackedQLayout{<:AbstractStridedLayout,<:AbstractStridedLayout},<:AbstractVecOrMat{T},<:AbstractMatrix{T}}) where T<:BlasReal =
    (B = M.B.parent; LAPACK.ormqr!('R','T',B.factors,B.τ,M.A))
materialize!(M::Rmul{<:AbstractStridedLayout,<:AdjQRPackedQLayout{<:AbstractStridedLayout,<:AbstractStridedLayout},<:AbstractVecOrMat{T},<:AbstractMatrix{T}}) where T<:BlasComplex =
    (B = M.B.parent; LAPACK.ormqr!('R','C',B.factors,B.τ,M.A))
materialize!(M::Rmul{<:AbstractStridedLayout,<:AdjQRCompactWYQLayout{<:AbstractStridedLayout,<:AbstractStridedLayout},<:AbstractVecOrMat{T},<:AbstractMatrix{T}}) where T<:BlasReal =
    (B = M.B.parent; LAPACK.gemqrt!('R','T',B.factors,B.T,M.A))
materialize!(M::Rmul{<:AbstractStridedLayout,<:AdjQRCompactWYQLayout{<:AbstractStridedLayout,<:AbstractStridedLayout},<:AbstractVecOrMat{T},<:AbstractMatrix{T}}) where T<:BlasComplex =
    (B = M.B.parent; LAPACK.gemqrt!('R','C',B.factors,B.T,M.A))
function materialize!(M::Rmul{<:Any,<:AdjQRPackedQLayout})
    A,adjQ = M.A,M.B
    Q = adjQ.parent
    mQ, nQ = size(Q.factors)
    mA, nA = size(A,1), size(A,2)
    if nA != mQ
        throw(DimensionMismatch("matrix A has dimensions ($mA,$nA) but matrix Q has dimensions ($mQ, $nQ)"))
    end
    Qfactors = Q.factors
    @inbounds begin
        for k = min(mQ,nQ):-1:1
            for i = 1:mA
                vAi = A[i,k]
                for j = k+1:mQ
                    vAi += A[i,j]*Qfactors[j,k]
                end
                vAi = vAi*conj(Q.τ[k])
                A[i,k] -= vAi
                for j = k+1:nA
                    A[i,j] -= vAi*conj(Qfactors[j,k])
                end
            end
        end
    end
    A
end

_qr(layout, axes, A; kwds...) = Base.invoke(qr, Tuple{AbstractMatrix{eltype(A)}}, A; kwds...)
_qr(layout, axes, A, pivot::P; kwds...) where P = Base.invoke(qr, Tuple{AbstractMatrix{eltype(A)},P}, A, pivot; kwds...)
_lu(layout, axes, A; kwds...) = Base.invoke(lu, Tuple{AbstractMatrix{eltype(A)}}, A; kwds...)
_lu(layout, axes, A, pivot::P; kwds...) where P = Base.invoke(lu, Tuple{AbstractMatrix{eltype(A)},P}, A, pivot; kwds...)
_qr!(layout, axes, A, args...; kwds...) = error("Overload _qr!(::$(typeof(layout)), axes, A)")
_lu!(layout, axes, A, args...; kwds...) = error("Overload _lu!(::$(typeof(layout)), axes, A)")
_factorize(layout, axes, A) = qr(A) # Default to QR

_inv_eye(_, ::Type{T}, axs::NTuple{2,Base.OneTo{Int}}) where T = Matrix{T}(I, map(length,axs)...)
function _inv_eye(A, ::Type{T}, (rows,cols)) where T
    dest = zero!(similar(A, T, (cols,rows)))
    dest[diagind(dest)] .= one(T)
    dest
end

function _inv(layout, axes, A)
    T = eltype(A)
    (rows,cols) = axes
    n = checksquare(A)
    S = typeof(zero(T)/one(T))      # dimensionful
    S0 = typeof(zero(T)/oneunit(T)) # dimensionless
    dest = _inv_eye(A, S0, (cols,rows))
    ldiv!(factorize(convert(AbstractMatrix{S}, A)), dest)
end

macro _layoutfactorizations(Typ)
    esc(quote
        LinearAlgebra.qr(A::$Typ, args...; kwds...) = ArrayLayouts._qr(ArrayLayouts.MemoryLayout(A), axes(A), A, args...; kwds...)
        LinearAlgebra.qr!(A::$Typ, args...; kwds...) = ArrayLayouts._qr!(ArrayLayouts.MemoryLayout(A), axes(A), A, args...; kwds...)
        LinearAlgebra.lu(A::$Typ, pivot::Union{Val{false}, Val{true}}; kwds...) = ArrayLayouts._lu(ArrayLayouts.MemoryLayout(A), axes(A), A, pivot; kwds...)
        LinearAlgebra.lu(A::$Typ{T}; kwds...) where T = ArrayLayouts._lu(ArrayLayouts.MemoryLayout(A), axes(A), A; kwds...)
        LinearAlgebra.lu!(A::$Typ, args...; kwds...) = ArrayLayouts._lu!(ArrayLayouts.MemoryLayout(A), axes(A), A, args...; kwds...)
        LinearAlgebra.factorize(A::$Typ) = ArrayLayouts._factorize(ArrayLayouts.MemoryLayout(A), axes(A), A)
        LinearAlgebra.inv(A::$Typ) = ArrayLayouts._inv(ArrayLayouts.MemoryLayout(A), axes(A), A)
    end)
end

macro layoutfactorizations(Typ)
    esc(quote
        ArrayLayouts.@_layoutfactorizations $Typ
        ArrayLayouts.@_layoutfactorizations SubArray{<:Any,2,<:$Typ}
    end)
end


text:
# Julia program to Right Rotate a number by a specific bits.

# Right Rotate 'cnt' number of bits of the given number 'n'
function right_rotate_bits(n, cnt)
    cnt = cnt % 31
    while(cnt > 0)
        # Store the current LSB in a temporary variable
        lsb = n & 1
        # Right rotate the given number by one bit and drop its LSB
        n = (n >> 1) & (~(1 << 31))
        # Set the dropped LSB as the new MSB
        n = n | (lsb << 31)
        # Decrement cnt
        cnt = cnt - 1
    end
    return n
end

print("Enter the number? ")
num = readline()
num = parse(Int, num)
print("How many bits do you want to rotate? ")
cnt = readline()
cnt = parse(Int, cnt)
left = right_rotate_bits(num, cnt)
print("The Right-rotated number is: $left")

"""
Time Complexity: O(n)
Space Complexity: O(1)

SAMPLE INPUT AND OUTPUT

Enter the number? 39
How many bits do you want to rotate? 17
The Right-rotated number is: 1277952
"""

max_issues_repo_name: adarshnjena/NeoAlgo

max_forks_repo_name: adarshnjena/NeoAlgo


text:
import Optim

export ## Types
       Duffing,
       Resonator,
       TunableDuffingTransmon,
       TunableFullTransmon,
       TunableTransmon,
       Qubit,
       ## Methods
       hamiltonian

#Resonator
mutable struct Resonator <: QSystem
    label::AbstractString
    freq::Float64
    dim::Int
end
hamiltonian(r::Resonator) = r.freq*number(r)

#Fixed frequency transmon
mutable struct Transmon <: QSystem
    label::AbstractString
    E_C::Float64
    E_J::Float64
    dim::Int
end

"""
Transmon(label::AbstractString, ν::Float64, α::Float64; dim=3)

#Arguments
* ν - 0 ↔ 1 frequency
* α - anharmonicity
"""
# function Transmon(label::AbstractString, ν::Float64, α::Float64; dim=3)
#     E_C = abs(α)
#     E_J = (ν + E_C)^2 / (8 * E_C)
#     return Transmon(label, E_C, E_J, dim)
# end

# hamiltonian(t::Transmon) = sqrt(8*t.E_J*t.E_C)*number(t) - t.E_C/12*(X(t)^4)
function hamiltonian(t::Transmon)
  N = floor(Int, dim(t)/2)
  4 * t.E_C * diagm(0 => (-N:N).^2) - t.E_J  * 0.5 * (diagm(-1 => ones(dim(t)-1)) + diagm(1 => ones(dim(t)-1)))
end

"""
Fit  E_C and E_J for a fixed frequency transmon given f_01 and anharmonicity
"""
function fit_fixed_transmon(f_01, α, dim)
    # helper function for least squares difference between measured and calculated
    function f(params)
      E_C = params[1]
      E_J  = params[2]
      t = QSimulator.Transmon("dummy", E_C, E_J, dim)
      levels = eigvals(hamiltonian(t))
      test_f_01 = levels[2] - levels[1]
      test_f_12 = levels[3] - levels[2]
      ( test_f_01 - f_01)^2 + ( test_f_12 -  test_f_01 - α)^2
    end
    res = Optim.optimize(f, [abs(α), (f_01+abs(α))^2 / (8 * abs(α))] )
    return res.minimizer[1], res.minimizer[2]
end

#Tunable transmon
abstract type TunableTransmon <: QSystem end

mutable struct TunableFullTransmon <: TunableTransmon
    label::AbstractString
    E_C::Float64
    E_J::Float64 #sum of junction E_J's
    d::Float64 #asymmetry parameter
    dim::Int
    flux::Float64 #total flux in units of Phi_0
end

#Helper function to calculate effective EJ for a transmon
# scale_EJ(flux::Float64, d::Float64) = cos(pi*flux)*sqrt(1 + d^2*(tan(pi*flux)^2))

scale_EJ(flux::Float64, d::Float64) = sqrt(cos(π*flux)^2 + d^2*(sin(π*flux)^2))

function hamiltonian(tt::TunableFullTransmon, t::Float64=0.0)
    myE_J = tt.E_J*scale_EJ(tt.flux, tt.d)
    return (sqrt(8*tt.E_C*myE_J)*number(tt) - (1.0/12)*tt.E_C*(X(tt)^4))
end

#Tunable Duffing approx. transmon
mutable struct TunableDuffingTransmon <: TunableTransmon
    label::AbstractString
    E_C::Float64
    E_J::Float64 #sum of junction E_J's
    d::Float64 #asymmetry parameter
    dim::Int
    flux::Float64 #total flux in units of Phi_0
end

function hamiltonian(tt::TunableDuffingTransmon, t::Float64=0.0)
    myE_J = tt.E_J*scale_EJ(tt.flux, tt.d)
    omega_p = sqrt(8*tt.E_C*myE_J)
    omega = [(omega_p-tt.E_C/2)*i-tt.E_C/2*i^2 for i in 0:(tt.dim-1)]
    return diagm(0 => omega)
end

#Basic two-level qubit
mutable struct Qubit <: QSystem
    label::AbstractString
    freq::Float64
end
dim(q::Qubit) = 2
hamiltonian(q::Qubit) = q.freq*number(q)

#Duffing oscillator
mutable struct Duffing <: QSystem
    label::AbstractString
    freq::Float64
    alpha::Float64
    dim::Int
end
hamiltonian(s::Duffing) = (s.freq - 0.5*s.alpha)*number(s) + 0.5*s.alpha * number(s)^2

"""
Fit E_C, E_J and d for a tunable transmon given f_max, f_min and anharmonicty
"""
function fit_tunable_transmon()
  0
end


text:
"""
    DifferentialRKHSMethod1d(epsilon, lambda)

Differential loss reproducing kernel Hilbert space (RKHS) method from [1], Section III.

### Parameters

- `epsilon`: variance of the Gaussian kernel, Eq. 17 in [1]
- `lambda`: regularization parameter, Eq. 20 in [1]

### Reference

[1] Radhakrishnan, A. & and Meyn, S. (2018). Feedback particle filter design using a differential-loss reproducing kernel Hilbert space. 2018 Annual American Control Conference (ACC). IEEE. https://doi.org/10.23919/ACC.2018.8431689
"""
struct DifferentialRKHSMethod <: GainEstimationMethod
    epsilon::Float64
    lambda::Float64
end

# Helper functions for DifferentialRKHSMethod1d
function DifferentialRKHS_MakeMatrices(vec::AbstractVector, eps::Number)
    # constructs kernel matrices for DifferentialRKHSMethod1d
    L = length(vec)
    eps1 = 2*eps
    eps2 = eps^2

    mat = zeros(Float64, L, L, 3)

    @inbounds @simd for i in 1:L
                mat[i,i,1] = 1
                mat[i,i,3] = 1/eps
                @simd for j in i+1:L
                    diff = vec[j] - vec[i]
                    sq = diff*diff
                    K = exp(-sq/eps1)
                    Kx = diff*K/eps
                    Kxy = (eps - sq) * K / eps2
                    mat[i,j,1] = mat[j,i,1] = K
                    mat[i,j,2] = Kx
                    mat[j,i,2] = -Kx
                    mat[i,j,3] = mat[j,i,3] = Kxy
                end
    end

    return view(mat, :, :, 1), view(mat, :, :, 2), view(mat, :, :, 3)
end

function DifferentialRKHS_M_and_b!(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix, v::AbstractVector, lamb::Number)
    # constructs matrix M and vector b for DifferentialRKHSMethod1d
    L = size(A,1)

    b = zeros(Float64, 2*L)
    LinearAlgebra.mul!(view(b,   1:L  ), A, v)
    LinearAlgebra.mul!(view(b, L+1:2*L), B, v)

    B2 = B*B
    BC = B*C
    C2 = C*C
    LinearAlgebra.lmul!(lamb, A)
    LinearAlgebra.lmul!(lamb, B)
    LinearAlgebra.lmul!(lamb, C)

    M = zeros(Float64, 2*L, 2*L)
    @simd for j in 1:L
        @simd for i in 1:L
            @inbounds M[i  ,j  ] =  A[i,j] - B2[i,j]
            @inbounds M[i+L,j  ] =  B[i,j] + BC[j,i]
            @inbounds M[i  ,j+L] = -B[i,j] - BC[i,j]
            @inbounds M[i+L,j+L] =  C[i,j] + C2[i,j]
        end
    end

    return M, b
end

function solve!(eq::PoissonEquation, method::DifferentialRKHSMethod)
    eps = method.epsilon
    N = no_of_particles(eq)
    lamb = N * method.lambda

    H = Htilde(eq)[1,:]

    # Gaussian kernel and partial derivative matrices
    K, Kx, Kxy = DifferentialRKHS_MakeMatrices(eq.positions[1,:], eps)

    # compute M and b
    M, b = DifferentialRKHS_M_and_b!(K, Kx, Kxy, H, lamb) # warning: this function multiplies K, Kx, and Kxy by lamb

    # solve linear system
    beta = M \ b

    # write potential and gain
    eq.potential[1,:] = ( K*view(beta, 1:N) - Kx*view(beta, N+1:2*N) ) / lamb
    broadcast!(-, eq.potential, eq.potential, Statistics.mean(eq.potential, dims=2))
    eq.gain[1,:,1] = ( Kx * view(beta, 1:N) + Kxy * view(beta, N+1:2*N) ) / lamb
end


text:
import MathOptInterface
import SparseArrays

const MOI = MathOptInterface
const CleverDicts = MOI.Utilities.CleverDicts

"""
    HighsBasisStatus

https://github.com/ERGO-Code/HiGHS/blob/4a5dd7499522f1fa730a31c59bba419b2bcc6839/src/lp_data/HConst.h#L151-L159
"""
@enum(HighsBasisStatus, kLower = 0, kBasic, kUpper, kZero, kNonbasic)

"""
    HighsMatrixFormat

https://github.com/ERGO-Code/HiGHS/blob/4a5dd7499522f1fa730a31c59bba419b2bcc6839/src/lp_data/HConst.h#L88
"""
@enum(HighsMatrixFormat, kNone = 0, kColwise, kRowwise)

"""
    HighsModelStatus

https://github.com/ERGO-Code/HiGHS/blob/4a5dd7499522f1fa730a31c59bba419b2bcc6839/src/lp_data/HConst.h#L126-L148
"""
@enum(
    HighsModelStatus,
    kNotset = 0,
    kLoadError,
    kModelError,
    kPresolveError,
    kSolveError,
    kPostsolveError,
    kModelEmpty,
    kOptimal,
    kInfeasible,
    kUnboundedOrInfeasible,
    kUnbounded,
    kObjectiveBound,
    kObjectiveTarget,
    kTimeLimit,
    kIterationLimit,
    kUnknown,
)

"""
    HighsObjSense

https://github.com/ERGO-Code/HiGHS/blob/4a5dd7499522f1fa730a31c59bba419b2bcc6839/src/lp_data/HConst.h#L86
"""
@enum(HighsObjSense, kMinimize = 1, kMaximize = -1)

"""
    HighsVartype

    https://github.com/ERGO-Code/HiGHS/blob/4a5dd7499522f1fa730a31c59bba419b2bcc6839/src/lp_data/HConst.h#L69-L73
"""
@enum(HighsVartype, kContinuous = 0, kInteger = 1, kImplicitInteger = 2)

@enum(
    _RowType,
    _ROW_TYPE_LESSTHAN,
    _ROW_TYPE_GREATERTHAN,
    _ROW_TYPE_INTERVAL,
    _ROW_TYPE_EQUAL_TO,
)

_row_type(::MOI.GreaterThan{Float64}) = _ROW_TYPE_GREATERTHAN
_row_type(::MOI.LessThan{Float64}) = _ROW_TYPE_LESSTHAN
_row_type(::MOI.EqualTo{Float64}) = _ROW_TYPE_EQUAL_TO
_row_type(::MOI.Interval{Float64}) = _ROW_TYPE_INTERVAL

_bounds(s::MOI.GreaterThan{Float64}) = s.lower, Inf
_bounds(s::MOI.LessThan{Float64}) = -Inf, s.upper
_bounds(s::MOI.EqualTo{Float64}) = s.value, s.value
_bounds(s::MOI.Interval{Float64}) = s.lower, s.upper

@enum(
    _BoundEnum,
    _BOUND_NONE,
    _BOUND_LESS_THAN,
    _BOUND_GREATER_THAN,
    _BOUND_LESS_AND_GREATER_THAN,
    _BOUND_INTERVAL,
    _BOUND_EQUAL_TO,
)

@enum(_TypeEnum, _TYPE_CONTINUOUS, _TYPE_INTEGER, _TYPE_BINARY,)

const _SCALAR_SETS = Union{
    MOI.GreaterThan{Float64},
    MOI.LessThan{Float64},
    MOI.EqualTo{Float64},
    MOI.Interval{Float64},
}

"""
_VariableInfo

A struct to store information about the variables.
"""
mutable struct _VariableInfo
    # We need to keep the index here because sometimes we call `LinearIndex`.
    index::MOI.VariableIndex
    # The variable name.
    name::String
    # The zero-indexed column in the HiGHS object.
    column::Cint
    # Storage to keep track of the variable bounds.
    bound::_BoundEnum
    lower::Float64
    upper::Float64
    # Track integrality
    type::_TypeEnum
    function _VariableInfo(
        index::MOI.VariableIndex,
        column::Cint,
        bound::_BoundEnum = _BOUND_NONE,
    )
        return new(index, "", column, bound, -Inf, Inf, _TYPE_CONTINUOUS)
    end
end

function _update_info(info::_VariableInfo, s::MOI.GreaterThan{Float64})
    _throw_if_existing_lower(info, s)
    if info.bound == _BOUND_LESS_THAN
        info.bound = _BOUND_LESS_AND_GREATER_THAN
    else
        info.bound = _BOUND_GREATER_THAN
    end
    info.lower = s.lower
    return
end

function _update_info(info::_VariableInfo, s::MOI.LessThan{Float64})
    _throw_if_existing_upper(info, s)
    if info.bound == _BOUND_GREATER_THAN
        info.bound = _BOUND_LESS_AND_GREATER_THAN
    else
        info.bound = _BOUND_LESS_THAN
    end
    info.upper = s.upper
    return
end

function _update_info(info::_VariableInfo, s::MOI.EqualTo{Float64})
    _throw_if_existing_lower(info, s)
    _throw_if_existing_upper(info, s)
    info.bound = _BOUND_EQUAL_TO
    info.lower = s.value
    info.upper = s.value
    return
end

function _update_info(info::_VariableInfo, s::MOI.Interval{Float64})
    _throw_if_existing_lower(info, s)
    _throw_if_existing_upper(info, s)
    info.bound = _BOUND_INTERVAL
    info.lower = s.lower
    info.upper = s.upper
    return
end

function _variable_info_dict()
    return CleverDicts.CleverDict{MOI.VariableIndex,_VariableInfo}(
        x -> x.value,
        i -> MOI.VariableIndex(i),
    )
end

"""
    _ConstraintInfo

A struct to store information about the affine constraints.
"""
mutable struct _ConstraintInfo
    # The constraint name.
    name::String
    # The zero-indexed row in the HiGHS object.
    row::Cint
    # Storage to keep track of the constraint bounds.
    set::_RowType
    lower::Float64
    upper::Float64
end

function _ConstraintInfo(set::_SCALAR_SETS)
    lower, upper = _bounds(set)
    return _ConstraintInfo("", 0, _row_type(set), lower, upper)
end

struct _ConstraintKey
    value::Int64
end

function _constraint_info_dict()
    return CleverDicts.CleverDict{_ConstraintKey,_ConstraintInfo}(
        c -> c.value,
        i -> _ConstraintKey(i),
    )
end

"""
    _set(c::_ConstraintInfo)

Return the set associated with a constraint.
"""
function _set(c::_ConstraintInfo)
    if c.set == _ROW_TYPE_LESSTHAN
        return MOI.LessThan(c.upper)
    elseif c.set == _ROW_TYPE_GREATERTHAN
        return MOI.GreaterThan(c.lower)
    elseif c.set == _ROW_TYPE_INTERVAL
        return MOI.Interval(c.lower, c.upper)
    else
        @assert c.set == _ROW_TYPE_EQUAL_TO
        return MOI.EqualTo(c.lower)
    end
end

@enum(_OptimizeStatus, _OPTIMIZE_NOT_CALLED, _OPTIMIZE_OK, _OPTIMIZE_ERRORED)

"""
    _Solution

A struct to store the vector solution from HiGHS because it doesn't support
accessing them element-wise.
"""
mutable struct _Solution
    status::_OptimizeStatus
    colvalue::Vector{Cdouble}
    coldual::Vector{Cdouble}
    colstatus::Vector{Cint}
    rowvalue::Vector{Cdouble}
    rowdual::Vector{Cdouble}
    rowstatus::Vector{Cint}
    has_primal_solution::Bool
    has_dual_solution::Bool
    has_primal_ray::Bool
    has_dual_ray::Bool
    function _Solution()
        return new(
            _OPTIMIZE_NOT_CALLED,
            Cdouble[],
            Cdouble[],
            Cint[],
            Cdouble[],
            Cdouble[],
            Cint[],
            false,
            false,
            false,
            false,
        )
    end
end

function Base.empty!(x::_Solution)
    x.status = _OPTIMIZE_NOT_CALLED
    empty!(x.colvalue)
    empty!(x.coldual)
    empty!(x.colstatus)
    empty!(x.rowvalue)
    empty!(x.rowdual)
    empty!(x.rowstatus)
    x.has_primal_solution = false
    x.has_dual_solution = false
    x.has_dual_ray = false
    x.has_primal_ray = false
    return x
end

Base.isempty(x::_Solution) = x.status == _OPTIMIZE_NOT_CALLED

mutable struct Optimizer <: MOI.AbstractOptimizer
    # A pointer to the underlying HiGHS optimizer.
    inner::Ptr{Cvoid}

    options::Dict{String,Any}

    # Storage for `MOI.Name`.
    name::String

    # A flag to keep track of MOI.FEASIBILITY_SENSE, since HiGHS only stores
    # MIN_SENSE or MAX_SENSE. This allows us to differentiate between MIN_SENSE
    # and FEASIBILITY_SENSE.
    is_feasibility::Bool
    is_objective_function_set::Bool
    is_objective_sense_set::Bool

    # HiGHS doesn't have special support for binary variables. Cache them here
    # to modify bounds on solve.
    binaries::Set{_VariableInfo}

    variable_info::typeof(_variable_info_dict())
    affine_constraint_info::typeof(_constraint_info_dict())

    # Mappings from variable and constraint names to their indices. These are
    # lazily built on-demand, so most of the time, they are `nothing`.
    name_to_variable::Union{
        Nothing,
        Dict{String,Union{Nothing,MOI.VariableIndex}},
    }
    name_to_constraint_index::Union{
        Nothing,
        Dict{String,Union{Nothing,MOI.ConstraintIndex}},
    }

    # HiGHS just returns a single solution struct :(
    solution::_Solution

    """
        Optimizer()

    Create a new Optimizer object.
    """
    function Optimizer()
        model = new(
            C_NULL,
            Dict{String,Any}(),
            "",
            true,
            false,
            false,
            Set{_VariableInfo}(),
            _variable_info_dict(),
            _constraint_info_dict(),
            nothing,
            nothing,
            _Solution(),
        )
        MOI.empty!(model)
        finalizer(Highs_destroy, model)
        return model
    end
end

Base.cconvert(::Type{Ptr{Cvoid}}, model::Optimizer) = model
Base.unsafe_convert(::Type{Ptr{Cvoid}}, model::Optimizer) = model.inner

function _check_ret(ret::Cint)
    if ret != Cint(0)
        error(
            "Encountered an error in HiGHS (Status $(ret)). Check the log " *
            "for details.",
        )
    end
    return
end

function Base.show(io::IO, model::Optimizer)
    return print(
        io,
        "A HiGHS model with $(Highs_getNumCols(model)) columns and " *
        "$(Highs_getNumRows(model)) rows.",
    )
end

function MOI.empty!(model::Optimizer)
    if model.inner != C_NULL
        Highs_destroy(model)
    end
    model.inner = Highs_create()
    for (key, value) in model.options
        MOI.set(model, MOI.RawOptimizerAttribute(key), value)
    end
    model.name = ""
    model.is_feasibility = true
    model.is_objective_function_set = false
    model.is_objective_sense_set = false
    empty!(model.binaries)
    empty!(model.variable_info)
    empty!(model.affine_constraint_info)
    model.name_to_variable = nothing
    model.name_to_constraint_index = nothing
    empty!(model.solution)
    ret = Highs_changeObjectiveOffset(model, 0.0)
    _check_ret(ret)
    return
end

function MOI.is_empty(model::Optimizer)
    offset = Ref{Cdouble}(0.0)
    ret = Highs_getObjectiveOffset(model, offset)
    _check_ret(ret)
    return Highs_getNumCols(model) == 0 &&
           Highs_getNumRows(model) == 0 &&
           model.is_feasibility &&
           !model.is_objective_function_set &&
           !model.is_objective_sense_set &&
           isempty(model.binaries) &&
           isempty(model.variable_info) &&
           isempty(model.affine_constraint_info) &&
           model.name_to_variable === nothing &&
           model.name_to_constraint_index === nothing &&
           isempty(model.solution) &&
           iszero(offset[])
end

MOI.get(::Optimizer, ::MOI.SolverName) = "HiGHS"

MOI.get(model::Optimizer, ::MOI.RawSolver) = model

function MOI.get(::Optimizer, ::MOI.ListOfVariableAttributesSet)
    return MOI.AbstractVariableAttribute[MOI.VariableName()]
end

function MOI.get(model::Optimizer, ::MOI.ListOfModelAttributesSet)
    attributes = MOI.AbstractModelAttribute[]
    if model.is_objective_sense_set
        push!(attributes, MOI.ObjectiveSense())
    end
    if model.is_objective_function_set
        push!(
            attributes,
            MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        )
    end
    if !isempty(model.name)
        push!(attributes, MOI.Name())
    end
    return attributes
end

function MOI.get(
    ::Optimizer,
    ::MOI.ListOfConstraintAttributesSet{F,S},
) where {F,S}
    attributes = MOI.AbstractConstraintAttribute[]
    if F != MOI.VariableIndex
        return push!(attributes, MOI.ConstraintName())
    end
    return attributes
end

function MOI.get(model::Optimizer, ::MOI.NumberOfConstraints{F,S}) where {F,S}
    # TODO: this could be more efficient.
    return length(MOI.get(model, MOI.ListOfConstraintIndices{F,S}()))
end

function MOI.get(model::Optimizer, ::MOI.ListOfConstraintTypesPresent)
    constraints = Set{Tuple{Type,Type}}()
    for info in values(model.variable_info)
        if info.bound == _BOUND_NONE
        elseif info.bound == _BOUND_LESS_THAN
            push!(constraints, (MOI.VariableIndex, MOI.LessThan{Float64}))
        elseif info.bound == _BOUND_GREATER_THAN
            push!(constraints, (MOI.VariableIndex, MOI.GreaterThan{Float64}))
        elseif info.bound == _BOUND_LESS_AND_GREATER_THAN
            push!(constraints, (MOI.VariableIndex, MOI.LessThan{Float64}))
            push!(constraints, (MOI.VariableIndex, MOI.GreaterThan{Float64}))
        elseif info.bound == _BOUND_EQUAL_TO
            push!(constraints, (MOI.VariableIndex, MOI.EqualTo{Float64}))
        else
            @assert info.bound == _BOUND_INTERVAL
            push!(constraints, (MOI.VariableIndex, MOI.Interval{Float64}))
        end
        if info.type == _TYPE_INTEGER
            push!(constraints, (MOI.VariableIndex, MOI.Integer))
        elseif info.type == _TYPE_BINARY
            push!(constraints, (MOI.VariableIndex, MOI.ZeroOne))
        end
    end
    for info in values(model.affine_constraint_info)
        push!(
            constraints,
            (MOI.ScalarAffineFunction{Float64}, typeof(_set(info))),
        )
    end
    return collect(constraints)
end

###
### MOI.RawOptimizerAttribute
###

function MOI.supports(model::Optimizer, param::MOI.RawOptimizerAttribute)
    typeP = Ref{Cint}()
    return Highs_getOptionType(model, param.name, typeP) == 0
end

function _check_option_status(ret::Cint)
    if ret != 0
        error("Encountered an error in HiGHS: Check the log for details.")
    end
    return
end

function _set_option(model::Optimizer, option::String, value::Bool)
    return Highs_setBoolOptionValue(model, option, Cint(value))
end

function _set_option(model::Optimizer, option::String, value::Integer)
    return Highs_setIntOptionValue(model, option, Cint(value))
end

function _set_option(model::Optimizer, option::String, value::AbstractFloat)
    return Highs_setDoubleOptionValue(model, option, Cdouble(value))
end

function _set_option(model::Optimizer, option::String, value::String)
    return Highs_setStringOptionValue(model, option, value)
end

function MOI.set(model::Optimizer, param::MOI.RawOptimizerAttribute, value)
    if !MOI.supports(model, param)
        throw(MOI.UnsupportedAttribute(param))
    end
    model.options[param.name] = value
    ret = _set_option(model, param.name, value)
    return _check_option_status(ret)
end

### MOI.get

function _get_bool_option(model::Optimizer, option::String)
    value = Ref{Cint}(0)
    ret = Highs_getBoolOptionValue(model, option, value)
    _check_option_status(ret)
    return Bool(value[])
end

function _get_int_option(model::Optimizer, option::String)
    value = Ref{Cint}(0)
    ret = Highs_getIntOptionValue(model, option, value)
    _check_option_status(ret)
    return value[]
end

function _get_double_option(model::Optimizer, option::String)
    value = Ref{Cdouble}()
    ret = Highs_getDoubleOptionValue(model, option, value)
    _check_option_status(ret)
    return value[]
end

function _get_string_option(model::Optimizer, option::String)
    buffer = Vector{Cchar}(undef, 1024)
    bufferP = pointer(buffer)
    GC.@preserve buffer begin
        ret = Highs_getStringOptionValue(model, option, bufferP)
        _check_option_status(ret)
        return unsafe_string(bufferP)
    end
end

function MOI.get(model::Optimizer, param::MOI.RawOptimizerAttribute)
    if !(param.name isa String)
        throw(MOI.UnsupportedAttribute(param))
    end
    typeP = Ref{Cint}()
    ret = Highs_getOptionType(model, param.name, typeP)
    if ret != 0
        throw(MOI.UnsupportedAttribute(param))
    elseif typeP[] == 0
        return _get_bool_option(model, param.name)::Bool
    elseif typeP[] == 1
        return _get_int_option(model, param.name)::Cint
    elseif typeP[] == 2
        return _get_double_option(model, param.name)::Cdouble
    else
        @assert typeP[] == 3
        return _get_string_option(model, param.name)::String
    end
end

###
### MOI.TimeLimitSec
###

MOI.supports(::Optimizer, ::MOI.TimeLimitSec) = true

function MOI.set(model::Optimizer, ::MOI.TimeLimitSec, ::Nothing)
    return MOI.set(model, MOI.RawOptimizerAttribute("time_limit"), Inf)
end

function MOI.set(model::Optimizer, ::MOI.TimeLimitSec, limit::Real)
    return MOI.set(
        model,
        MOI.RawOptimizerAttribute("time_limit"),
        Float64(limit),
    )
end

function MOI.get(model::Optimizer, ::MOI.TimeLimitSec)
    return MOI.get(model, MOI.RawOptimizerAttribute("time_limit"))
end

###
### MOI.Silent
###

MOI.supports(::Optimizer, ::MOI.Silent) = true

function MOI.get(model::Optimizer, ::MOI.Silent)
    return !MOI.get(model, MOI.RawOptimizerAttribute("output_flag"))
end

function MOI.set(model::Optimizer, ::MOI.Silent, flag::Bool)
    MOI.set(model, MOI.RawOptimizerAttribute("output_flag"), !flag)
    return
end

###
### MOI.Name
###

MOI.supports(::Optimizer, ::MOI.Name) = true

MOI.get(model::Optimizer, ::MOI.Name) = model.name

MOI.set(model::Optimizer, ::MOI.Name, name::String) = (model.name = name)

###
### Variables
###

function MOI.get(model::Optimizer, ::MOI.NumberOfVariables)
    return length(model.variable_info)
end

function MOI.get(model::Optimizer, ::MOI.ListOfVariableIndices)
    return sort!(collect(keys(model.variable_info)), by = x -> x.value)
end

function _info(model::Optimizer, key::MOI.VariableIndex)
    if haskey(model.variable_info, key)
        return model.variable_info[key]
    end
    return throw(MOI.InvalidIndex(key))
end

"""
    column(model::Optimizer, x::MOI.VariableIndex)

Return the 0-indexed column associated with `x` in `model`.
"""
column(model::Optimizer, x::MOI.VariableIndex) = _info(model, x).column

function MOI.add_variable(model::Optimizer)
    # Initialize `_VariableInfo` with a dummy `VariableIndex` and a column,
    # because we need `add_item` to tell us what the `VariableIndex` is.
    index = CleverDicts.add_item(
        model.variable_info,
        _VariableInfo(MOI.VariableIndex(0), Cint(0)),
    )
    info = _info(model, index)
    # Now, set `.index` and `.column`.
    info.index = index
    info.column = Cint(length(model.variable_info) - 1)
    ret = Highs_addCol(model, 0.0, -Inf, Inf, 0, C_NULL, C_NULL)
    _check_ret(ret)
    return index
end

function MOI.is_valid(model::Optimizer, v::MOI.VariableIndex)
    return haskey(model.variable_info, v)
end

function MOI.delete(model::Optimizer, v::MOI.VariableIndex)
    col = column(model, v)
    ret = Highs_deleteColsByRange(model, col, col)
    _check_ret(ret)
    delete!(model.variable_info, v)
    for other_info in values(model.variable_info)
        if other_info.column > col
            other_info.column -= 1
        end
    end
    model.name_to_variable = nothing
    model.name_to_constraint_index = nothing
    return
end

#
# Variable names
#

function MOI.supports(
    ::Optimizer,
    ::MOI.VariableName,
    ::Type{MOI.VariableIndex},
)
    return true
end

function MOI.get(model::Optimizer, ::Type{MOI.VariableIndex}, name::String)
    if model.name_to_variable === nothing
        _rebuild_name_to_variable(model)
    end
    if haskey(model.name_to_variable, name)
        variable = model.name_to_variable[name]
        if variable === nothing
            error("Duplicate name detected: $(name)")
        end
        return variable
    end
    return nothing
end

function MOI.get(model::Optimizer, ::MOI.VariableName, v::MOI.VariableIndex)
    return _info(model, v).name
end

function MOI.set(
    model::Optimizer,
    ::MOI.VariableName,
    v::MOI.VariableIndex,
    name::String,
)
    info = _info(model, v)
    info.name = name
    model.name_to_variable = nothing
    return
end

function _rebuild_name_to_variable(model::Optimizer)
    model.name_to_variable = Dict{String,Union{Nothing,MOI.VariableIndex}}()
    for (index, info) in model.variable_info
        if isempty(info.name)
            continue
        end
        if haskey(model.name_to_variable, info.name)
            model.name_to_variable[info.name] = nothing
        else
            model.name_to_variable[info.name] = index
        end
    end
    return
end

###
### Objectives
###

MOI.supports(::Optimizer, ::MOI.ObjectiveSense) = true

function MOI.set(
    model::Optimizer,
    ::MOI.ObjectiveSense,
    sense::MOI.OptimizationSense,
)
    x = sense == MOI.MAX_SENSE ? Cint(kMaximize) : Cint(kMinimize)
    ret = Highs_changeObjectiveSense(model, x)
    _check_ret(ret)
    if sense == MOI.FEASIBILITY_SENSE
        model.is_feasibility = true
        # TODO(odow): cache the mask.
        n = MOI.get(model, MOI.NumberOfVariables())
        ret =
            Highs_changeColsCostByMask(model, ones(Cint, n), zeros(Cdouble, n))
        _check_ret(ret)
        ret = Highs_changeObjectiveOffset(model, 0.0)
        _check_ret(ret)
    else
        model.is_feasibility = false
    end
    return
end

function MOI.get(model::Optimizer, ::MOI.ObjectiveSense)
    if model.is_feasibility
        return MOI.FEASIBILITY_SENSE
    end
    senseP = Ref{Cint}()
    ret = Highs_getObjectiveSense(model, senseP)
    _check_ret(ret)
    return senseP[] == Cint(kMinimize) ? MOI.MIN_SENSE : MOI.MAX_SENSE
end

function MOI.supports(
    ::Optimizer,
    ::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}},
)
    return true
end

function MOI.get(::Optimizer, ::MOI.ObjectiveFunctionType)
    return MOI.ScalarAffineFunction{Float64}
end

function MOI.set(
    model::Optimizer,
    ::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}},
    f::MOI.ScalarAffineFunction{Float64},
)
    num_vars = length(model.variable_info)
    obj = zeros(Float64, num_vars)
    for term in f.terms
        col = column(model, term.variable)
        obj[col+1] += term.coefficient
    end
    # TODO(odow): cache the mask.
    mask = ones(Cint, num_vars)
    ret = Highs_changeColsCostByMask(model, mask, obj)
    _check_ret(ret)
    ret = Highs_changeObjectiveOffset(model, f.constant)
    _check_ret(ret)
    return
end

function MOI.get(
    model::Optimizer,
    ::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}},
)
    ncols = MOI.get(model, MOI.NumberOfVariables())
    offset = Ref{Cdouble}(0.0)
    ret = Highs_getObjectiveOffset(model, offset)
    _check_ret(ret)
    if ncols == 0
        return MOI.ScalarAffineFunction{Float64}(
            MOI.ScalarAffineTerm{Float64}[],
            offset[],
        )
    end
    num_colsP, nnzP = Ref{Cint}(0), Ref{Cint}(0)
    costs = Vector{Cdouble}(undef, ncols)
    ret = Highs_getColsByRange(
        model,
        0,
        ncols - 1,
        num_colsP,
        costs,
        C_NULL,
        C_NULL,
        nnzP,
        C_NULL,
        C_NULL,
        C_NULL,
    )
    return MOI.ScalarAffineFunction{Float64}(
        MOI.ScalarAffineTerm{Float64}[
            MOI.ScalarAffineTerm(costs[info.column+1], index) for
            (index, info) in model.variable_info if
            !iszero(costs[info.column+1])
        ],
        offset[],
    )
end

function MOI.get(model::Optimizer, ::MOI.ObjectiveFunction{F}) where {F}
    obj = MOI.get(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
    )
    return convert(F, obj)
end

function MOI.modify(
    model::Optimizer,
    ::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}},
    chg::MOI.ScalarConstantChange{Float64},
)
    ret = Highs_changeObjectiveOffset(model, chg.new_constant)
    _check_ret(ret)
    model.is_objective_function_set = true
    return
end

function MOI.modify(
    model::Optimizer,
    ::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}},
    chg::MOI.ScalarCoefficientChange{Float64},
)
    ret = Highs_changeColCost(
        model,
        column(model, chg.variable),
        chg.new_coefficient,
    )
    _check_ret(ret)
    model.is_objective_function_set = true
    return
end

###
### VariableIndex-in-Set constraints.
###

function MOI.supports_constraint(
    ::Optimizer,
    ::Type{MOI.VariableIndex},
    ::Type{<:_SCALAR_SETS},
)
    return true
end

function _bound_enums(::Type{MOI.LessThan{Float64}})
    return (_BOUND_LESS_THAN, _BOUND_LESS_AND_GREATER_THAN)
end

function _bound_enums(::Type{MOI.GreaterThan{Float64}})
    return (_BOUND_GREATER_THAN, _BOUND_LESS_AND_GREATER_THAN)
end

_bound_enums(::Type{MOI.Interval{Float64}}) = (_BOUND_INTERVAL,)

_bound_enums(::Type{MOI.EqualTo{Float64}}) = (_BOUND_EQUAL_TO,)

function MOI.get(
    model::Optimizer,
    ::MOI.ListOfConstraintIndices{MOI.VariableIndex,S},
) where {S<:_SCALAR_SETS}
    indices = MOI.ConstraintIndex{MOI.VariableIndex,S}[
        MOI.ConstraintIndex{MOI.VariableIndex,S}(key.value) for
        (key, info) in model.variable_info if info.bound in _bound_enums(S)
    ]
    return sort!(indices, by = x -> x.value)
end

function _info(
    model::Optimizer,
    c::MOI.ConstraintIndex{MOI.VariableIndex,<:Any},
)
    var_index = MOI.VariableIndex(c.value)
    if haskey(model.variable_info, var_index)
        return _info(model, var_index)
    end
    return throw(MOI.InvalidIndex(c))
end

"""
    column(
        model::Optimizer,
        c::MOI.ConstraintIndex{MOI.VariableIndex,<:Any},
    )

Return the 0-indexed column associated with the variable bounds `c` in `model`.
"""
function column(
    model::Optimizer,
    c::MOI.ConstraintIndex{MOI.VariableIndex,<:Any},
)
    return _info(model, c).column
end

function MOI.is_valid(
    model::Optimizer,
    c::MOI.ConstraintIndex{MOI.VariableIndex,MOI.LessThan{Float64}},
)
    if haskey(model.variable_info, MOI.VariableIndex(c.value))
        info = _info(model, c)
        return info.bound == _BOUND_LESS_THAN ||
               info.bound == _BOUND_LESS_AND_GREATER_THAN
    end
    return false
end

function MOI.is_valid(
    model::Optimizer,
    c::MOI.ConstraintIndex{MOI.VariableIndex,MOI.GreaterThan{Float64}},
)
    if haskey(model.variable_info, MOI.VariableIndex(c.value))
        info = _info(model, c)
        return info.bound == _BOUND_GREATER_THAN ||
               info.bound == _BOUND_LESS_AND_GREATER_THAN
    end
    return false
end

function MOI.is_valid(
    model::Optimizer,
    c::MOI.ConstraintIndex{MOI.VariableIndex,MOI.Interval{Float64}},
)
    return haskey(model.variable_info, MOI.VariableIndex(c.value)) &&
           _info(model, c).bound == _BOUND_INTERVAL
end

function MOI.is_valid(
    model::Optimizer,
    c::MOI.ConstraintIndex{MOI.VariableIndex,MOI.EqualTo{Float64}},
)
    return haskey(model.variable_info, MOI.VariableIndex(c.value)) &&
           _info(model, c).bound == _BOUND_EQUAL_TO
end

function MOI.get(
    model::Optimizer,
    ::MOI.ConstraintFunction,
    c::MOI.ConstraintIndex{MOI.VariableIndex,<:Any},
)
    MOI.throw_if_not_valid(model, c)
    return MOI.VariableIndex(c.value)
end

function MOI.set(
    ::Optimizer,
    ::MOI.ConstraintFunction,
    ::MOI.ConstraintIndex{MOI.VariableIndex,<:Any},
    ::MOI.VariableIndex,
)
    return throw(MOI.SettingVariableIndexNotAllowed())
end

function _throw_if_existing_lower(
    info::_VariableInfo,
    ::S,
) where {S<:MOI.AbstractSet}
    if info.bound == _BOUND_LESS_AND_GREATER_THAN
        throw(MOI.LowerBoundAlreadySet{MOI.GreaterThan{Float64},S}(info.index))
    elseif info.bound == _BOUND_GREATER_THAN
        throw(MOI.LowerBoundAlreadySet{MOI.GreaterThan{Float64},S}(info.index))
    elseif info.bound == _BOUND_INTERVAL
        throw(MOI.LowerBoundAlreadySet{MOI.Interval{Float64},S}(info.index))
    elseif info.bound == _BOUND_EQUAL_TO
        throw(MOI.LowerBoundAlreadySet{MOI.EqualTo{Float64},S}(info.index))
    end
    return
end

function _throw_if_existing_upper(
    info::_VariableInfo,
    ::S,
) where {S<:MOI.AbstractSet}
    if info.bound == _BOUND_LESS_AND_GREATER_THAN
        throw(MOI.UpperBoundAlreadySet{MOI.LessThan{Float64},S}(info.index))
    elseif info.bound == _BOUND_LESS_THAN
        throw(MOI.UpperBoundAlreadySet{MOI.LessThan{Float64},S}(info.index))
    elseif info.bound == _BOUND_INTERVAL
        throw(MOI.UpperBoundAlreadySet{MOI.Interval{Float64},S}(info.index))
    elseif info.bound == _BOUND_EQUAL_TO
        throw(MOI.UpperBoundAlreadySet{MOI.EqualTo{Float64},S}(info.index))
    end
    return
end

function MOI.add_constraint(
    model::Optimizer,
    f::MOI.VariableIndex,
    s::S,
) where {S<:_SCALAR_SETS}
    info = _info(model, f)
    _update_info(info, s)
    index = MOI.ConstraintIndex{MOI.VariableIndex,typeof(s)}(f.value)
    MOI.set(model, MOI.ConstraintSet(), index, s)
    return index
end

function MOI.delete(
    model::Optimizer,
    c::MOI.ConstraintIndex{MOI.VariableIndex,MOI.LessThan{Float64}},
)
    MOI.throw_if_not_valid(model, c)
    info = _info(model, c)
    ret = Highs_changeColBounds(model, info.column, info.lower, Inf)
    _check_ret(ret)
    info.upper = Inf
    if info.bound == _BOUND_LESS_AND_GREATER_THAN
        info.bound = _BOUND_GREATER_THAN
    else
        info.bound = _BOUND_NONE
    end
    model.name_to_constraint_index = nothing
    return
end

function MOI.delete(
    model::Optimizer,
    c::MOI.ConstraintIndex{MOI.VariableIndex,MOI.GreaterThan{Float64}},
)
    MOI.throw_if_not_valid(model, c)
    info = _info(model, c)
    ret = Highs_changeColBounds(model, info.column, -Inf, info.upper)
    _check_ret(ret)
    info.lower = -Inf
    if info.bound == _BOUND_LESS_AND_GREATER_THAN
        info.bound = _BOUND_LESS_THAN
    else
        info.bound = _BOUND_NONE
    end
    model.name_to_constraint_index = nothing
    return
end

function MOI.delete(
    model::Optimizer,
    c::MOI.ConstraintIndex{MOI.VariableIndex,MOI.Interval{Float64}},
)
    MOI.throw_if_not_valid(model, c)
    info = _info(model, c)
    ret = Highs_changeColBounds(model, info.column, -Inf, Inf)
    _check_ret(ret)
    info.lower, info.upper = -Inf, Inf
    info.bound = _BOUND_NONE
    model.name_to_constraint_index = nothing
    return
end

function MOI.delete(
    model::Optimizer,
    c::MOI.ConstraintIndex{MOI.VariableIndex,MOI.EqualTo{Float64}},
)
    MOI.throw_if_not_valid(model, c)
    info = _info(model, c)
    ret = Highs_changeColBounds(model, info.column, -Inf, Inf)
    _check_ret(ret)
    info.lower, info.upper = -Inf, Inf
    info.bound = _BOUND_NONE
    model.name_to_constraint_index = nothing
    return
end

function MOI.get(
    model::Optimizer,
    ::MOI.ConstraintSet,
    c::MOI.ConstraintIndex{MOI.VariableIndex,MOI.GreaterThan{Float64}},
)
    MOI.throw_if_not_valid(model, c)
    return MOI.GreaterThan(_info(model, c).lower)
end

function MOI.get(
    model::Optimizer,
    ::MOI.ConstraintSet,
    c::MOI.ConstraintIndex{MOI.VariableIndex,MOI.LessThan{Float64}},
)
    MOI.throw_if_not_valid(model, c)
    return MOI.LessThan(_info(model, c).upper)
end

function MOI.get(
    model::Optimizer,
    ::MOI.ConstraintSet,
    c::MOI.ConstraintIndex{MOI.VariableIndex,MOI.EqualTo{Float64}},
)
    MOI.throw_if_not_valid(model, c)
    return MOI.EqualTo(_info(model, c).lower)
end

function MOI.get(
    model::Optimizer,
    ::MOI.ConstraintSet,
    c::MOI.ConstraintIndex{MOI.VariableIndex,MOI.Interval{Float64}},
)
    MOI.throw_if_not_valid(model, c)
    info = _info(model, c)
    return MOI.Interval(info.lower, info.upper)
end

function MOI.set(
    model::Optimizer,
    ::MOI.ConstraintSet,
    c::MOI.ConstraintIndex{MOI.VariableIndex,S},
    s::S,
) where {S<:_SCALAR_SETS}
    MOI.throw_if_not_valid(model, c)
    lower, upper = _bounds(s)
    info = _info(model, c)
    if S == MOI.LessThan{Float64}
        ret = Highs_changeColBounds(model, info.column, info.lower, upper)
        _check_ret(ret)
        info.upper = upper
    elseif S == MOI.GreaterThan{Float64}
        ret = Highs_changeColBounds(model, info.column, lower, info.upper)
        _check_ret(ret)
        info.lower = lower
    else
        ret = Highs_changeColBounds(model, info.column, lower, upper)
        _check_ret(ret)
        info.lower = lower
        info.upper = upper
    end
    return
end

function MOI.supports(
    ::Optimizer,
    ::MOI.ConstraintName,
    ::Type{MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64},S}},
) where {S}
    return true
end

###
### ScalarAffineFunction-in-Set
###

function MOI.supports_constraint(
    ::Optimizer,
    ::Type{MOI.ScalarAffineFunction{Float64}},
    ::Type{<:_SCALAR_SETS},
)
    return true
end

function MOI.get(
    model::Optimizer,
    ::MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{Float64},S},
) where {S<:_SCALAR_SETS}
    indices = MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64},S}[
        MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64},S}(key.value)
        for (key, info) in model.affine_constraint_info if _set(info) isa S
    ]
    return sort!(indices; by = x -> x.value)
end

function _info(
    model::Optimizer,
    c::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64},<:_SCALAR_SETS},
)
    key = _ConstraintKey(c.value)
    if haskey(model.affine_constraint_info, key)
        return model.affine_constraint_info[key]
    end
    return throw(MOI.InvalidIndex(c))
end

function row(
    model::Optimizer,
    c::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64},<:_SCALAR_SETS},
)
    return _info(model, c).row
end

function MOI.is_valid(
    model::Optimizer,
    c::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64},S},
) where {S<:_SCALAR_SETS}
    key = _ConstraintKey(c.value)
    info = get(model.affine_constraint_info, key, nothing)
    if info === nothing
        return false
    end
    return _set(info) isa S
end

function _indices_and_coefficients(
    indices::Vector{Cint},
    coefficients::Vector{Float64},
    model::Optimizer,
    f::MOI.ScalarAffineFunction{Float64},
)
    i = 1
    for term in f.terms
        indices[i] = column(model, term.variable)
        coefficients[i] = term.coefficient
        i += 1
    end
    return indices, coefficients
end

function _indices_and_coefficients(
    model::Optimizer,
    f::MOI.ScalarAffineFunction{Float64},
)
    f_canon = MOI.Utilities.canonical(f)
    nnz = length(f_canon.terms)
    indices, coefficients = zeros(Cint, nnz), zeros(Cdouble, nnz)
    _indices_and_coefficients(indices, coefficients, model, f_canon)
    return indices, coefficients
end

function MOI.add_constraint(
    model::Optimizer,
    f::MOI.ScalarAffineFunction{Float64},
    s::_SCALAR_SETS,
)
    if !iszero(f.constant)
        throw(
            MOI.ScalarFunctionConstantNotZero{Float64,typeof(f),typeof(s)}(
                f.constant,
            ),
        )
    end
    key = CleverDicts.add_item(model.affine_constraint_info, _ConstraintInfo(s))
    model.affine_constraint_info[key].row =
        Cint(length(model.affine_constraint_info) - 1)
    indices, coefficients = _indices_and_coefficients(model, f)
    lower, upper = _bounds(s)
    ret = Highs_addRow(
        model,
        lower,
        upper,
        length(indices),
        indices,
        coefficients,
    )
    _check_ret(ret)
    return MOI.ConstraintIndex{typeof(f),typeof(s)}(key.value)
end

function MOI.delete(
    model::Optimizer,
    c::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64},<:_SCALAR_SETS},
)
    r = row(model, c)
    ret = Highs_deleteRowsByRange(model, r, r)
    _check_ret(ret)
    for info in values(model.affine_constraint_info)
        if info.row > r
            info.row -= 1
        end
    end
    key = _ConstraintKey(c.value)
    delete!(model.affine_constraint_info, key)
    model.name_to_constraint_index = nothing
    return
end

function MOI.get(
    model::Optimizer,
    ::MOI.ConstraintSet,
    c::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64},S},
) where {S<:_SCALAR_SETS}
    return _set(_info(model, c))
end

function MOI.set(
    model::Optimizer,
    ::MOI.ConstraintSet,
    c::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64},S},
    s::S,
) where {S<:_SCALAR_SETS}
    lower, upper = _bounds(s)
    info = _info(model, c)
    ret = Highs_changeRowBounds(model, info.row, lower, upper)
    _check_ret(ret)
    info.lower, info.upper = lower, upper
    return
end

function MOI.get(
    model::Optimizer,
    ::MOI.ConstraintFunction,
    c::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64},S},
) where {S<:_SCALAR_SETS}
    r = row(model, c)
    num_row = Ref{Cint}(0)
    num_nz = Ref{Cint}(0)
    matrix_start = Ref{Cint}(0)
    lower, upper = Ref{Cdouble}(), Ref{Cdouble}()
    _ = Highs_getRowsByRange(
        model,
        r,
        r,
        num_row,
        lower,
        upper,
        num_nz,
        matrix_start,
        C_NULL,
        C_NULL,
    )
    matrix_index = Vector{Cint}(undef, num_nz[])
    matrix_value = Vector{Cdouble}(undef, num_nz[])
    ret = Highs_getRowsByRange(
        model,
        r,
        r,
        num_row,
        lower,
        upper,
        num_nz,
        matrix_start,
        matrix_index,
        matrix_value,
    )
    _check_ret(ret)
    return MOI.ScalarAffineFunction(
        MOI.ScalarAffineTerm{Float64}[
            MOI.ScalarAffineTerm(
                val,
                model.variable_info[CleverDicts.LinearIndex(col + 1)].index,
            ) for (col, val) in zip(matrix_index, matrix_value) if !iszero(val)
        ],
        0.0,
    )
end

function MOI.supports(
    ::Optimizer,
    ::MOI.ConstraintName,
    ::Type{
        <:MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64},<:_SCALAR_SETS},
    },
)
    return true
end

function MOI.get(
    model::Optimizer,
    ::MOI.ConstraintName,
    c::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64},<:Any},
)
    return _info(model, c).name
end

function MOI.set(
    model::Optimizer,
    ::MOI.ConstraintName,
    c::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64},<:Any},
    name::String,
)
    info = _info(model, c)
    info.name = name
    model.name_to_constraint_index = nothing
    return
end

function MOI.get(model::Optimizer, ::Type{MOI.ConstraintIndex}, name::String)
    if model.name_to_constraint_index === nothing
        _rebuild_name_to_constraint_index(model)
    end
    if haskey(model.name_to_constraint_index, name)
        constr = model.name_to_constraint_index[name]
        if constr === nothing
            error("Duplicate constraint name detected: $(name)")
        end
        return constr
    end
    return nothing
end

function MOI.get(
    model::Optimizer,
    C::Type{MOI.ConstraintIndex{F,S}},
    name::String,
) where {F,S}
    index = MOI.get(model, MOI.ConstraintIndex, name)
    if index isa C
        return index::MOI.ConstraintIndex{F,S}
    end
    return nothing
end

function _rebuild_name_to_constraint_index(model::Optimizer)
    model.name_to_constraint_index =
        Dict{String,Union{Nothing,MOI.ConstraintIndex}}()
    for (key, info) in model.affine_constraint_info
        if isempty(info.name)
            continue
        end
        S = typeof(_set(info))
        _set_name_to_constraint_index(
            model,
            info.name,
            MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64},S}(key.value),
        )
    end
    return
end

function _set_name_to_constraint_index(
    model::Optimizer,
    name::String,
    index::MOI.ConstraintIndex,
)
    if haskey(model.name_to_constraint_index, name)
        model.name_to_constraint_index[name] = nothing
    else
        model.name_to_constraint_index[name] = index
    end
    return
end

function MOI.modify(
    model::Optimizer,
    c::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64},S},
    chg::MOI.ScalarCoefficientChange{Float64},
) where {S<:_SCALAR_SETS}
    ret = Highs_changeCoeff(
        model,
        row(model, c),
        column(model, chg.variable),
        chg.new_coefficient,
    )
    _check_ret(ret)
    return
end

function MOI.set(
    model::Optimizer,
    ::MOI.ConstraintFunction,
    c::MOI.ConstraintIndex{F,S},
    f::F,
) where {F<:MOI.ScalarAffineFunction{Float64},S<:_SCALAR_SETS}
    if !iszero(f.constant)
        throw(MOI.ScalarFunctionConstantNotZero{Float64,F,S}(f.constant))
    end
    # This loop is a bit slow because we query the existing function first, and
    # then set each coefficient one-by-one. It's probably okay until someone
    # complains but it will require some upstream changes to introduce a faster
    # way of modifying a list of coefficients.
    old = MOI.get(model, MOI.ConstraintFunction(), c)
    terms = Dict(x.variable => 0.0 for x in old.terms)
    for term in f.terms
        terms[term.variable] = term.coefficient
    end
    r = row(model, c)
    for (k, v) in terms
        ret = Highs_changeCoeff(model, r, column(model, k), v)
        _check_ret(ret)
    end
    return
end

###
### Optimize methods.
###

"""
    _store_solution(model::Optimizer, ret::Cint)

Get the solution from a run of HiGHS.

`ret` is the `HighsStatus` enum {OK = 0, Warning, Error} from `Highs_run`.
"""
function _store_solution(model::Optimizer, ret::Cint)
    x = model.solution
    x.status = ret == 0 ? _OPTIMIZE_OK : _OPTIMIZE_ERRORED
    x.has_primal_solution = false
    x.has_dual_solution = false
    x.has_dual_ray = false
    x.has_primal_ray = false
    numCols = Highs_getNumCols(model)
    numRows = Highs_getNumRows(model)
    resize!(x.colvalue, numCols)
    resize!(x.coldual, numCols)
    resize!(x.colstatus, numCols)
    resize!(x.rowvalue, numRows)
    resize!(x.rowdual, numRows)
    resize!(x.rowstatus, numRows)
    status = HighsModelStatus(Highs_getModelStatus(model))
    statusP = Ref{Cint}()
    if status == kInfeasible
        ret = Highs_getDualRay(model, statusP, x.rowdual)
        _check_ret(ret)
        x.has_dual_ray = statusP[] == 1
    elseif status == kUnbounded
        ret = Highs_getPrimalRay(model, statusP, x.colvalue)
        _check_ret(ret)
        x.has_primal_ray = statusP[] == 1
    else
        Highs_getIntInfoValue(model, "primal_solution_status", statusP)
        x.has_primal_solution = statusP[] == 2
        Highs_getIntInfoValue(model, "dual_solution_status", statusP)
        x.has_dual_solution = statusP[] == 2
        if x.has_primal_solution
            Highs_getSolution(
                model,
                x.colvalue,
                x.coldual,
                x.rowvalue,
                x.rowdual,
            )
            Highs_getBasis(model, x.colstatus, x.rowstatus)
        end
    end
    return
end

function MOI.optimize!(model::Optimizer)
    for info in model.binaries
        Highs_changeColBounds(
            model,
            info.column,
            max(0.0, info.lower),
            min(1.0, info.upper),
        )
    end
    ret = Highs_run(model)
    _store_solution(model, ret)
    for info in model.binaries
        Highs_changeColBounds(model, info.column, info.lower, info.upper)
    end
    return
end

function MOI.get(model::Optimizer, ::MOI.TerminationStatus)
    if model.solution.status == _OPTIMIZE_NOT_CALLED
        return MOI.OPTIMIZE_NOT_CALLED
    elseif model.solution.status == _OPTIMIZE_ERRORED
        return MOI.OTHER_ERROR
    end
    status = HighsModelStatus(Highs_getModelStatus(model))
    if status == kNotset
        return MOI.OTHER_ERROR
    elseif status == kLoadError
        return MOI.OTHER_ERROR
    elseif status == kModelError
        return MOI.INVALID_MODEL
    elseif status == kPresolveError
        return MOI.OTHER_ERROR
    elseif status == kSolveError
        return MOI.OTHER_ERROR
    elseif status == kPostsolveError
        return MOI.OTHER_ERROR
    elseif status == kModelEmpty
        return MOI.INVALID_MODEL
    elseif status == kOptimal
        return MOI.OPTIMAL
    elseif status == kInfeasible
        return MOI.INFEASIBLE
    elseif status == kUnboundedOrInfeasible
        return MOI.INFEASIBLE_OR_UNBOUNDED
    elseif status == kUnbounded
        return MOI.DUAL_INFEASIBLE
    elseif status == kObjectiveBound
        return MOI.OBJECTIVE_LIMIT
    elseif status == kObjectiveTarget
        return MOI.OBJECTIVE_LIMIT
    elseif status == kTimeLimit
        return MOI.TIME_LIMIT
    elseif status == kIterationLimit
        return MOI.ITERATION_LIMIT
    else
        @assert status == kUnknown
        return MOI.OTHER_ERROR
    end
end

function MOI.get(model::Optimizer, ::MOI.ResultCount)
    x = model.solution
    if x.has_primal_solution || x.has_dual_ray || x.has_primal_ray
        return 1
    else
        return 0
    end
end

function MOI.get(model::Optimizer, ::MOI.RawStatusString)
    status = HighsModelStatus(Highs_getModelStatus(model))
    return string(status)
end

function MOI.get(model::Optimizer, attr::MOI.PrimalStatus)
    if attr.result_index != 1
        return MOI.NO_SOLUTION
    elseif model.solution.has_primal_solution
        return MOI.FEASIBLE_POINT
    elseif model.solution.has_primal_ray
        return MOI.INFEASIBILITY_CERTIFICATE
    end
    return MOI.NO_SOLUTION
end

function MOI.get(model::Optimizer, attr::MOI.DualStatus)
    if attr.result_index != 1
        return MOI.NO_SOLUTION
    elseif model.solution.has_dual_solution
        return MOI.FEASIBLE_POINT
    elseif model.solution.has_dual_ray
        return MOI.INFEASIBILITY_CERTIFICATE
    end
    return MOI.NO_SOLUTION
end

function MOI.get(model::Optimizer, attr::MOI.ObjectiveValue)
    MOI.check_result_index_bounds(model, attr)
    return Highs_getObjectiveValue(model)
end

function MOI.get(model::Optimizer, attr::MOI.DualObjectiveValue)
    MOI.check_result_index_bounds(model, attr)
    return MOI.Utilities.get_fallback(model, attr, Float64)
end

function MOI.get(model::Optimizer, ::MOI.ObjectiveBound)
    # TODO(odow): there is a bug in HiGHS where it reports incorrect values for
    # mip_dual_bound in the LP case. As a work-around, just return the most
    # optimistic of the primal and dual values.
    p = Ref{Cdouble}()
    Highs_getDoubleInfoValue(model, "mip_dual_bound", p)
    sense = _sense_corrector(model)
    primal = Highs_getObjectiveValue(model)
    return sense == -1 ? max(primal, -p[]) : min(primal, p[])
end

function MOI.get(model::Optimizer, ::MOI.SolveTimeSec)
    return Highs_getRunTime(model)
end

function MOI.get(model::Optimizer, ::MOI.SimplexIterations)
    return Highs_getSimplexIterationCount(model)
end

function MOI.get(
    model::Optimizer,
    attr::MOI.VariablePrimal,
    x::MOI.VariableIndex,
)
    MOI.check_result_index_bounds(model, attr)
    return model.solution.colvalue[column(model, x)+1]
end

function MOI.get(
    model::Optimizer,
    attr::MOI.ConstraintPrimal,
    c::MOI.ConstraintIndex{MOI.VariableIndex,<:_SCALAR_SETS},
)
    MOI.check_result_index_bounds(model, attr)
    return MOI.get(model, MOI.VariablePrimal(), MOI.VariableIndex(c.value))
end

function MOI.get(
    model::Optimizer,
    attr::MOI.ConstraintPrimal,
    c::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64},<:_SCALAR_SETS},
)
    MOI.check_result_index_bounds(model, attr)
    return model.solution.rowvalue[row(model, c)+1]
end

"""
    _sense_corrector(model::Optimizer)

Return `-1` if `MAX_SENSE`, and `1` is `MIN_SENSE`. Useful for correcting
solution attributes which are sense-dependent.
"""
function _sense_corrector(model::Optimizer)
    senseP = Ref{Cint}()
    ret = Highs_getObjectiveSense(model, senseP)
    _check_ret(ret)
    return senseP[]
end

"""
    _farkas_variable_dual(model::Optimizer, col::Cint)

Return a Farkas dual associated with the variable bounds of `col`. Given a dual
ray:

    ā * x = λ' * A * x <= λ' * b = -β + sum(āᵢ * Uᵢ | āᵢ < 0) + sum(āᵢ * Lᵢ | āᵢ > 0)

The Farkas dual of the variable is ā, and it applies to the upper bound if ā < 0,
and it applies to the lower bound if ā > 0.
"""
function _farkas_variable_dual(model::Optimizer, col::Cint)
    num_nz, num_cols = Ref{Cint}(0), Ref{Cint}(0)
    # TODO(odow): how does getColsByRangeWork???
    m = Highs_getNumRows(model)
    matrix_start = zeros(Cint, 2)
    matrix_index = Vector{Cint}(undef, m)
    matrix_value = Vector{Cdouble}(undef, m)
    ret = Highs_getColsByRange(
        model,
        col,
        col,
        num_cols,
        C_NULL,
        C_NULL,
        C_NULL,
        num_nz,
        matrix_start,
        matrix_index,
        matrix_value,
    )
    _check_ret(ret)
    dual = 0.0
    for i in 1:num_nz[]
        dual += -model.solution.rowdual[matrix_index[i]+1] * matrix_value[i]
    end
    return dual
end

"""
    _signed_dual(dual::Float64, ::Type{Set})

A heuristic for determining whether the dual of an interval constraint applies
to the lower or upper bound. It can be wrong by at most the solver's tolerance.
"""
_signed_dual(dual::Float64, ::Type{MOI.LessThan{Float64}}) = min(dual, 0.0)
_signed_dual(dual::Float64, ::Type{MOI.GreaterThan{Float64}}) = max(dual, 0.0)
_signed_dual(dual::Float64, ::Any) = dual

"""
    _signed_dual(dual::Float64, ::Type{Set}, status::HighsBasisStatus)

Determine whether the dual of an interval constraint applies to the lower or
upper bound using the basis status reported by HiGHS.
"""
function _signed_dual(
    dual::Float64,
    ::Type{MOI.LessThan{Float64}},
    status::HighsBasisStatus,
)
    return status == kUpper ? dual : 0.0
end

function _signed_dual(
    dual::Float64,
    ::Type{MOI.GreaterThan{Float64}},
    status::HighsBasisStatus,
)
    return status == kLower ? dual : 0.0
end

_signed_dual(dual::Float64, ::Any, ::HighsBasisStatus) = dual

function MOI.get(
    model::Optimizer,
    attr::MOI.ConstraintDual,
    c::MOI.ConstraintIndex{MOI.VariableIndex,S},
) where {S<:_SCALAR_SETS}
    MOI.check_result_index_bounds(model, attr)
    col = column(model, c)
    if model.solution.has_dual_ray[] == 1
        return _signed_dual(_farkas_variable_dual(model, col), S)
    end
    dual = _sense_corrector(model) * model.solution.coldual[col+1]
    stat = HighsBasisStatus(model.solution.colstatus[col+1])
    return _signed_dual(dual, S, stat)
end

function MOI.get(
    model::Optimizer,
    attr::MOI.ConstraintDual,
    c::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64},S},
) where {S<:_SCALAR_SETS}
    MOI.check_result_index_bounds(model, attr)
    r = row(model, c) + 1
    dual = model.solution.rowdual[r]
    if model.solution.has_dual_ray[] == 1
        return _signed_dual(dual, S)
    end
    stat = HighsBasisStatus(model.solution.rowstatus[r])
    return _signed_dual(_sense_corrector(model) * dual, S, stat)
end

###
### MOI.ConstraintBasisStatus
###

# HiGHS only reports a single basis status for each ranged constraint. Therefore
# if the status is kLower or kUpper, we must distinguish whether or not it is
# refering to the given constraint.
function _nonbasic_status(code::HighsBasisStatus, ::Type{<:MOI.Interval})
    return code == kLower ? MOI.NONBASIC_AT_LOWER : MOI.NONBASIC_AT_UPPER
end
function _nonbasic_status(code::HighsBasisStatus, ::Type{<:MOI.LessThan})
    return code == kLower ? MOI.BASIC : MOI.NONBASIC
end
function _nonbasic_status(code::HighsBasisStatus, ::Type{<:MOI.GreaterThan})
    return code == kLower ? MOI.NONBASIC : MOI.BASIC
end
_nonbasic_status(::HighsBasisStatus, ::Type{<:MOI.EqualTo}) = MOI.NONBASIC

function MOI.get(
    model::Optimizer,
    attr::MOI.ConstraintBasisStatus,
    c::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64},S},
) where {S<:_SCALAR_SETS}
    MOI.check_result_index_bounds(model, attr)
    stat = HighsBasisStatus(model.solution.rowstatus[row(model, c)+1])
    if stat == kLower
        return _nonbasic_status(stat, S)
    elseif stat == kBasic
        return MOI.BASIC
    elseif stat == kUpper
        return _nonbasic_status(stat, S)
    end
    @assert stat == kZero || stat == kNonbasic
    return MOI.NONBASIC
end

function MOI.get(
    model::Optimizer,
    attr::MOI.VariableBasisStatus,
    x::MOI.VariableIndex,
)
    MOI.check_result_index_bounds(model, attr)
    stat = HighsBasisStatus(model.solution.colstatus[column(model, x)+1])
    if stat == kLower
        return MOI.NONBASIC_AT_LOWER
    elseif stat == kBasic
        return MOI.BASIC
    elseif stat == kUpper
        return MOI.NONBASIC_AT_UPPER
    end
    @assert stat == kZero || stat == kNonbasic
    return MOI.NONBASIC
end

###
### Integrality
###

function MOI.supports_constraint(
    model::Optimizer,
    ::Type{MOI.VariableIndex},
    ::Type{<:Union{MOI.ZeroOne,MOI.Integer}},
)
    solver = MOI.get(model, MOI.RawOptimizerAttribute("solver"))
    return solver != "simplex" && solver != "ipm"
end

function MOI.is_valid(
    model::Optimizer,
    ci::MOI.ConstraintIndex{MOI.VariableIndex,MOI.Integer},
)
    return _info(model, ci).type == _TYPE_INTEGER
end

function MOI.is_valid(
    model::Optimizer,
    ci::MOI.ConstraintIndex{MOI.VariableIndex,MOI.ZeroOne},
)
    return _info(model, ci).type == _TYPE_BINARY
end

function MOI.get(
    model::Optimizer,
    ::MOI.ListOfConstraintIndices{MOI.VariableIndex,MOI.Integer},
)
    indices = MOI.ConstraintIndex{MOI.VariableIndex,MOI.Integer}[
        MOI.ConstraintIndex{MOI.VariableIndex,MOI.Integer}(key.value) for
        (key, info) in model.variable_info if info.type == _TYPE_INTEGER
    ]
    return sort!(indices, by = x -> x.value)
end

function MOI.get(
    model::Optimizer,
    ::MOI.ListOfConstraintIndices{MOI.VariableIndex,MOI.ZeroOne},
)
    indices = MOI.ConstraintIndex{MOI.VariableIndex,MOI.ZeroOne}[
        MOI.ConstraintIndex{MOI.VariableIndex,MOI.ZeroOne}(key.value) for
        (key, info) in model.variable_info if info.type == _TYPE_BINARY
    ]
    return sort!(indices, by = x -> x.value)
end

function MOI.add_constraint(
    model::Optimizer,
    f::MOI.VariableIndex,
    ::MOI.Integer,
)
    info = _info(model, f)
    info.type = _TYPE_INTEGER
    ci = MOI.ConstraintIndex{MOI.VariableIndex,MOI.Integer}(f.value)
    Highs_changeColIntegrality(model, info.column, kInteger)
    return ci
end

function MOI.add_constraint(
    model::Optimizer,
    f::MOI.VariableIndex,
    ::MOI.ZeroOne,
)
    info = _info(model, f)
    info.type = _TYPE_BINARY
    ci = MOI.ConstraintIndex{MOI.VariableIndex,MOI.ZeroOne}(f.value)
    Highs_changeColIntegrality(model, info.column, kInteger)
    push!(model.binaries, info)
    return ci
end

function MOI.get(
    model::Optimizer,
    ::MOI.ConstraintSet,
    ci::MOI.ConstraintIndex{MOI.VariableIndex,MOI.Integer},
)
    MOI.throw_if_not_valid(model, ci)
    return MOI.Integer()
end

function MOI.get(
    model::Optimizer,
    ::MOI.ConstraintSet,
    ci::MOI.ConstraintIndex{MOI.VariableIndex,MOI.ZeroOne},
)
    MOI.throw_if_not_valid(model, ci)
    return MOI.ZeroOne()
end

function MOI.delete(
    model::Optimizer,
    ci::MOI.ConstraintIndex{MOI.VariableIndex,MOI.Integer},
)
    MOI.throw_if_not_valid(model, ci)
    info = _info(model, ci)
    info.type = _TYPE_CONTINUOUS
    Highs_changeColIntegrality(model, info.column, kContinuous)
    return
end

function MOI.delete(
    model::Optimizer,
    ci::MOI.ConstraintIndex{MOI.VariableIndex,MOI.ZeroOne},
)
    MOI.throw_if_not_valid(model, ci)
    info = _info(model, ci)
    info.type = _TYPE_CONTINUOUS
    Highs_changeColIntegrality(model, info.column, kContinuous)
    delete!(model.binaries, info)
    return
end

###
### MOI.copy_to
###

function _add_bounds(::Vector{Float64}, ub, i, s::MOI.LessThan{Float64})
    ub[i] = s.upper
    return
end

function _add_bounds(lb, ::Vector{Float64}, i, s::MOI.GreaterThan{Float64})
    lb[i] = s.lower
    return
end

function _add_bounds(lb, ub, i, s::MOI.EqualTo{Float64})
    lb[i], ub[i] = s.value, s.value
    return
end

function _add_bounds(lb, ub, i, s::MOI.Interval{Float64})
    lb[i], ub[i] = s.lower, s.upper
    return
end

_constraints(src, F, S) = MOI.get(src, MOI.ListOfConstraintIndices{F,S}())

function _extract_bound_data(
    dest::Optimizer,
    src::MOI.ModelLike,
    mapping,
    collower::Vector{Float64},
    colupper::Vector{Float64},
    ::Type{S},
) where {S}
    for c_index in _constraints(src, MOI.VariableIndex, S)
        f = MOI.get(src, MOI.ConstraintFunction(), c_index)
        s = MOI.get(src, MOI.ConstraintSet(), c_index)
        new_f = mapping[f]
        info = _info(dest, new_f)
        _add_bounds(collower, colupper, info.column + 1, s)
        _update_info(info, s)
        mapping[c_index] = MOI.ConstraintIndex{MOI.VariableIndex,S}(new_f.value)
    end
    return
end

function _copy_to_columns(dest::Optimizer, src::MOI.ModelLike, mapping)
    x_src = MOI.get(src, MOI.ListOfVariableIndices())
    numcols = Cint(length(x_src))
    for (i, x) in enumerate(x_src)
        index = CleverDicts.add_item(
            dest.variable_info,
            _VariableInfo(MOI.VariableIndex(0), Cint(0)),
        )
        info = _info(dest, index)
        info.name = MOI.get(dest, MOI.VariableName(), x)
        info.index = index
        info.column = Cint(i - 1)
        mapping[x] = index
    end
    fobj =
        MOI.get(src, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
    c = fill(0.0, numcols)
    for term in fobj.terms
        i = mapping[term.variable].value
        c[i] += term.coefficient
    end
    return numcols, c, fobj.constant
end

_add_sizehint!(vec, n) = sizehint!(vec, length(vec) + n)

function _extract_row_data(
    dest::Optimizer,
    src::MOI.ModelLike,
    mapping,
    rowlower::Vector{Float64},
    rowupper::Vector{Float64},
    I::Vector{Cint},
    J::Vector{Cint},
    V::Vector{Float64},
    ::Type{S},
) where {S}
    row = length(I) == 0 ? 1 : I[end] + 1
    list = _constraints(src, MOI.ScalarAffineFunction{Float64}, S)
    numrows = length(list)
    _add_sizehint!(rowlower, numrows)
    _add_sizehint!(rowupper, numrows)
    n_terms = 0
    fs = Array{MOI.ScalarAffineFunction{Float64}}(undef, numrows)
    for (i, c_index) in enumerate(list)
        f = MOI.get(src, MOI.ConstraintFunction(), c_index)
        fs[i] = f
        set = MOI.get(src, MOI.ConstraintSet(), c_index)
        l, u = _bounds(set)
        push!(rowlower, l - f.constant)
        push!(rowupper, u - f.constant)
        n_terms += length(f.terms)
        key = CleverDicts.add_item(
            dest.affine_constraint_info,
            _ConstraintInfo(set),
        )
        dest.affine_constraint_info[key].row =
            Cint(length(dest.affine_constraint_info) - 1)
        mapping[c_index] =
            MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64},S}(key.value)
    end
    _add_sizehint!(I, n_terms)
    _add_sizehint!(J, n_terms)
    _add_sizehint!(V, n_terms)
    for f in fs
        for term in f.terms
            push!(I, row)
            push!(J, Cint(mapping[term.variable].value))
            push!(V, term.coefficient)
        end
        row += 1
    end
    return
end

function _check_input_data(dest::Optimizer, src::MOI.ModelLike)
    for (F, S) in MOI.get(src, MOI.ListOfConstraintTypesPresent())
        if !MOI.supports_constraint(dest, F, S)
            throw(
                MOI.UnsupportedConstraint{F,S}(
                    "HiGHS does not support constraints of type $F-in-$S.",
                ),
            )
        end
    end
    fobj_type = MOI.get(src, MOI.ObjectiveFunctionType())
    if !MOI.supports(dest, MOI.ObjectiveFunction{fobj_type}())
        throw(MOI.UnsupportedAttribute(MOI.ObjectiveFunction(fobj_type)))
    end
    return
end

MOI.supports_incremental_interface(::Optimizer) = true

function MOI.copy_to(dest::Optimizer, src::MOI.ModelLike)
    MOI.empty!(dest)
    _check_input_data(dest, src)
    mapping = MOI.Utilities.IndexMap()
    numcol, colcost, offset = _copy_to_columns(dest, src, mapping)
    collower, colupper = fill(-Inf, numcol), fill(Inf, numcol)
    rowlower, rowupper = Float64[], Float64[]
    I, J, V = Cint[], Cint[], Float64[]
    for S in (
        MOI.GreaterThan{Float64},
        MOI.LessThan{Float64},
        MOI.EqualTo{Float64},
        MOI.Interval{Float64},
    )
        _extract_bound_data(dest, src, mapping, collower, colupper, S)
        _extract_row_data(dest, src, mapping, rowlower, rowupper, I, J, V, S)
    end
    numrow = Cint(length(rowlower))
    A = SparseArrays.sparse(I, J, V, numrow, numcol)
    # Extract integrality constraints
    integrality = fill(kContinuous, numcol)
    for ci in _constraints(src, MOI.VariableIndex, MOI.ZeroOne)
        integrality[_info(dest, ci).column+1] = kInteger
        new_x = mapping[MOI.VariableIndex(ci.value)]
        mapping[ci] = typeof(ci)(new_x.value)
        push!(dest.binaries, _info(dest, mapping[ci]))
    end
    for ci in _constraints(src, MOI.VariableIndex, MOI.Integer)
        integrality[_info(dest, ci).column+1] = kInteger
        new_x = mapping[MOI.VariableIndex(ci.value)]
        mapping[ci] = typeof(ci)(new_x.value)
    end
    # Build the model
    is_max = MOI.get(src, MOI.ObjectiveSense()) == MOI.MAX_SENSE
    ret = Highs_passModel(
        dest,
        numcol,
        numrow,
        length(A.nzval),
        0,          # q_num_nz,
        kColwise,   # a_format,
        0,          # q_format,
        is_max ? kMaximize : kMinimize,
        offset,
        colcost,
        collower,
        colupper,
        rowlower,
        rowupper,
        A.colptr .- Cint(1),
        A.rowval .- Cint(1),
        A.nzval,
        Cint[0],   # qstart,
        Cint[],    # qindex,
        Cdouble[], # qvalue,
        integrality,
    )
    _check_ret(ret)
    return mapping
end


text:
module Production

##################Dependencies####################
using Revise

if  pwd() ∉ LOAD_PATH
  push!(LOAD_PATH,pwd())
end

#=if "$(pwd())\\data" ∉ LOAD_PATH
  push!(LOAD_PATH,"$(pwd())\\data")
end=#

using Finometrics #my personal module
using DataFrames, Distributions, Dates, LinearAlgebra, Statistics,
   Formatting, StatsBase, Serialization,  CSV, Random, Plots
  #, Measures, Gadfly, Fontconfig, Cairo NOTE: uncomment for graphs

############################Exported functions
export Firm, Problem, Simulation, graphscript!

#############################Constants

const TOL_CONSTRAINT = 1.0e-10
const N_TOL = 1.0e-10

const OUT_NAME = "test"
const OUT_PAPER = "C:\\Users\\Clinton\\Dropbox\\Apps\\Overleaf\\productionproject\\sub"
const OUT_PATH = "output"

const PENALTY = -10_000.0

##########################code files##########################
include("Numerical.jl")
include("Firm.jl")
include("Problem.jl")
include("Simulation.jl")
include("Graphscripts.jl")

end # module

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
module ReactionNetworks

using Random
export ReactionChannels, AbstractReactionNetwork, simulate

struct ReactionChannel
    X::Vector{Float64}
    function ReactionChannel()
        new(randexp(1))
    end
end

function Base.reset(r::ReactionChannel)
    empty!(r.X)
    push!(r.X, randexp())
    return nothing
end
function Base.iterate(r::ReactionChannel, i::Int=1)
    length(r.X)==i && push!(r.X, randexp())
    r.X[i], i+1
end
Base.IteratorSize(::Type{ReactionChannel}) = Base.IsInfinite()
Base.eltype(::Type{ReactionChannel}) = Float64

abstract type AbstractReactionNetwork{N, M} end
function stoichiometric_matrix(::AbstractReactionNetwork) end
function reaction_rates!(v::AbstractVector{Float64}, RN::AbstractReactionNetwork{N,M}, x::AbstractVector{Int}) where {N,M} where T
    RN(v, x)
    return nothing
end
function reaction_rates(RN::AbstractReactionNetwork{N,M}, x::AbstractVector{Int}) where {N,M} where T
    v = zeros(M)
    reaction_rates!(v, RN, x)
    return v
end

ReactionChannels{M} = NTuple{M, Union{Nothing,ReactionChannel}}
ReactionChannels{M}() where M = tuple(Iterators.repeated(nothing, M)...)

function simulate(
    RN::AbstractReactionNetwork{N,M},
    x0::AbstractVector{Int},
    tspan::Tuple{Float64,Float64},
    channels::ReactionChannels{M} = ReactionChannels{M}(),
) where {N,M}

    length(x0)==N || error("Incorrect dimensions")
    S = stoichiometric_matrix(RN)
    t0, t1 = tspan
    t0<t1 || error("Need t0<t1 in tspan")

    x = copy(x0)
    t = t0

    xsav = Vector{Int}()
    tsav = Vector{Float64}()
    csav = tuple(((isnothing(chn) ? ReactionChannel() : chn) for chn in channels)...)

    distance_generators = broadcast(Iterators.Stateful, csav)
    distance_to_next_jump = [first(dg) for dg in distance_generators]
    speeds = reaction_rates(RN, x)
    times_to_next_jump = distance_to_next_jump ./ speeds

    while true
        push!(tsav, t)
        append!(xsav, x)

        dt, jumper_idx = findmin(times_to_next_jump)
        t += dt
        t>t1 && break
        @. distance_to_next_jump -= dt*speeds
        distance_to_next_jump[jumper_idx] = first(distance_generators[jumper_idx])

        x .+= S[jumper_idx]

        reaction_rates!(speeds, RN, x)
        @. times_to_next_jump = distance_to_next_jump / speeds
    end

    return tsav, reshape(xsav,N,:), csav
end

include("EnzymeModels.jl")
export EnzymeModels

end # module


text:
#' ---
#' title : Example Case 3 - Year Planning
#' author : Andrew Rosemberg
#' date : 15th Feb 2019
#' ---

#' # Introduction

#' This an example of the HydroPowerModels package for solving a simple stochastic case with the following specifications:
#'    - 3 Buses
#'    - 3 Lines
#'    - 2 Generators
#'    - 1 Reservoir and Hydrogenerator
#'    - 3 Scenarios
#'    - 12 Stages

#' # Case

#' ## Importing package and optimizer
using GLPK
using HydroPowerModels

#' ## Initialization
#'+ results =  "hidden"
if !@isdefined plot_bool
    plot_bool = true
end
using Random
seed = 1221

#' ## Load Case Specifications

#' Data
alldata = HydroPowerModels.parse_folder(joinpath(WEAVE_ARGS[:testcases_dir],"case3"));

#' Plot power grid graph
if plot_bool == true
    Random.seed!(seed)
    HydroPowerModels.plot_grid(alldata[1],node_label=false)
end

#' Parameters
params = create_param(  stages = 12,
                        model_constructor_grid  = DCPPowerModel,
                        post_method             = PowerModels.post_opf,
                        optimizer               = with_optimizer(GLPK.Optimizer));

#' ## Build Model
#'+ results =  "hidden"
m = hydrothermaloperation(alldata, params);

#' ## Train
#'+ results =  "hidden"
Random.seed!(seed)
start_time = time()
HydroPowerModels.train(m,iteration_limit = 100,stopping_rules= [SDDP.Statistical(num_replications = 20,iteration_period=20)]);
end_time = time() - start_time

#' Termination Status and solve time (s)
(SDDP.termination_status(m.policygraph), end_time)

#' Bounds
if plot_bool == true
    HydroPowerModels.plot_bound(m)
end

#' ## Simulation
import Random
Random.seed!(seed)
results = HydroPowerModels.simulate(m, 100);
results

#' ## Testing Results
using Test
#' Bound
@test isapprox(SDDP.calculate_bound(m.policygraph), 11298.218756240378, atol=10)
#' Number of Simulations
@test length(results[:simulations]) == 100

#' ## Plot Aggregated Results
if plot_bool == true
    HydroPowerModels.plot_aggregated_results(results)
end

#' # Annex 1: Case Summary
if plot_bool == true
    PowerModels.print_summary(alldata[1]["powersystem"])
end

#' # Annex 2: Plot Results
if plot_bool == true
    HydroPowerModels.plotresults(results)
end


text:
include("WeightGeneration.jl")

using .WeightGeneration

using NCDatasets
using Distributed
using SharedArrays
using Formatting

s_map_file = "b.e11.B1850C5CN.f09_g16.005.pop.h.SST.100001-109912.nc"
d_map_file = "domain.ocn.gx3v7.120323.nc"

wgt_file = "wgt_gx1v6_to_gx3v7.nc"

const NNN_max = 9
const missing_value = 1e20

Dataset(s_map_file, "r") do ds
        global s_mask, s_Nx, s_Ny, s_N, s_lat, s_lon

        s_mask = replace(reshape(ds["REGION_MASK"][:], :), missing=>NaN)
        s_Nx = ds.dim["nlon"]
        s_Ny = ds.dim["nlat"]
        s_N  = s_Nx * s_Ny
        s_lon = replace(reshape(ds["TLONG"][:], :), missing=>NaN)
        s_lat = replace(reshape(ds["TLAT"][:], :), missing=>NaN)

        global gi_s = WeightGeneration.GridInfo(
            gc_lon = s_lon,
            gc_lat = s_lat,
            area   = copy(s_lon),
            mask   = s_mask,
            unit_of_angle = :deg,
            dims = [s_Nx, s_Ny],
        )

end

Dataset(d_map_file, "r") do ds
        global d_mask, d_Nx, d_Ny, d_N, d_lat, d_lon

        d_mask = replace(reshape(ds["mask"][:], :), missing=>NaN)
        d_Nx = ds.dim["ni"]
        d_Ny = ds.dim["nj"]
        d_N = d_Nx * d_Ny
        d_lon = replace(reshape(ds["xc"][:], :), missing=>NaN)
        d_lat = replace(reshape(ds["yc"][:], :), missing=>NaN)

        global gi_d = WeightGeneration.GridInfo(
            gc_lon = d_lon,
            gc_lat = d_lat,
            area   = copy(d_lon),
            mask   = d_mask,
            unit_of_angle = :deg,
            dims = [d_Nx, d_Ny],
        )
end

WeightGeneration.genWeight_NearestNeighbors(wgt_file, gi_s, gi_d, NNN_max)

max_stars_repo_name: meteorologytoday/CESM-diagnostic

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM