text:
"""
IntersectingBoundaries()

Returns
-------
IntVert::AbstractArray{Float64, 2}            Contains all the intersecting points between all the boundaries of the simplices. Each row represents some intersecting point. Matrix of size D-by-N, where D = total number of intersecting points. Each row is an intersecting point.
ConvexExpIntVert::AbstractArray{Float64, 2}   Dimension D-by-(2N + 2), where D = total number of intersecting points. (2N + 2) corresponds to the number of vertices in each simplex * 2. The first n+1 columns correspond to the convex expansion coefficients of the intersecting points in terms of the vertices generating 's₁'. The remaining n+1 to (2N+2) columns correspond to the convex expansion coefficients of the intersecting points in terms of the vertices generating 's₂'. The faces of each simplex are numbered according to the column labels.
"""
function IntersectionOfBoundaries(s₁::AbstractArray{Float64, 2}, s₂::AbstractArray{Float64, 2},
                                    convexexp1in2::AbstractArray{Float64, 2},
                                    convexexp2in1::AbstractArray{Float64, 2},
                                    ordered_vertices1::AbstractVector{Int},
                                    ordered_vertices2::AbstractVector{Int},
                                    num1in2::Int,
                                    num2in1::Int,
                                    Ncomm::Int,
                                    tolerance::Float64)

    n::Int = size(s₁, 1)
    IntVert = Vector{Float64}(0)
    ConvexExpIntVert = Vector{Float64}(0)
    Z = zeros(Float64, 2*n+2, 1)

    if num1in2 <= num2in1
        # The last one is excluded because it corresponds to the whole simplex
        Labels = collect(2^num1in2:(2^(n+1)-2))

        #The one vertex boundaries correspond to the labels 2.^(0:n)
        #but since we discard all labels up to 2^(Nin(1))-1, then the one
        #vertex boundaries we still need to discard correspond to the labels
        #2.^(Nin(1):n), which correspond to the vertices with indices Nin(1)+1:n+1
        OneIndexLabels = transpose(2 .^ (num1in2:n)) -  2^num1in2 + 1
        Labels[OneIndexLabels] .= 0

        # Now all the boundaries are generated by 2 or more vertices
        # Notice that the first component of BoundaryBinaryLabels corresponds to the label 2^(Nin(1))+1
        BBL₁ = Binary(Labels[findall(Labels .> 0)], n)
        StartingIndexBound2 = 2^num2in1 - 2^num1in2 + 1 - (num2in1 - num1in2)
        BBL₂ = view(BBL₁, StartingIndexBound2:size(BBL₁, 1), :)
    else
        # The last one is excluded because it corresponds to the whole simplex
        Labels = transpose(collect(2^num2in1:(2^(n+1) - 2)))

        #The one vertex boundaries correspond to the labels 2.^(0:n)
        #but since we discard all labels up to 2^(Nin(1))-1, then the one
        #vertex boundaries we still need to discard correspond to the labels
        #2.^(Nin(1):n), which correspond to the vertices with indices Nin(1)+1:n+1
        OneIndexLabels = transpose(2 .^ (num2in1:n)) -  2^num2in1 + 1
        Labels[OneIndexLabels'] = zeros(size(OneIndexLabels))
        Labels = vec(Labels[findall(Labels .> 0)])

        # Now all the boundaries are generated by 2 or more vertices. The first component of
        # BoundaryBinaryLabels corresponds to the label 2^(Nin(1))+1
        BBL₂ = Binary(Labels, n)
        StartingIndexBound1 = 2^num1in2 - 2^num2in1 + 1 - (num1in2 - num2in1)
        BBL₁ = view(BBL₂, StartingIndexBound1:size(BBL₂, 1), :)
    end

    N = size(BBL₁, 1) * size(BBL₂, 1)
    N2 = size(BBL₂, 1)

    Indices = 1:n+1

    @inbounds for i = 1:N
        index1 = ceil(Int, i/N2)
        index2 = i - (index1 - 1) * N2
        b1 = view(BBL₁, index1, :) # column vector
        b2 = view(BBL₂, index2, :)# column vector

        num_vert = sum(b1 + b2) # total number of vertices in both boundaries (including repeatitions)
        no_common_vert = true

        #b1[1:Ncomm] + b2[1:Ncomm]   contains, 0,1 or 2
        #If there is any vertex shared by the boundaries, a 2 will appear somewhere
        if Ncomm > 0
            no_common_vert = maximum(b1[1:Ncomm] + b2[1:Ncomm]) < 2
        end

        if num_vert <= n+2 && no_common_vert
            dim1 = sum(b1)
            dim2 = sum(b2)
            if dim1 >= dim2
                r = dim1
                s = dim2
                TargetVertices = view(s₂,:,:)

                refboundary = view(ordered_vertices1, Indices[findall(b1 .> 0)]) # column vec
                targetboundary = view(ordered_vertices2, Indices[findall(b2 .> 0)]) # column vec
                β = view(convexexp2in1, :, :)
                Γ = view(convexexp2in1, setdiff(1:n+1, refboundary), targetboundary)
                Rank= rank(Γ)
                Rank0 = rank([Γ;ones(1, s)])
                no_vanishing_column = minimum(maximum(abs.(Γ), 1))
                switch = 0
            else
                s = dim1
                r = dim2
                TargetVertices = view(s₁, :, :)
                refboundary = view(ordered_vertices2, Indices[findall(b2 .> 0)]) # column vec
                targetboundary = view(ordered_vertices1, Indices[findall(b1 .> 0)]) # column vec
                β = view(convexexp1in2, :, :)
                Γ = view(convexexp1in2, setdiff(1:n+1, refboundary), targetboundary)
                Rank= rank(Γ)
                Rank0 = rank([Γ; ones(1, s)])
                no_vanishing_column = minimum(maximum(abs.(Γ), 1))
                switch = 1
            end

            if Rank0-Rank == 1 && Rank == s-1 && no_vanishing_column > 0
                λ = QR(Γ,tolerance)

                α = [1 - ones(Int, 1, r-1) * view(β, refboundary[2:r], targetboundary);
                        view(β, refboundary[2:r], targetboundary)] * λ
                α[abs.(α) .<= tolerance] = 0

                #if some of these coefficients are negative the boundaries simply do
                #not intersect.
                #if some are zero, then the boundaries are not minimal and the intersecting point has already been
                #computed or will be computed as the intersection of minimal
                #boundaries.
                #This also rules out duplication of points
                if min(minimum(α), minimum(λ)) > 0 && sum(findall(α .> 0)) > 1
                    # Filtering out non minimal boundaries
                    newpoint = view(TargetVertices, :, targetboundary)*λ
                    Z .= 0.0
                    if switch == 0
                        Z[refboundary] .= view(α, :)
                        Z[targetboundary .+ (n + 1)] .= view(λ, :)
                    else
                        Z[refboundary .+ (n + 1)] .= view(α, :)
                        Z[targetboundary] .= view(λ, :)
                    end
                    append!(IntVert,newpoint)
                    append!(ConvexExpIntVert, Z)
                end
            end
        end
    end

    return reshape(IntVert, n, transpose(div(length(IntVert), n))),
           reshape(ConvexExpIntVert, 2*n+2, transpose(div(length(ConvexExpIntVert), 2*n+2)))
end


text:
using DataStructures
using Base.Test

##############
# DefaultDicts
##############

# construction
@test_throws ErrorException DefaultDict()
@test_throws ErrorException DefaultDict(String, Int)

if VERSION >= v"0.4.0-dev+980"
    @test typeof(DefaultDict(0.0, 1 => 1.0)) == DefaultDict{Int,Float64,Float64}
end

# empty dictionary
d = DefaultDict(Char, Int, 1)
@test length(d) == 0
@test isempty(d)
@test d['c'] == 1
@test !isempty(d)
empty!(d)
@test isempty(d)

# access, modification
@test (d['a'] += 1) == 2
@test 'a' in keys(d)
@test haskey(d, 'a')
@test get(d, 'b', 0) == 0
@test !('b' in keys(d))
@test !haskey(d, 'b')
@test pop!(d, 'a') == 2
@test isempty(d)

for c in 'a':'z'
    d[c] = c-'a'+1
end

@test d['z'] == 26
@test d['@'] == 1
@test length(d) == 27
delete!(d, '@')
@test length(d) == 26

for (k,v) in d
    @test v == k-'a'+1
end

@test sort(collect(keys(d))) == collect('a':'z')
@test sort(collect(values(d))) == collect(1:26)

# Starting from an existing dictionary
# Note: dictionary is copied upon construction
e = Dict([('a',1), ('b',3), ('c',5)])
f = DefaultDict(0, e)
@test f['d'] == 0
@test_throws KeyError e['d']
e['e'] = 9
@test e['e'] == 9
@test f['e'] == 0

s = similar(d)
@test typeof(s) == typeof(d)
@test s.d.default == d.d.default

#####################
# DefaultOrderedDicts
#####################

# construction
@test_throws ErrorException DefaultOrderedDict()
@test_throws ErrorException DefaultOrderedDict(String, Int)

# empty dictionary
d = DefaultOrderedDict(Char, Int, 1)
@test length(d) == 0
@test isempty(d)
@test d['c'] == 1
@test !isempty(d)
empty!(d)
@test isempty(d)

# access, modification
@test (d['a'] += 1) == 2
@test 'a' in keys(d)
@test haskey(d, 'a')
@test get(d, 'b', 0) == 0
@test !('b' in keys(d))
@test !haskey(d, 'b')
@test pop!(d, 'a') == 2
@test isempty(d)

for c in 'a':'z'
    d[c] = c-'a'+1
end

@test d['z'] == 26
@test d['@'] == 1
@test length(d) == 27
delete!(d, '@')
@test length(d) == 26

for (k,v) in d
    @test v == k-'a'+1
end

@test collect(keys(d)) == collect('a':'z')
@test collect(values(d)) == collect(1:26)

s = similar(d)
@test typeof(s) == typeof(d)
@test s.d.default == d.d.default

ext: jl


text:
function find_prev_tag(tag)
    project_root = (@__DIR__) |> dirname |> abspath
    tags = readlines(`git -C $project_root tag`)
    sort!(tags)
    idx = indexin([tag], tags)[1]
    if idx == nothing return "NO-PREV-TAG" end
    return get(tags, idx - 1, "NO-PREV-TAG")
end

# modify build-tmp.jl to only check correct version libs
function install_products_filter(build_file)
    prod_filter = raw"""products = filter(products) do prod
    endswith(prod.libnames[1], "$(VERSION.major)_$(VERSION.minor)")
end
"""
    lines = open(build_file) do io
        read(io, String) |> x -> split(x, "\n")
    end
    prod_in, prod_out, filter_written = false, false, false
    open(build_file, "w") do io
        for line in lines
            if occursin("products = [", line)
                prod_in = true
            end
            if prod_in && line == "]"
                prod_out = true
            end
            write(io, line * "\n")
            if prod_out && !filter_written
                write(io, prod_filter * "\n")
                filter_written = true
            end
        end
    end
end

function include_build_script(version_str, try_prev=false)
    build_script_url = "https://github.com/TuringLang/Libtask.jl/releases/download/v$(version_str)/build_LibtaskDylib.v$(version_str).jl"
    build_script = joinpath(@__DIR__, "tmp-build.jl")
    build_script = try download(build_script_url, build_script) catch end
    if build_script == nothing && try_prev # no such file
        version_str = find_prev_tag("v$version_str") |> strip |> (x) -> lstrip(x, ['v'])
        return include_build_script(version_str, false)
    end
    install_products_filter(build_script)
    include(build_script)
end

function get_version_str()
    path = joinpath(@__DIR__, "../Project.toml")
    version_reg = r"version\s*=\s*\"(.*)\""
    open(path) do file
        lines = readlines(file)
        for line in lines
            m = match(version_reg, line)
            if isa(m, RegexMatch) return m.captures[1] end
        end
    end
end

version_str = get_version_str() |> strip |> (x) -> lstrip(x, ['v'])
include_build_script(version_str, true)


text:
using JuMP
using Ipopt

function Solve_Linearized_Dist_SOC(buses, lines, generators, fp0, fq0, l0)
info("Starting Linearized DistSOC Model")

fp0 = vcat([0], fp0)
fq0 = vcat([0], fq0)
l0 = vcat([0], l0)

println(fp0)

bus_set = collect(keys(buses))
line_set = collect(keys(lines))
gen_set = collect(keys(generators))
gen_bus = []
for g in keys(generators)
    push!(gen_bus, generators[g].bus_idx)
end
root_bus = 0
for b in keys(buses)
    if buses[b].is_root
        root_bus = b
    end
end

lines_to = Dict()
for l in keys(lines)
    lines_to[lines[l].to_node] = lines[l]
end

v_root = 1

m = Model(solver=IpoptSolver())

@variable(m, v[bus_set] >= 0) #variable for voltage square
@variable(m, fp[bus_set]) #variable for active power flow
@variable(m, fq[bus_set]) # variable for reactive power flow
@variable(m, gp[bus_set]) #variable for active power generation
@variable(m, gq[bus_set]) #variable for reactive power generation
@variable(m, l[bus_set] >= 0) # variable for quadratic current

@objective(m, Min, sum((gp[b])*buses[b].generator.cost for b in gen_bus))

for b in gen_bus
# All buses with generation
    @constraint(m, gp[b] <= buses[b].generator.g_P_max) #Upper limit constraint for active power generation
    @constraint(m, gp[b] >= 0) #Lower limit constraint for active power generation
    @constraint(m, gq[b] <= buses[b].generator.g_Q_max) #Upper limit constraint for reactive power generation
    @constraint(m, gq[b] >= -buses[b].generator.g_Q_max) #Lower limit constraint for reactive power generation
end

for b in setdiff(bus_set, gen_bus)
# All buses without generation
    @constraint(m, gp[b] == 0)
    @constraint(m, gq[b] == 0)
end

for b in bus_set
# All buses
    @constraint(m, buses[b].d_P - gp[b] + sum((fp[k] + l[k]*lines_to[k].r) for k in buses[b].children) == fp[b])
    @constraint(m, buses[b].d_Q - gq[b] + sum((fq[k] + l[k]*lines_to[k].x) for k in buses[b].children) == fq[b])

    @constraint(m, v[b] <= buses[b].v_max)
    @constraint(m, v[b] >= buses[b].v_min)
end

for b in setdiff(bus_set, [root_bus])
# All buses without root
    b_ancestor = buses[b].ancestor[1]
    R = lines_to[b].r
    X = lines_to[b].x
    term_p = (fp0[b]+l0[b]*R)
    term_q = (fq0[b]+l0[b]*X)

    @constraint(m, -fp0[b]^2 - fq0[b]^2 + 2fp0[b]*fp[b] + 2fq0[b]*fq[b] <= lines_to[b].s_max^2)
    @constraint(m, term_p^2 + term_q^2 +
                     2*(fp[b] - fp0[b])*term_p +
                     2*(fq[b] - fq0[b])*term_q +
                     2*(l[b] - l0[b])*(term_p*R + term_q*X) <= lines_to[b].s_max^2)
    # Same as above but rearranged
    # @constraint(m, -1*(fp0[b] + l0[b]*R)^2 - (fq0[b] + l0[b]*X)^2 +
    #                 2*(fp0[b] + l0[b]*R)*fp[b] + 2*(fq0[b] + l0[b]*X)*fq[b] +
    #                 2*(l0[b]*R^2 + fp0[b]*R + l0[b]*X^2 + fq0[b]*X)*l[b] <= lines_to[b].s_max^2)

    @constraint(m, v[b_ancestor] - 2*(R*fp[b] + X*fq[b]) - l[b]*(R^2 + X^2) == v[b])
    @constraint(m, ((2*fp0[b])/l0[b]) * fp[b] + ((2*fq0[b])/l0[b]) * fq[b] - ((fp0[b]^2 + fq0[b]^2)/l0[b]^2) * l[b] <= v[b])
end

#Voltage and flow constraints for root node
@constraint(m, v[root_bus] == v_root)
@constraint(m, l[root_bus] == 0)
@constraint(m, fp[root_bus] == 0)
@constraint(m, fq[root_bus] == 0)

solve(m)

# Prepare Results
bus_results = DataFrame(bus = Any[], d_P = Any[], d_Q = Any[], g_P = Any[], g_Q = Any[], v_squared = Any[])
line_results = DataFrame(line = Any[], from = Any[], to = Any[], f_P = Any[], f_q = Any[], a_squared = Any[])
for b in bus_set
    res = [b, buses[b].d_P, buses[b].d_Q, getvalue(gp[b]), getvalue(gq[b]), getvalue(v[b])]
    push!(bus_results, res)

    # Calc square current on line
    a_res = 0
    fp_res = getvalue(fp[b])
    fq_res = getvalue(fq[b])
    if b != root_bus
        a_res = getvalue(l[b])
        lres = [lines_to[b].index, buses[b].ancestor[1], b, fp_res, fq_res, a_res]
        push!(line_results, lres)
    end
end

sort!(bus_results, cols=:bus)
sort!(line_results, cols=:to)

return getobjectivevalue(m), bus_results, line_results

end


text:
using DiagonalSBP
using SparseArrays: sparsevec, spzeros, sparse
using LinearAlgebra: I, eigvals
using Logging: @info
using Printf: @sprintf
using PGFPlotsX: @pgf, Plot, Axis, Table, pgfsave, SemiLogYAxis, LogLogAxis
using ForwardDiff: derivative

function noncharacteristic_operator(p, N, R)
  # Total number of points
  Nq = N + 1

  # Get the various derivative operators
  (D, S0, SN, HI, H, x) = DiagonalSBP.D2(p, N; xc = (0, 1))

  # Form the SBP inner product matrix
  A = -H * D + SN - S0

  L¹ = sparsevec([ 1], [1], Nq)'
  L² = sparsevec([Nq], [1], Nq)'

  α = -(1 - R) / (1 + R)

  # Solving the system in first order form with [u; v]
  Ouu = spzeros(Nq, Nq)
  Ouv = I
  Ovu = -A
  Ovv = α * L²' * L² + α * L¹' * L¹

  # Form the operator
  O = [Ouu Ouv;
       HI*Ovu HI*Ovv]

  return (O, x, H)
end

function characteristic_operator(p, N, R)
  # Total number of points
  Nq = N + 1

  # Get the various derivative operators
  (D, S0, SN, HI, H, x) = DiagonalSBP.D2(p, N; xc = (0, 1))
  B¹, B² = S0[1:1, :], SN[Nq:Nq, :]

  # Form the SBP inner product matrix
  A = -H * D + SN - S0

  L¹ = sparse([1], [ 1], [1], 1, Nq)
  L² = sparse([1], [Nq], [1], 1, Nq)
  n¹ = -1
  n² = +1

  θ_R = DiagonalSBP.D2_remainder_parameters(p).θ_R
  θ_H = DiagonalSBP.D2_remainder_parameters(p).θ_H
  Γ = N * (1 / θ_R + 1 / θ_H)

  # Solving the system in first order form with [u; v; ϕ¹; ϕ²]
  Ouu = spzeros(Nq, Nq)
  Ouv = I
  Ouϕ¹ = spzeros(Nq, 1)
  Ouϕ² = spzeros(Nq, 1)

  Ovu  = -A
  Ovv  = spzeros(Nq, Nq)
  Ovϕ¹ = spzeros(Nq, 1)
  Ovϕ² = spzeros(Nq, 1)

  Oϕ¹u  = spzeros(1, Nq)
  Oϕ¹v  = spzeros(1, Nq)
  Oϕ¹ϕ¹ = spzeros(1, 1)
  Oϕ¹ϕ² = spzeros(1, 1)

  Oϕ²u  = spzeros(1, Nq)
  Oϕ²v  = spzeros(1, Nq)
  Oϕ²ϕ¹ = spzeros(1, 1)
  Oϕ²ϕ² = spzeros(1, 1)

  Ovu  += ((1 - R) / 2) * L¹' * (n¹ * B¹ - Γ * L¹) + n¹ * B¹' * L¹
  Ovv  += -((1 - R) / 2) * L¹' * L¹
  Ovϕ¹ += ((1 - R) / 2) * Γ * L¹' - n¹ * B¹'

  Oϕ¹u  += ((1 + R) / 2) * (Γ * L¹ - n¹ * B¹)
  Oϕ¹v  += ((1 + R) / 2) * L¹
  Oϕ¹ϕ¹ += sparse([1], [1], [-Γ * (1 + R) / 2], 1, 1)

  Ovu  +=  ((1 - R) / 2) * L²' * (n² * B² - Γ * L²) + n² * B²' * L²
  Ovv  += -((1 - R) / 2) * L²' * L²
  Ovϕ² +=  ((1 - R) / 2) * Γ * L²' - n² * B²'

  Oϕ²u  += ((1 + R) / 2) * (Γ * L² - n² * B²)
  Oϕ²v  += ((1 + R) / 2) * L²
  Oϕ²ϕ² += sparse([1], [1], [-Γ * (1 + R) / 2], 1, 1)

  # Form the operator
  O = [    Ouu  Ouv  Ouϕ¹  Ouϕ²
       HI*[Ovu  Ovv  Ovϕ¹  Ovϕ²]
           Oϕ¹u Oϕ¹v Oϕ¹ϕ¹ Oϕ¹ϕ²
           Oϕ²u Oϕ²v Oϕ²ϕ¹ Oϕ²ϕ²]

  return (O, x, H)
end

function convergence_test(p, Ns, R, t_final = 1.0)
  @info "sbp order = $p"

  # Pulse center
  μ = 1 / 2

  # Pulse width
  w = 1 / 15

  # Analytic solution
  # f(x) = exp(-((x - μ) / w)^2) / 2
  f(x) = 0 < x < 1 ? sin(2 * pi * x)^6 : 0

  ue(x, t) = (f(x - t) + f(x + t)) + R * (f(2 - x - t) + f(-x + t))
  ∂t_ue(x, t) = derivative(t->ue(x, t), t)

  # storage for error
  err_NC  = zeros(length(Ns))
  err_C   = zeros(length(Ns))

  # loop over all N
  for (iter, N) = enumerate(Ns)
    @info @sprintf " level %d with N = %4d" iter N

    (NC, x, H) = noncharacteristic_operator(p, N, R)
    q0 = [ue.(x, 0); ∂t_ue.(x, 0)]

    qf = exp(t_final * Matrix(NC)) * q0

    Δu = ue.(x, t_final) - qf[1:N+1]
    Δv = ∂t_ue.(x, t_final) - qf[(N+2):(2N+2)]

    err_NC[iter] = sqrt(Δu' * H * Δu)
    if iter > 1
      rate = ((log(err_NC[iter-1]) - log(err_NC[iter])) /
              (log(Ns[iter]) - log(Ns[iter-1])))
      @info @sprintf """ Non-Characteristic
      error = %.2e
      rate  = %.2e""" err_NC[iter] rate
    else
      @info @sprintf """ Non-Characteristic
      error = %.2e""" err_NC[iter]
    end

    (C, x, H) = characteristic_operator(p, N, R)
    q0 = [ue.(x, 0); ∂t_ue.(x, 0); 0; 0]

    qf = exp(t_final * Matrix(C)) * q0

    Δu = ue.(x, t_final) - qf[1:N+1]
    Δv = ∂t_ue.(x, t_final) - qf[(N+2):(2N+2)]

    err_C[iter] = sqrt(Δu' * H * Δu)
    if iter > 1
      rate = ((log(err_C[iter-1]) - log(err_C[iter])) /
              (log(Ns[iter]) - log(Ns[iter-1])))
      @info @sprintf """ Characteristic
      error = %.2e
      rate  = %.2e""" err_C[iter] rate
    else
      @info @sprintf """ Characteristic
      error = %.2e""" err_C[iter]
    end
  end
  return err_C, err_NC
end

function make_spectrum_plots(p, N, R)
  h = 1 / N
  (NC, _...) = noncharacteristic_operator(p, N, R)
  (C, _...)  = characteristic_operator(p, N, R)

  ev_NC = h * eigvals(Matrix(NC))
  ev_C  = h * eigvals(Matrix(C))

  ymin, ymax, xmin, xmax = -2.5, 2.5, -3.35, 1e-1

  pts = nothing
  if !all(xmin .< real.(ev_NC) .< xmax)
    pts = ev_NC[findall(.!(xmin .< real.(ev_NC) .< xmax))]
    @assert length(pts) == 2
    @assert pts[1] ≈ pts[2]
    @assert abs(imag(pts[1])) < 1e-14
    ev_NC = [ev_NC[xmin .< real.(ev_NC) .< xmax]; -3.1]
  end
  @assert all(xmin .< real.(ev_C ) .< xmax)
  @assert all(ymin .< imag.(ev_NC) .< ymax)
  @assert all(ymin .< imag.(ev_C ) .< ymax)

  xmin = xmin - 0.2

  @pgf a = Axis(
            {
             ymin = ymin,
             ymax = ymax,
             xmin = xmin,
             xmax = xmax,
             legend_entries = {"non-characteristic", "characteristic"},
             title = "R = $R, N = $N",
             xlabel = "\$\\textrm{Re}(h \\lambda)\$",
             ylabel = "\$\\textrm{Im}(h \\lambda)\$",
            }
           )
  @pgf push!(a, Plot(
                     {
                      only_marks,
                      color = "red",
                      mark = "x",
                     },
                     Table(real.(ev_NC), imag.(ev_NC)))
            )
  @pgf push!(a, Plot(
                     {
                      only_marks,
                      color = "blue",
                      mark = "+",
                     },
                     Table(real.(ev_C), imag.(ev_C))),
            )
  if !isnothing(pts)
    @pgf push!(a,
               [@sprintf """
                \\node[anchor=south] () at (axis cs:-3.1,0){\$%.2e\$};
                \\node[anchor=west] (source) at (axis cs:-3.1,0){};
                \\node (destination) at (axis cs:-3.35,0){};
                \\draw[->](source)--(destination);
                """ real(pts[1])])
  end

  pgfsave("bc_spectra_$R.tex", a)
  pgfsave("bc_spectra_$R.pdf", a)
end

function compare_eigs(p, N)
  h = 1 / N
  Rs = -0.95:0.05:0.95
  real_C = zeros(length(Rs))
  real_NC = zeros(length(Rs))
  imag_C = zeros(length(Rs))
  imag_NC = zeros(length(Rs))
  for (iter, R) = enumerate(Rs)
    (NC, _...) = noncharacteristic_operator(p, N, R)
    (C, _...)  = characteristic_operator(p, N, R)

    ev_NC = h * eigvals(Matrix(NC))
    ev_C  = h * eigvals(Matrix(C))

    real_C[iter] = minimum(real.(ev_C))
    real_NC[iter] = minimum(real.(ev_NC))
    imag_C[iter] = maximum(imag.(ev_C))
    imag_NC[iter] = maximum(imag.(ev_NC))
  end

  @pgf a = SemiLogYAxis(
                {
                 xmin = -1,
                 xmax =  1,
                 xlabel = "\$R\$",
                 ylabel = "max\$(-\\textrm{Re}(h \\lambda)\$",
                 legend_entries = {"non-characteristic", "characteristic"},
                }
               )
  @pgf push!(a, Plot(
                     {
                      only_marks,
                      color = "red",
                      mark = "x",
                     },
                     Table(Rs, -real_NC)
                    )
            )
  @pgf push!(a, Plot(
                     {
                      only_marks,
                      color = "blue",
                      mark = "+",
                     },
                     Table(Rs, -real_C)
                    )
            )
  pgfsave("real_eig.tex", a)
  pgfsave("real_eig.pdf", a)
end

#=
for R in (-0.99, -0.9, 0, 0.9, 1.0)
  make_spectrum_plots(4, 50, R)
end
compare_eigs(4, 50)
=#

let
  Ns = 17 * 2 .^ (0:5)
  hs = 1 ./ Ns

  tri_x = [hs[end-1], hs[end], hs[end-1], hs[end-1]]
  tri_y = [hs[end-1], hs[end], hs[end  ], hs[end-1]] ./ hs[end-1]
  scale = Dict()

  scale[(-0.99, 2)] = 0.001
  scale[( 0.0 , 2)] = 0.0002
  scale[( 0.99, 2)] = 0.001

  scale[(-0.99, 4)] = 1e-6
  scale[( 0.0 , 4)] = 5e-7
  scale[( 0.99, 4)] = 2e-6

  scale[(-0.99, 6)] = 2e-7
  scale[( 0.0 , 6)] = 1e-7
  scale[( 0.99, 6)] = 1e-8
  for R in (-0.99, 0.0, 0.99)
    C2, NC2 = convergence_test(2, Ns, R, 0.9)
    C4, NC4 = convergence_test(4, Ns, R, 0.9)
    C6, NC6 = convergence_test(6, Ns, R, 0.9)

    @pgf a = LogLogAxis(
                        {
                         ymax = 1,
                         ymin = 1e-11,
                         xmin = hs[end]/2,
                         xmax = 2hs[1],
                         xlabel = "\$h\$",
                         ylabel = "\$\\|\\Delta u\\|_{H}\$",
                        }
                       )
    tc2 = 0.6*NC2[end-1] * tri_y.^2
    @pgf push!(a, Plot(
                       {
                        color = "black",
                       },
                       Table(tri_x, tc2)
                      )
              )
    tc4 = 0.6*NC4[end-1] * tri_y.^4
    @pgf push!(a, Plot(
                       {
                        color = "black",
                       },
                       Table(tri_x, 0.6*NC4[end-1] * tri_y.^4)
                      )
              )
    tc6 = 0.6*NC6[end-1] * tri_y.^6
    @pgf push!(a, Plot(
                       {
                        color = "black",
                       },
                       Table(tri_x, 0.6*NC6[end-1] * tri_y.^5)
                      )
              )
    #=
    @pgf push!(a,
              [@sprintf("""
               \\node[anchor=north] () at (axis cs:%e, %e){\$s = 2\$};
               \\node[anchor=north] () at (axis cs:%e, %e){\$s = 4\$};
               \\node[anchor=north] () at (axis cs:%e, %e){\$s = 5\$};
               """,
              (hs[end-1] + hs[end]) / 2, tc2[2],
              (hs[end-1] + hs[end]) / 2, tc4[2],
              (hs[end-1] + hs[end]) / 2, tc6[2],
               )])
    =#
    @pgf push!(a,
               [@sprintf("""
                \\node[anchor=west] () at (axis cs:%e, %e){\$2\$};
                \\node[anchor=west] () at (axis cs:%e, %e){\$4\$};
                \\node[anchor=west] () at (axis cs:%e, %e){\$5\$};
                """,
                hs[end-1], exp((log(tc2[1]) + log(tc2[2]))/2),
                hs[end-1], exp((log(tc4[1]) + log(tc4[2]))/2),
                hs[end-1], exp((log(tc6[1]) + log(tc6[2]))/2),
               )])

    @pgf push!(a, Plot(
                       {
                        color = "red",
                        mark = "o",
                       },
                       Table(hs, C2)
                      )
              )
    @pgf push!(a, Plot(
                       {
                        color = "red",
                        mark = "x",
                       },
                       Table(hs, NC2)
                      )
              )
    @pgf push!(a, Plot(
                       {
                        color = "blue",
                        mark = "o",
                       },
                       Table(hs, C4)
                      )
              )
    @pgf push!(a, Plot(
                       {
                        color = "blue",
                        mark = "x",
                       },
                       Table(hs, NC4)
                      )
              )
    @pgf push!(a, Plot(
                       {
                        color = "green",
                        mark = "o",
                       },
                       Table(hs, C6)
                      )
              )
    @pgf push!(a, Plot(
                       {
                        color = "green",
                        mark = "x",
                       },
                       Table(hs, NC6)
                      )
              )
    pgfsave("convergence_$R.tex", a)
    pgfsave("convergence_$R.pdf", a)
  end
end


text:
module TestSets

using SparseArrays
using Test
using VecMathOptInterface

const MOI = VecMathOptInterface
const MOIU = MOI.Utilities

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$(name)", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    return
end

function test_diagonal_element()
    k = 0
    for j in 1:10
        for i in 1:j
            k += 1
            @test MOIU.is_diagonal_vectorized_index(k) == (i == j)
        end
    end
    for side_dim in 1:10
        set = MOI.PositiveSemidefiniteConeTriangle(side_dim)
        vec_dim = MOI.dimension(set)
        @test MOIU.side_dimension_for_vectorized_dimension(vec_dim) == side_dim
    end
    return
end

function test_side_dimension()
    for side_dim in 1:10
        set = MOI.PositiveSemidefiniteConeTriangle(side_dim)
        vec_dim = MOI.dimension(set)
        @test MOIU.side_dimension_for_vectorized_dimension(vec_dim) == side_dim
    end
    return
end

function test_shifts()
    @test MOIU.supports_shift_constant(MOI.EqualTo{Int})
    @test MOIU.shift_constant(MOI.EqualTo(3), 1) == MOI.EqualTo(4)
    @test MOIU.supports_shift_constant(MOI.GreaterThan{Int})
    @test MOIU.shift_constant(MOI.GreaterThan(6), -1) == MOI.GreaterThan(5)
    @test MOIU.supports_shift_constant(MOI.LessThan{Int})
    @test MOIU.shift_constant(MOI.LessThan(2), 2) == MOI.LessThan(4)
    @test MOIU.supports_shift_constant(MOI.Interval{Int})
    @test MOIU.shift_constant(MOI.Interval(-2, 3), 1) == MOI.Interval(-1, 4)
    @test MOIU.supports_shift_constant(MOI.ZeroOne) == false
    @test_throws MethodError MOIU.shift_constant(MOI.ZeroOne(), 1.0)
    return
end

function test_set_dot()
    vec = zeros(6)
    @test MOIU.set_dot(vec, vec, MOI.SecondOrderCone(6)) == 0
    @test MOIU.set_dot(vec, vec, MOI.PositiveSemidefiniteConeTriangle(3)) == 0
    @test MOIU.set_dot(vec, vec, MOI.PositiveSemidefiniteConeSquare(3)) == 0
    vec[5] = 1
    @test MOIU.set_dot(vec, vec, MOI.SecondOrderCone(6)) == 1
    @test MOIU.set_dot(vec, vec, MOI.PositiveSemidefiniteConeTriangle(3)) == 2
    @test MOIU.set_dot(vec, vec, MOI.PositiveSemidefiniteConeSquare(3)) == 1
    vec[5] = 0
    vec[3] = 1
    @test MOIU.set_dot(vec, vec, MOI.SecondOrderCone(6)) == 1
    @test MOIU.set_dot(vec, vec, MOI.PositiveSemidefiniteConeTriangle(3)) == 1
    @test MOIU.set_dot(vec, vec, MOI.PositiveSemidefiniteConeSquare(3)) == 1

    vec = zeros(7)
    @test MOIU.set_dot(vec, vec, MOI.RootDetConeTriangle(3)) == 0
    vec[5] = 1
    @test MOIU.set_dot(vec, vec, MOI.RootDetConeTriangle(3)) == 2
    vec = zeros(8)
    @test MOIU.set_dot(vec, vec, MOI.LogDetConeTriangle(3)) == 0
    vec[5] = 1
    @test MOIU.set_dot(vec, vec, MOI.LogDetConeTriangle(3)) == 1

    sp_vec = spzeros(6)
    @test MOIU.set_dot(sp_vec, sp_vec, MOI.SecondOrderCone(6)) == 0
    @test MOIU.set_dot(
        sp_vec,
        sp_vec,
        MOI.PositiveSemidefiniteConeTriangle(3),
    ) == 0
    @test MOIU.set_dot(sp_vec, sp_vec, MOI.PositiveSemidefiniteConeSquare(3)) ==
          0
    sp_vec[5] = 1
    @test MOIU.set_dot(sp_vec, sp_vec, MOI.SecondOrderCone(6)) == 1
    @test MOIU.set_dot(
        sp_vec,
        sp_vec,
        MOI.PositiveSemidefiniteConeTriangle(3),
    ) == 2
    @test MOIU.set_dot(sp_vec, sp_vec, MOI.PositiveSemidefiniteConeSquare(3)) ==
          1
    sp_vec[5] = 0
    sp_vec[3] = 1
    @test MOIU.set_dot(sp_vec, sp_vec, MOI.SecondOrderCone(6)) == 1
    @test MOIU.set_dot(
        sp_vec,
        sp_vec,
        MOI.PositiveSemidefiniteConeTriangle(3),
    ) == 1
    @test MOIU.set_dot(sp_vec, sp_vec, MOI.PositiveSemidefiniteConeSquare(3)) ==
          1

    sp_vec = spzeros(7)
    @test MOIU.set_dot(sp_vec, sp_vec, MOI.RootDetConeTriangle(3)) == 0
    sp_vec[5] = 1
    @test MOIU.set_dot(sp_vec, sp_vec, MOI.RootDetConeTriangle(3)) == 2
    sp_vec = spzeros(8)
    @test MOIU.set_dot(sp_vec, sp_vec, MOI.LogDetConeTriangle(3)) == 0
    sp_vec[5] = 1
    @test MOIU.set_dot(sp_vec, sp_vec, MOI.LogDetConeTriangle(3)) == 1
end

function test_dot_coefficients()
    vec = zeros(6)
    @test MOIU.dot_coefficients(vec, MOI.SecondOrderCone(6)) == vec
    @test MOIU.dot_coefficients(vec, MOI.PositiveSemidefiniteConeTriangle(3)) ==
          vec
    vec[5] = 1
    @test MOIU.dot_coefficients(vec, MOI.SecondOrderCone(6)) == vec
    @test MOIU.dot_coefficients(vec, MOI.PositiveSemidefiniteConeTriangle(3)) ==
          vec ./ 2
    vec[5] = 0
    vec[3] = 1
    @test MOIU.dot_coefficients(vec, MOI.SecondOrderCone(6)) == vec
    @test MOIU.dot_coefficients(vec, MOI.PositiveSemidefiniteConeTriangle(3)) ==
          vec

    vec = zeros(7)
    @test MOIU.dot_coefficients(vec, MOI.RootDetConeTriangle(3)) == vec
    vec[5] = 1
    @test MOIU.dot_coefficients(vec, MOI.RootDetConeTriangle(3)) == vec ./ 2
    vec = zeros(8)
    @test MOIU.dot_coefficients(vec, MOI.LogDetConeTriangle(3)) == vec
    vec[5] = 1
    @test MOIU.dot_coefficients(vec, MOI.LogDetConeTriangle(3)) == vec

    sp_vec = spzeros(6)
    @test MOIU.dot_coefficients(sp_vec, MOI.SecondOrderCone(6)) == sp_vec
    @test MOIU.dot_coefficients(
        sp_vec,
        MOI.PositiveSemidefiniteConeTriangle(3),
    ) == sp_vec
    sp_vec[5] = 1
    @test MOIU.dot_coefficients(sp_vec, MOI.SecondOrderCone(6)) == sp_vec
    @test MOIU.dot_coefficients(
        sp_vec,
        MOI.PositiveSemidefiniteConeTriangle(3),
    ) == sp_vec ./ 2
    sp_vec[5] = 0
    sp_vec[3] = 1
    @test MOIU.dot_coefficients(sp_vec, MOI.SecondOrderCone(6)) == sp_vec
    @test MOIU.dot_coefficients(
        sp_vec,
        MOI.PositiveSemidefiniteConeTriangle(3),
    ) == sp_vec

    sp_vec = spzeros(7)
    @test MOIU.dot_coefficients(sp_vec, MOI.RootDetConeTriangle(3)) == sp_vec
    sp_vec[5] = 1
    @test MOIU.dot_coefficients(sp_vec, MOI.RootDetConeTriangle(3)) ==
          sp_vec ./ 2
    sp_vec = spzeros(8)
    @test MOIU.dot_coefficients(sp_vec, MOI.LogDetConeTriangle(3)) == sp_vec
    sp_vec[5] = 1
    @test MOIU.dot_coefficients(sp_vec, MOI.LogDetConeTriangle(3)) == sp_vec
end

function test_trimap()
    @test MOIU.trimap(1, 1) == 1
    @test MOIU.trimap(1, 2) == 2
    @test MOIU.trimap(2, 1) == 2
    @test MOIU.trimap(2, 2) == 3
    @test MOIU.trimap(3, 1) == 4
    @test MOIU.trimap(1, 3) == 4
    @test MOIU.trimap(3, 2) == 5
    @test MOIU.trimap(2, 3) == 5
    @test MOIU.trimap(3, 3) == 6
    return
end

end  # module

TestSets.runtests()


text:
@testset "dual_sets.jl" begin
    @test Dualization._dual_set_type(MOI.Integer) === nothing
    @test Dualization._dual_set_type(MOI.ZeroOne) === nothing

    greater_than_Float64 = MOI.GreaterThan(0.0)
    greater_than_Int64 = MOI.GreaterThan(0)
    less_than_Float64 = MOI.LessThan(0.0)
    less_than_Int64 = MOI.LessThan(0)
    nonpositives_3 = MOI.Nonpositives(3)
    nonnegatives_3 = MOI.Nonnegatives(3)
    zeros_3 = MOI.Zeros(3)
    nonpositives_4 = MOI.Nonpositives(4)
    nonnegatives_4 = MOI.Nonnegatives(4)
    zeros_4 = MOI.Zeros(4)

    # GreaterThan
    @test Dualization._dual_set(MOI.GreaterThan(0.0)) == greater_than_Float64
    @test Dualization._dual_set(MOI.GreaterThan(1.0)) == greater_than_Float64
    @test Dualization._dual_set(MOI.GreaterThan(-1.0)) == greater_than_Float64
    @test Dualization._dual_set(MOI.GreaterThan(0)) == greater_than_Int64
    @test Dualization._dual_set(MOI.GreaterThan(1)) == greater_than_Int64
    @test Dualization._dual_set(MOI.GreaterThan(-1)) == greater_than_Int64

    # LessThan
    @test Dualization._dual_set(MOI.LessThan(0.0)) == less_than_Float64
    @test Dualization._dual_set(MOI.LessThan(1.0)) == less_than_Float64
    @test Dualization._dual_set(MOI.LessThan(-1.0)) == less_than_Float64
    @test Dualization._dual_set(MOI.LessThan(0)) == less_than_Int64
    @test Dualization._dual_set(MOI.LessThan(1)) == less_than_Int64
    @test Dualization._dual_set(MOI.LessThan(-1)) == less_than_Int64

    # EqualTo
    @test Dualization._dual_set(MOI.EqualTo(0.0)) == nothing
    @test Dualization._dual_set(MOI.EqualTo(1.0)) == nothing
    @test Dualization._dual_set(MOI.EqualTo(-1.0)) == nothing
    @test Dualization._dual_set(MOI.EqualTo(0)) == nothing
    @test Dualization._dual_set(MOI.EqualTo(1)) == nothing
    @test Dualization._dual_set(MOI.EqualTo(-1)) == nothing

    # Nonpositives
    @test Dualization._dual_set(nonpositives_3) == nonpositives_3
    @test Dualization._dual_set(nonpositives_4) == nonpositives_4

    # Nonnegatives
    @test Dualization._dual_set(nonnegatives_3) == nonnegatives_3
    @test Dualization._dual_set(nonnegatives_4) == nonnegatives_4

    # Zeros
    @test Dualization._dual_set(zeros_3) == MOI.Reals(3)
    @test Dualization._dual_set(zeros_4) == MOI.Reals(4)

    #SOC
    soc = MOI.SecondOrderCone(2)
    soc3 = MOI.SecondOrderCone(3)
    @test Dualization._dual_set(soc) == soc
    @test Dualization._dual_set(soc) != soc3
    @test Dualization._dual_set(soc3) == soc3

    #RSOC
    rsoc = MOI.RotatedSecondOrderCone(2)
    rsoc3 = MOI.RotatedSecondOrderCone(3)
    @test Dualization._dual_set(rsoc) == rsoc
    @test Dualization._dual_set(rsoc) != rsoc3
    @test Dualization._dual_set(rsoc3) == rsoc3

    #PSD
    @test Dualization._dual_set_type(MOI.PositiveSemidefiniteConeTriangle) == MOI.PositiveSemidefiniteConeTriangle
    psd = MOI.PositiveSemidefiniteConeTriangle(2)
    psd3 = MOI.PositiveSemidefiniteConeTriangle(3)
    @test Dualization._dual_set(psd) == psd
    @test Dualization._dual_set(psd) != psd3
    @test Dualization._dual_set(psd3) == psd3

    # Exponential
    exp = MOI.ExponentialCone()
    dual_exp = MOI.DualExponentialCone()
    @test Dualization._dual_set(exp) == dual_exp
    @test Dualization._dual_set(exp) != exp
    @test Dualization._dual_set(dual_exp) == exp
    @test Dualization._dual_set(dual_exp) != dual_exp

    # Power
    pow = MOI.PowerCone(0.3)
    pow4 = MOI.PowerCone(0.4)
    dual_pow = MOI.DualPowerCone(0.3)
    @test Dualization._dual_set(pow) == dual_pow
    @test Dualization._dual_set(pow) != pow
    @test Dualization._dual_set(dual_pow) == pow
    @test Dualization._dual_set(dual_pow) != pow4
    @test Dualization._dual_set(dual_pow) != dual_pow
end


text:
#
# This is an example of estimation of a portfolio model
# as in Hansen and Singleton, 1982. The data comes from
# Tauchen, JBES, 1986 (available by ftp).
# The notation follows my lecture notes. In particular
# gamma is the coefficient of RRA. Try out different sets
# of intruments and different lag lengths. You will see
# that beta is stable across  models, but that gamma is not.
#
# Michael Creel, Oct. 24, 2000
# revised 7 Nov. 2002
# translated to Octave 9/9/2003
# translated to Julia 07 Aug. 2017
# michael.creel@uab.es

# The main thing you have to to is define the moments for estimation,
# as in the function that follows immediately, and ensure that the
# efficient weight matrix is estimated appropriately. The rest
# is already programmed.

# the following function defines the moment conditions, given the parameter
# vector and the instruments. It returns the individual contributions
# to allow estimation of the covariance matrix of the moments.
using DelimitedFiles, Statistics, Econometrics
function portfolio_moments(theta, data)
	# parameters
	beta = theta[1]
	gam = theta[2]
	# data items
	c = data[:,1]
	r = data[:,2]
    inst = data[:,3:end]
	#  form error function
	# note that c = c_t / c_t-1 (for stationarity) was done in data preparation
	e = 1.0 .- beta*(1.0 .+ r) .* (c .^ (-gam))
	# cross with instruments
	m = e.*inst
end

function main()
data = readdlm("tauchen.data")
c = data[:,1]
p = data[:,2]
d = data[:,3]
# form net return and stationary consumption
r = (p + d) ./ lag(p,1) .- 1.0
c = c ./ lag(c,1); # ensure stationarity
# choose maximal lag of instruments.
max_lag = 1
inst = [c r d p]
inst = lags(inst, max_lag)
inst, junk, junk = stnorm(inst)
inst = [ones(size(inst,1),1) inst]
# drop rows with missing values
data = [c r inst]
data = data[max_lag+1:end,:]
# do estimation
moments = θ -> portfolio_moments(θ, data)
weight = 1.0
names = ["beta","gamma"]
# initial consistent estimate
theta = [0.9, 1]
thetahat, obj_value, D, ms, convergence = gmm(moments, theta, weight)
# second step with efficient weight matrix
omega = 50.0*cov(ms)
weight = inv(omega)
title = "Two step GMM estimation of portfolio model"
results = gmmresults(moments, thetahat, weight, title, names)
# CUE GMM
title = "CUE GMM estimation of portfolio model"
gmmresults(moments, results[1], "", title, names)
return
end
main()

max_stars_repo_name: nilshg/Econometrics

max_issues_repo_name: nilshg/Econometrics

max_forks_repo_name: nilshg/Econometrics


text:
cd(@__DIR__)

function normalMean(x)
    N    = length(x)
    sumx = 0.0
    for i in 1:N
       sumx += x[i]
    end
    return sumx/N
end

function onlineMean(oldMean,oldN,newX)
     oldSum  = oldMean * oldN
     newN    = oldN+1
     newMean = (oldSum + newX) / newN
     return (newMean, newN)
end

X = [1,2,3,4,5]

nMean = normalMean(X)

function computeMeanInLargeFile(filename)
    currentMean = 0.0
    currentN    = 0.0
    open(filename,"r") do f
        for ln in eachline(f)
            if ln != ""
                (currentMean, currentN) =  onlineMean(currentMean,currentN,parse(Float64,ln))
                println("Current mean: $currentMean")
            end
        end
    end

    return currentMean
end

olMean = computeMeanInLargeFile("x.txt")

max_stars_repo_name: sylvaticus/IntroSPMLJuliaCourse

max_issues_repo_name: sylvaticus/IntroSPMLJuliaCourse

max_forks_repo_name: sylvaticus/IntroSPMLJuliaCourse


text:
# Contains code helpful for the generation of transients

function τₛ(model::TcspcFitModel)
    # Returns lifetime values with amplitudes as a (2 * N) x 3 array

    τₘ = zeros(2 * nlcomp(model),3)
    τₛ!(τₘ,model)
    return τₘ
end

function τₛ!(τₘ::AbstractMatrix{<:AbstractFloat},
    model::TcspcFitModel)
    τᵢ = 1
    for τ ∈ model.τ
        τₘ[τᵢ,:] = [τ.α.val,τ.α.fixed,τ.α.glob]
        τₘ[τᵢ + 1,:] = [τ.τ.val,τ.τ.fixed,τ.τ.glob]
        τᵢ += 2
    end
end

function θₛ(model::TcspcFitModel)
    # Returns anisotropy values with amplitudes as a (2 * N * M) x 6 array

    θₘ = zeros(2 * (nacomp(model) * nlcomp(model)),6)
    θₛ!(θₘ,model)
    return θₘ
end

function θₛ!(θₘ::AbstractMatrix{<:AbstractFloat},
    model::TcspcFitModel)
    θᵢ = 1
    for θ ∈ model.θ
        for τ ∈ model.τ
            θₘ[θᵢ,:] = vcat([τ.α.val,τ.α.fixed,τ.α.glob],
                [θ.α.val,θ.α.fixed,θ.α.glob])
            θₘ[θᵢ + 1,:] = vcat([τ.τ.val,τ.τ.fixed,τ.τ.glob],
                [θ.τ.val,θ.τ.fixed,θ.τ.glob])
            θᵢ += 2
        end
    end
end

function taumat(model::TcspcFitModel,tcal::AbstractFloat,
    nbins::Integer)
    # Column for each lifetime in model
    t = -t_vector(tcal,nbins)
    nlcomp = length(model.τ)
    A = zeros(nbins,nlcomp)
    τₘ = τₛ(model)
    taumat!(A,t,vec(τₘ[2:2:end,1]))
    return A
end

function w_taumat(model::TcspcFitModel,tcal::AbstractFloat,
    nbins::Integer)
    # Weights the matrix by the amplitude
    t = -t_vector(tcal,nbins)
    nlcomp = length(model.τ)
    A = zeros(nbins,nlcomp)
    τₘ = τₛ(model)
    α = τₘ[1:2:end,1]
    τ = τₘ[2:2:end,1]
    w_taumat!(A,t,α,τ)
    return A
end

function taumat!(A::AbstractVecOrMat{<:AbstractFloat},
    t::AbstractVector{<:AbstractFloat},
    τ::AbstractVector{<:AbstractFloat})
    # Evaluates in place for everything
    Ai = 1
    for τᵢ ∈ τ
        A[:,Ai] = exp.(t / τᵢ)
        Ai += 1
    end
end

function taumat!(model::TcspcFitModel,
    τₘ::AbstractMatrix{<:AbstractFloat},
    A::AbstractVecOrMat{<:AbstractFloat},
    t::AbstractVector{<:AbstractFloat})
    # Evaluates in place with buffer inputs
    τₛ!(τₘ,model)
    taumat!(A,t,τₘ[2:2:end,1])
end

function w_taumat!(A::AbstractVecOrMat{<:AbstractFloat},
    t::AbstractVector{<:AbstractFloat},
    α::AbstractVector{<:AbstractFloat},
    τ::AbstractVector{<:AbstractFloat})
    # Evaluates weighted tau matrix in place
    taumat!(A,t,τ)
    A *= α'
end

function w_taumat!(model::TcspcFitModel,
    τₘ::AbstractMatrix{<:AbstractFloat},
    A::AbstractVecOrMat{<:AbstractFloat},
    t::AbstractVector{<:AbstractFloat})
    τₛ!(τₘ,model)
    taumat!(A,t,τₘ[1:2:end,1],τₘ[2:2:end,1])
end

function animat(model::TcspcFitModel,tcal::AbstractFloat,
    nbins::Integer)

    t = -t_vector(tcal,nbins)
    # Determine effective lifetimes from anisotropies
    θₘ = θₛ(model)
    τₑₘ = θₘ[2:2:end,[1, 4]]
    τₑ = ((1 ./ τₑₘ[:,1]) + (1 ./ τₑₘ[:,2])) .^ (-1)
    ncomp = length(τₑ)
    A = zeros(nbins,ncomp)
    taumat!(A,t,τₑ)
    return A
end

function animat!(model::TcspcFitModel,
    θₘ::AbstractMatrix{<:AbstractFloat},
    A::AbstractVecOrMat{<:AbstractFloat},
    t::AbstractVector{<:AbstractFloat})
    θₛ!(θₘ,model)
    τₑₘ = θₘ[2:2:end,[1, 4]]
    τₑ = ((1 ./ τₑₘ[:,1]) + (1 ./ τₑₘ[:,2])) .^ (-1)
    taumat!(A,t,τₑ)
end

function w_animat(model::TcspcFitModel,tcal::AbstractFloat,
    nbins::Integer)
    # Weights the matrix by the amplitude
    t = -t_vector(tcal,nbins)
    θₘ = θₛ(model)
    αₑₘ = θₘ[1:2:end,[1, 4]]
    αₑ = αₑₘ[:,1] .* αₑₘ[:,2]
    τₑₘ = θₘ[2:2:end,[1, 4]]
    τₑ = ((1 ./ τₑₘ[:,1]) + (1 ./ τₑₘ[:,2])) .^ (-1)
    ncomp = length(τₑ)
    A = zeros(nbins,ncomp)
    w_taumat!(A,t,αₑ,τₑ)
    return A
end

function w_animat!(model::TcspcFitModel,
    θₘ::AbstractMatrix{<:AbstractFloat},
    A::AbstractVecOrMat{<:AbstractFloat},
    t::AbstractVector{<:AbstractFloat})
    θₛ!(θₘ,model)
    αₑₘ = θₘ[1:2:end,[1, 4]]
    αₑ = αₑₘ[:,1] .* αₑₘ[:,2]
    τₑₘ = θₘ[2:2:end,[1, 4]]
    τₑ = ((1 ./ τₑₘ[:,1]) + (1 ./ τₑₘ[:,2])) .^ (-1)
    w_taumat!(A,t,αₑ,τₑ)
end

function convmat!(A::AbstractVecOrMat,irf::TcspcTransient)
    # Convolves matrix with instrument response function
    for i = 1:size(A,2)
        A[:,i] = abs.((ifft(fftshift(fft(A[:,i]) .* fft(irf.data)))))
    end
end

function model_vec(model::TcspcFitModel,irf::TcspcTransient,
    bck::TcspcTransient)
    # Returns model vector summed up
    A = taumat(model,irf.t.tcal,irf.t.nbins)
    # Convolve with irf
    convmat!(A,irf)
    τₘ = τₛ(model)
    m = A * τₘ[1:2:end,1]
    if nacomp(model) > 0
        m *= (1 / 3)
        Aₐ = animat(model,irf.t.tcal,irf.t.nbins)
        convmat!(Aₐ,irf)
        θₘ = θₛ(model)
        αₑₘ = θₘ[1:2:end,[1, 4]]
        αₑ = αₑₘ[:,1] .* αₑₘ[:,2]
        m += (η(irf.t) / 3) * (Aₐ * αₑ)
    end
    m += model.S.val * irf.data
    m += model.V.val * bck.data
    m .+= model.Z.val
    return m
end

function model_vec_gs(model::TcspcFitModel,irf::TcspcTransient,
    bck::TcspcTransient)
    # Returns model vector summed up
    A = taumat(model,irf.t.tcal,irf.t.nbins)
    fill!(A,0)
    # Convolve with grinvald-steinberg
    τₘ = τₛ(model)
    grin_stein_conv!(A,irf.t.tcal,irf,τₘ[1:2:end,1],τₘ[2:2:end,1])
    m = sum(A,dims=2)
    if nacomp(model) > 0
        m *= (1 / 3)
        Aₐ = animat(model,irf.t.tcal,irf.t.nbins)
        fill!(Aₐ,0)
        θₘ = θₛ(model)
        αₑₘ = θₘ[1:2:end,[1, 4]]
        αₑ = αₑₘ[:,1] .* αₑₘ[:,2]
        τₑₘ = θₘ[2:2:end,[1, 4]]
        τₑ = ((1 ./ τₑₘ[:,1]) + (1 ./ τₑₘ[:,2])) .^ (-1)
        grin_stein_conv!(Aₐ,irf.t.tcal,irf,αₑ,τₑ)
        m += (η(irf.t) / 3) * sum(Aₐ,dims=2)
    end
    m += model.S.val * irf.data
    m += model.V.val * bck.data
    m .+= model.Z.val
    return m
end

function model_vec!(model::TcspcFitModel,
    m::AbstractVector{<:AbstractFloat},
    τₘ::AbstractMatrix{<:AbstractFloat},
    A::AbstractVecOrMat{<:AbstractFloat},
    t::AbstractVector{<:AbstractFloat},
    irf::TcspcTransient,bck::TcspcTransient)
    # Evaluates everything in place with anisotropy
    τₛ!(τₘ,model)
    taumat!(A,t,τₘ[2:2:end,1])
    convmat!(A,irf)
    m[:] = A * τₘ[1:2:end,1]
    m += model.S.val * irf.data
    m += model.V.val * bck.data
    m += model.Z.val
end

function model_vec!(model::TcspcFitModel,
    m::AbstractVector{<:AbstractFloat},
    τₘ::AbstractMatrix{<:AbstractFloat},
    θₘ::AbstractMatrix{<:AbstractFloat},
    A::AbstractVecOrMat{<:AbstractFloat},
    t::AbstractVector{<:AbstractFloat},
    irf::TcspcTransient,bck::TcspcTransient)
    # Evaluates everything in place
    τₛ!(τₘ,model)
    taumat!(A[:,1:nlcomp(model)],t,τₘ[2:2:end,1])
    convmat!(A[:,1:nlcomp(model)],irf)
    m[:] = (1 / 3) * A[:,1:nlcomp(model)] * τₘ[1:2:end,1]
    θₛ!(θₘ,model)
    animat!(model,θₘ,A[:,(nlcomp(model) + 1):end],t)
    convmat!(A[:,(nlcomp(model) + 1):end],irf)
    αₑₘ = θₘ[1:2:end,[1, 4]]
    αₑ = αₑₘ[:,1] .* αₑₘ[:,2]
    m += (η(irf.t) / 3) * (A[:,(nlcomp(model) + 1):end] * αₑ)
    m += model.S.val * irf.data
    m += model.V.val * bck.data
    m += model.Z.val
end


text:
#+FindTag,Internal
tags(o::Extract_Tag_Result)::Array{String,1} = o._tags


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl
=#

insts = OrderedDict()
insts["minimal"] = [
    ((2, false, true, true, true),),
    ((2, false, false, true, true),),
    ((2, true, true, true, true),),
    ((2, true, false, true, true),),
    ((2, false, true, true, false),),
    ((2, false, false, true, false),),
    ((2, true, true, true, false),),
    ((2, true, false, true, false),),
    ((2, false, false, false, false),),
    ((2, true, false, false, false),),
    ]
insts["fast"] = [
    ((5, false, true, true, true),),
    ((5, false, false, true, false),),
    ((5, true, true, true, false),),
    ((5, true, false, true, true),),
    ((5, false, true, false, false),),
    ((5, false, false, false, false),),
    ((5, true, true, false, false),),
    ((5, true, false, false, false),),
    ((20, false, true, true, false),),
    ((20, false, false, true, true),),
    ((20, true, true, true, true),),
    ((20, true, false, true, false),),
    ((20, false, true, false, false),),
    ((20, false, false, false, false),),
    ((20, true, true, false, false),),
    ((20, true, false, false, false),),
    ((100, false, true, false, false),),
    ((100, false, false, false, false),),
    ]
insts["various"] = [
    ((50, false, true, true, true),),
    ((50, false, false, true, true),),
    ((50, true, true, true, false),),
    ((50, true, false, true, false),),
    ((50, false, true, false, false),),
    ((50, false, false, false, false),),
    ((50, true, true, false, false),),
    ((50, true, false, false, false),),
    ((100, false, true, true, true),),
    ((100, false, false, true, false),),
    ((100, true, true, true, false),),
    ((100, true, false, true, true),),
    ((100, false, true, false, false),),
    ((100, false, false, false, false),),
    ((100, true, true, false, false),),
    ((100, true, false, false, false),),
    ((200, false, true, false, false),),
    ((200, false, false, false, false),),
    ((200, false, true, true, false),),
    ((200, false, false, true, false),),
    ((200, true, true, true, true),),
    ((200, true, false, true, true),),
    ((400, false, true, false, false),),
    ((400, false, false, false, false),),
    ((400, false, true, true, false),),
    ((400, false, false, true, true),),
    ((400, true, true, true, true),),
    ((400, true, false, true, false),),
    ]
return (NearestPSDJuMP, insts)


text:
"""
Convert a type into its corresponding Vulkan type.

### Examples
```jldoctest
julia> to_vk(UInt32, v"1")
0x00400000

julia> to_vk(NTuple{6, UInt8}, "hello")
(0x68, 0x65, 0x6c, 0x6c, 0x6f, 0x00)
```
"""
function to_vk end

to_vk(T, x) = convert(T, x)
to_vk(T::Type{<:NTuple}, x) = to_vk.(eltype(T), x)
to_vk(T::Type{UInt32}, version::VersionNumber) = VK_MAKE_VERSION(version.major, version.minor, version.patch)
to_vk(T::Type{NTuple{N,UInt8}}, s::AbstractString) where {N} = T(s * '\0' ^ (N - length(s)))

"""
Convert a Vulkan type into its corresponding Julia type.

### Examples
```jldoctest
julia> from_vk(VersionNumber, UInt32(VK_MAKE_VERSION(1, 2, 3)))
v"1.2.3"

julia> from_vk(String, (0x68, 0x65, 0x6c, 0x6c, 0x6f, 0x00))
"hello"

julia> from_vk(Bool, UInt32(1))
true
```
"""
function from_vk end

from_vk(T::Type{<:VulkanStruct{false}}, x) = T(x)
from_vk(T::Type{<:VulkanStruct{true}}, x) = T(x, [])
from_vk(T, x) = convert(T, x)
from_vk(T::Type{<:NTuple}, x) = from_vk.(eltype(T), x)
from_vk(T::Type{VersionNumber}, version::UInt32) = T(VK_VERSION_MAJOR(version), VK_VERSION_MINOR(version), VK_VERSION_PATCH(version))

function from_vk(T::Type{S}, str::NTuple{N,UInt8}) where {N,S <: AbstractString}
    nullchar = findfirst(iszero, str)
    if isnothing(nullchar)
        T(collect(str))
    else
        T(collect(str[1:nullchar - 1]))
    end
end


text:
function typeinf(interp::CustomInterpreter, frame::InferenceState)
    interp.depth[] += 1
    ret = @invoke typeinf(interp::AbstractInterpreter, frame::InferenceState)
    interp.depth[] -= 1

    return ret
end

function _typeinf(interp::CustomInterpreter, frame::InferenceState)
    ret = @invoke _typeinf(interp::AbstractInterpreter, frame::InferenceState)

    return ret
end

function transform_result_for_cache(interp::CustomInterpreter, linfo::MethodInstance,
                                    @nospecialize(inferred_result))
    ret = @invoke transform_result_for_cache(interp::AbstractInterpreter, linfo::MethodInstance,
                                             inferred_result)

    return ret
end

function cache_result!(interp::CustomInterpreter, result::InferenceResult, valid_worlds::WorldRange)
    ret = @invoke cache_result!(interp::AbstractInterpreter, result::InferenceResult, valid_worlds::WorldRange)

    return ret
end

function finish(me::InferenceState, interp::CustomInterpreter)
    ret = @invoke finish(me::InferenceState, interp::AbstractInterpreter)

    return ret
end

function finish(src::CodeInfo, interp::CustomInterpreter)
    ret = @invoke finish(me::CodeInfo, interp::AbstractInterpreter)

    return ret
end

function is_same_frame(interp::CustomInterpreter, linfo::MethodInstance, frame::InferenceState)
    ret = @invoke is_same_frame(interp::AbstractInterpreter, linfo::MethodInstance, frame::InferenceState)

    return ret
end

function resolve_call_cycle!(interp::CustomInterpreter, linfo::MethodInstance, parent::InferenceState)
    ret = @invoke resolve_call_cycle!(interp::AbstractInterpreter, linfo::MethodInstance, parent::InferenceState)

    return ret
end

function typeinf_edge(interp::CustomInterpreter, method::Method, @nospecialize(atypes), sparams::SimpleVector, caller::InferenceState)
    ret = @invoke typeinf_edge(interp::AbstractInterpreter, method::Method, atypes, sparams::SimpleVector, caller::InferenceState)

    return ret
end

max_stars_repo_name: aviatesk/AbstractInterpreterPlayground

max_issues_repo_name: aviatesk/AbstractInterpreterPlayground

max_forks_repo_name: aviatesk/AbstractInterpreterPlayground


text:
module TestMOIwrapper

using CPLEX
using MathOptInterface
using Test

const MOI  = MathOptInterface
const MOIT = MOI.Test
const MOIB = MOI.Bridges

const CONFIG = MOIT.TestConfig(basis = true)

const OPTIMIZER = CPLEX.Optimizer()
MOI.set(OPTIMIZER, MOI.Silent(), true)
# Turn off presolve reductions so CPLEX will generate infeasibility
# certificates.
MOI.set(OPTIMIZER, MOI.RawParameter("CPX_PARAM_REDUCE"), 0)

const BRIDGED_OPTIMIZER = MOI.Bridges.full_bridge_optimizer(OPTIMIZER, Float64)

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

function test_unittest()
    MOIT.unittest(BRIDGED_OPTIMIZER, CONFIG, [
        # TODO(odow): bug! We can't delete a vector of variables  if one is in
        # a second order cone.
        "delete_soc_variables",
    ])
end

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

function test_contlineartest()
    MOIT.contlineartest(BRIDGED_OPTIMIZER, CONFIG)
end

function test_intlineartest()
    # interval somehow needed for indicator tests
    interval_optimizer = MOIB.LazyBridgeOptimizer(OPTIMIZER)
    MOIB.add_bridge(interval_optimizer, MOIB.Constraint.SplitIntervalBridge{Float64})
    MOIT.intlineartest(BRIDGED_OPTIMIZER, CONFIG)
    MOIT.intlineartest(interval_optimizer, CONFIG)
end

function test_contquadratictest()
    MOIT.contquadratictest(
        BRIDGED_OPTIMIZER,
        MOIT.TestConfig(atol = 1e-3, rtol = 1e-3),
        ["ncqcp"], # CPLEX doesn't support non-convex problems
    )
end

function test_CPXPARAM_OptimalityTarget()
    # Test setting CPXPARAM_OptimalityTarget because it changes the problem
    # type.
    # Max x^2
    # s.t. 1 <= x <= 4
    model = CPLEX.Optimizer()
    MOI.set(model, MOI.Silent(), true)
    x = MOI.add_variable(model)
    MOI.add_constraint(model, MOI.SingleVariable(x), MOI.Interval(1.0, 4.0))
    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}(),
        MOI.ScalarQuadraticFunction(
            MOI.ScalarAffineTerm{Float64}[],
            [MOI.ScalarQuadraticTerm(2.0, x, x)],
            0.0,
        )
    )
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    MOI.set(
        model,
        MOI.RawParameter("CPXPARAM_OptimalityTarget"),
        CPX_OPTIMALITYTARGET_OPTIMALGLOBAL,
    )
    MOI.optimize!(model)
    @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL
    @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
    @test MOI.get(model, MOI.ObjectiveValue()) ≈ 16.0 atol=1e-6
    @test MOI.get(model, MOI.VariablePrimal(), x) ≈ 4.0 atol=1e-6
end

function test_contconic()
    MOIT.lintest(BRIDGED_OPTIMIZER, CONFIG)

    soc_config = MOIT.TestConfig(atol=5e-3)

    # TODO(odow): investigate why infeasibility certificates not generated for
    # SOC.
    MOIT.soctest(BRIDGED_OPTIMIZER, soc_config, ["soc3"])
    MOIT.soc3test(
        BRIDGED_OPTIMIZER,
        MOIT.TestConfig(atol = 1e-3, infeas_certificates = false)
    )

    MOIT.rsoctest(BRIDGED_OPTIMIZER, soc_config, ["rotatedsoc2"])
    MOIT.rotatedsoc2test(
        BRIDGED_OPTIMIZER,
        # Need for `duals = false` fixed by https://github.com/jump-dev/MathOptInterface.jl/pull/1171
        # Remove in a future MOI 0.9.18+ release.
        MOIT.TestConfig(
            atol = 1e-3, infeas_certificates = false, duals = false
        ),
    )

    MOIT.geomeantest(BRIDGED_OPTIMIZER, soc_config)
end

function test_solvername()
    @test MOI.get(BRIDGED_OPTIMIZER, MOI.SolverName()) == "CPLEX"
end

function test_default_objective_test()
    MOIT.default_objective_test(BRIDGED_OPTIMIZER)
end

function test_default_status_test()
    MOIT.default_status_test(BRIDGED_OPTIMIZER)
end

function test_nametest()
    MOIT.nametest(BRIDGED_OPTIMIZER)
end

function test_validtest()
    MOIT.validtest(BRIDGED_OPTIMIZER)
end

function test_emptytest()
    MOIT.emptytest(BRIDGED_OPTIMIZER)
end

function test_orderedindicestest()
    MOIT.orderedindicestest(BRIDGED_OPTIMIZER)
end

function test_copytest()
    MOIT.copytest(
        BRIDGED_OPTIMIZER,
        MOI.Bridges.full_bridge_optimizer(CPLEX.Optimizer(), Float64)
    )
end

function test_scalar_function_constant_not_zero()
    MOIT.scalar_function_constant_not_zero(OPTIMIZER)
end

function test_start_values_test()
    model = CPLEX.Optimizer()
    x = MOI.add_variables(model, 2)
    @test MOI.supports(model, MOI.VariablePrimalStart(), MOI.VariableIndex)
    @test MOI.get(model, MOI.VariablePrimalStart(), x[1]) === nothing
    @test MOI.get(model, MOI.VariablePrimalStart(), x[2]) === nothing
    MOI.set(model, MOI.VariablePrimalStart(), x[1], 1.0)
    MOI.set(model, MOI.VariablePrimalStart(), x[2], nothing)
    @test MOI.get(model, MOI.VariablePrimalStart(), x[1]) == 1.0
    @test MOI.get(model, MOI.VariablePrimalStart(), x[2]) === nothing
    MOI.optimize!(model)
    @test MOI.get(model, MOI.ObjectiveValue()) == 0.0
end

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

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

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

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

function test_automatic_env()
    model_1 = CPLEX.Optimizer()
    model_2 = CPLEX.Optimizer()
    @test model_1.env.ptr !== model_2.env.ptr
end

function test_user_provided_env_empty()
    env = CPLEX.Env()
    model = CPLEX.Optimizer(env)
    @test model.env === env
    @test env.ptr != C_NULL
    MOI.empty!(model)
    @test model.env === env
    @test env.ptr != C_NULL
end

function test_automatic_env_empty()
    model = CPLEX.Optimizer()
    env = model.env
    MOI.empty!(model)
    @test model.env === env
    @test env.ptr != C_NULL
end

function test_manual_env()
    env = CPLEX.Env()
    model = CPLEX.Optimizer(env)
    finalize(env)
    @test env.finalize_called
    finalize(model)
    @test env.ptr == C_NULL
end

function test_cont_int_cont()
    atol = 1e-5
    rtol = 1e-5

    model = CPLEX.Optimizer()
    MOI.empty!(model)
    @test MOI.is_empty(model)

    # min -x
    # st   x + y <= 1.5   (x + y - 1.5 ∈ Nonpositives)
    #       x, y >= 0   (x, y ∈ Nonnegatives)

    v = MOI.add_variables(model, 2)
    @test MOI.get(model, MOI.NumberOfVariables()) == 2

    cf = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0,1.0], v), 0.0)
    c = MOI.add_constraint(model, cf, MOI.LessThan(1.5))
    @test MOI.get(model, MOI.NumberOfConstraints{MOI.ScalarAffineFunction{Float64},MOI.LessThan{Float64}}()) == 1

    MOI.add_constraint.(model, MOI.SingleVariable.(v), MOI.GreaterThan(0.0))
    @test MOI.get(model, MOI.NumberOfConstraints{MOI.SingleVariable,MOI.GreaterThan{Float64}}()) == 2

    objf = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([-1.0,0.0], v), 0.0)
    MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), objf)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    @test MOI.get(model, MOI.ObjectiveSense()) == MOI.MIN_SENSE
    @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED

    MOI.optimize!(model)

    @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL
    @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
    @test MOI.get(model, MOI.ObjectiveValue()) ≈ -1.5 atol=atol rtol=rtol
    @test MOI.get(model, MOI.VariablePrimal(), v) ≈ [1.5, 0] atol=atol rtol=rtol
    @test MOI.get(model, MOI.ConstraintPrimal(), c) ≈ 1.5 atol=atol rtol=rtol
    @test MOI.get(model, MOI.DualStatus()) == MOI.FEASIBLE_POINT
    @test MOI.get(model, MOI.ConstraintDual(), c) ≈ -1.0 atol=atol rtol=rtol

    # Add integrality constraints
    int = MOI.add_constraint.(model, MOI.SingleVariable.(v), MOI.Integer())
    MOI.optimize!(model)
    @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL
    @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
    @test MOI.get(model, MOI.ObjectiveValue()) ≈ -1.0 atol=atol rtol=rtol
    @test MOI.get(model, MOI.VariablePrimal(), v) ≈ [1.0, 0] atol=atol rtol=rtol
    @test MOI.get(model, MOI.ConstraintPrimal(), c) ≈ 1.0 atol=atol rtol=rtol
    @test MOI.get(model, MOI.DualStatus()) == MOI.NO_SOLUTION

    # Remove integrality constraints
    MOI.delete.(model, int)
    MOI.optimize!(model)
    @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL
    @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
    @test MOI.get(model, MOI.ObjectiveValue()) ≈ -1.5 atol=atol rtol=rtol
    @test MOI.get(model, MOI.VariablePrimal(), v) ≈ [1.5, 0] atol=atol rtol=rtol
    @test MOI.get(model, MOI.ConstraintPrimal(), c) ≈ 1.5 atol=atol rtol=rtol
    @test MOI.get(model, MOI.DualStatus()) == MOI.FEASIBLE_POINT
    @test MOI.get(model, MOI.ConstraintDual(), c) ≈ -1.0 atol=atol rtol=rtol
end

function test_conflict_bounds()
    # @testset "Variable bounds (SingleVariable and LessThan/GreaterThan)" begin
    # Test similar to ../C_API/iis.jl, but ported to MOI.
    model = CPLEX.Optimizer()
    x = MOI.add_variable(model)
    c1 = MOI.add_constraint(model, MOI.SingleVariable(x), MOI.GreaterThan(2.0))
    c2 = MOI.add_constraint(model, MOI.SingleVariable(x), MOI.LessThan(1.0))

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

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

function test_conflict_scalaraffine()
    # @testset "Variable bounds (ScalarAffine)" begin
    # Same test as ../C_API/iis.jl, but ported to MOI.
    model = CPLEX.Optimizer()
    x = MOI.add_variable(model)
    c1 = MOI.add_constraint(model, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0], [x]), 0.0), MOI.GreaterThan(2.0))
    c2 = MOI.add_constraint(model, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0], [x]), 0.0), MOI.LessThan(1.0))

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

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

function test_conflict_two_bound()
    model = CPLEX.Optimizer()
    x = MOI.add_variable(model)
    y = MOI.add_variable(model)
    b1 = MOI.add_constraint(model, MOI.SingleVariable(x), MOI.GreaterThan(0.0))
    b2 = MOI.add_constraint(model, MOI.SingleVariable(y), MOI.GreaterThan(0.0))
    cf1 = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 1.0], [x, y]), 0.0)
    c1 = MOI.add_constraint(model, cf1, MOI.LessThan(-1.0))
    cf2 = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, -1.0], [x, y]), 0.0)
    c2 = MOI.add_constraint(model, cf2, MOI.GreaterThan(1.0))

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

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

function test_conflict_two_equalto()
    model = CPLEX.Optimizer()
    x = MOI.add_variable(model)
    y = MOI.add_variable(model)
    b1 = MOI.add_constraint(model, MOI.SingleVariable(x), MOI.GreaterThan(0.0))
    b2 = MOI.add_constraint(model, MOI.SingleVariable(y), MOI.GreaterThan(0.0))
    cf1 = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 1.0], [x, y]), 0.0)
    c1 = MOI.add_constraint(model, cf1, MOI.EqualTo(-1.0))
    cf2 = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, -1.0], [x, y]), 0.0)
    c2 = MOI.add_constraint(model, cf2, MOI.GreaterThan(1.0))

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

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

function test_conflict_variables_outside()
    model = CPLEX.Optimizer()
    x = MOI.add_variable(model)
    y = MOI.add_variable(model)
    z = MOI.add_variable(model)
    b1 = MOI.add_constraint(model, MOI.SingleVariable(x), MOI.GreaterThan(0.0))
    b2 = MOI.add_constraint(model, MOI.SingleVariable(y), MOI.GreaterThan(0.0))
    b3 = MOI.add_constraint(model, MOI.SingleVariable(z), MOI.GreaterThan(0.0))
    cf1 = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 1.0], [x, y]), 0.0)
    c1 = MOI.add_constraint(model, cf1, MOI.LessThan(-1.0))
    cf2 = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, -1.0, 1.0], [x, y, z]), 0.0)
    c2 = MOI.add_constraint(model, cf2, MOI.GreaterThan(1.0))

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

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

function test_conflict_no_conflict()
    model = CPLEX.Optimizer()
    x = MOI.add_variable(model)
    c1 = MOI.add_constraint(model, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0], [x]), 0.0), MOI.GreaterThan(1.0))
    c2 = MOI.add_constraint(model, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0], [x]), 0.0), MOI.LessThan(2.0))

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

    # Once it's called, no problem.
    MOI.compute_conflict!(model)
    @test MOI.get(model, CPLEX.ConflictStatus()) == CPLEX.CPX_STAT_CONFLICT_FEASIBLE
    @test MOI.get(model, MOI.ConflictStatus()) == MOI.NO_CONFLICT_EXISTS
    @test MOI.get(model, MOI.ConstraintConflictStatus(), c1) == MOI.NOT_IN_CONFLICT
    @test MOI.get(model, MOI.ConstraintConflictStatus(), c2) == MOI.NOT_IN_CONFLICT
end

function test_ZeroOne_NONE()
    model = CPLEX.Optimizer()
    x = MOI.add_variable(model)
    c = MOI.add_constraint(model, MOI.SingleVariable(x), MOI.ZeroOne())
    tmp = Ref{Cdouble}()
    CPLEX.CPXgetlb(model.env, model.lp, tmp, 0, 0)
    @test tmp[] == 0.0
    CPLEX.CPXgetub(model.env, model.lp, tmp, 0, 0)
    @test tmp[] == 1.0
    MOI.delete(model, c)
    CPLEX.CPXgetlb(model.env, model.lp, tmp, 0, 0)
    @test tmp[] == -CPLEX.CPX_INFBOUND
    CPLEX.CPXgetub(model.env, model.lp, tmp, 0, 0)
    @test tmp[] == CPLEX.CPX_INFBOUND
end

function test_ZeroOne_LESS_THAN()
    model = CPLEX.Optimizer()
    x = MOI.add_variable(model)
    MOI.add_constraint(model, MOI.SingleVariable(x), MOI.LessThan(2.0))
    c = MOI.add_constraint(model, MOI.SingleVariable(x), MOI.ZeroOne())
    tmp = Ref{Cdouble}()
    CPLEX.CPXgetlb(model.env, model.lp, tmp, 0, 0)
    @test tmp[] == 0.0
    CPLEX.CPXgetub(model.env, model.lp, tmp, 0, 0)
    @test tmp[] == 2.0
    MOI.delete(model, c)
    CPLEX.CPXgetlb(model.env, model.lp, tmp, 0, 0)
    @test tmp[] == -CPLEX.CPX_INFBOUND
    CPLEX.CPXgetub(model.env, model.lp, tmp, 0, 0)
    @test tmp[] == 2.0
end

function test_ZeroOne_GREATER_THAN()
    model = CPLEX.Optimizer()
    x = MOI.add_variable(model)
    MOI.add_constraint(model, MOI.SingleVariable(x), MOI.GreaterThan(-2.0))
    c = MOI.add_constraint(model, MOI.SingleVariable(x), MOI.ZeroOne())
    tmp = Ref{Cdouble}()
    CPLEX.CPXgetlb(model.env, model.lp, tmp, 0, 0)
    @test tmp[] == -2.0
    CPLEX.CPXgetub(model.env, model.lp, tmp, 0, 0)
    @test tmp[] == 1.0
    MOI.delete(model, c)
    CPLEX.CPXgetlb(model.env, model.lp, tmp, 0, 0)
    @test tmp[] == -2.0
    CPLEX.CPXgetub(model.env, model.lp, tmp, 0, 0)
    @test tmp[] == CPLEX.CPX_INFBOUND
end

function test_ZeroOne_INTERVAL()
    model = CPLEX.Optimizer()
    x = MOI.add_variable(model)
    MOI.add_constraint(model, MOI.SingleVariable(x), MOI.Interval(-2.0, 2.0))
    c = MOI.add_constraint(model, MOI.SingleVariable(x), MOI.ZeroOne())
    tmp = Ref{Cdouble}()
    CPLEX.CPXgetlb(model.env, model.lp, tmp, 0, 0)
    @test tmp[] == -2.0
    CPLEX.CPXgetub(model.env, model.lp, tmp, 0, 0)
    @test tmp[] == 2.0
    MOI.delete(model, c)
    CPLEX.CPXgetlb(model.env, model.lp, tmp, 0, 0)
    @test tmp[] == -2.0
    CPLEX.CPXgetub(model.env, model.lp, tmp, 0, 0)
    @test tmp[] == 2.0
end

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

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

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

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

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

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

end  # module TestMOIwrapper

runtests(TestMOIwrapper)


text:
@with_kw struct Shader
    code::String
    stage::ShaderStageFlag
    bindings::Vector{DescriptorBinding} = DescriptorBinding[]
    entry_point::Symbol = :main
end

Shader(shader::ShaderFile, bindings; entry_point=:main) = Shader(code=spirv_code(shader), stage=shader.stage; bindings, entry_point)

has_bindings(shader::Shader) = !isempty(shader.bindings)

PipelineShaderStageCreateInfo(shader::Shader) = PipelineShaderStageCreateInfo(shader.stage, shader.mod, string(shader.entry_point))

function create_descriptor_set_layouts(shaders::AbstractVector{Shader})::Vector{DescriptorSetLayout}
    sets = DefaultOrderedDict{Int,Vector{_DescriptorSetLayoutBinding}}(() -> _DescriptorSetLayoutBinding[])
    for shader ∈ shaders
        for resource_binding ∈ shader.bindings
            push!(sets[resource_binding.set], _DescriptorSetLayoutBinding(resource_binding.binding, resource_binding.descriptor_type, shader.stage; descriptor_count=1))
        end
    end
    if !all(keys(sets) .== 0:(length(sets) - 1))
        error("Invalid layout description (non-contiguous sets from 0) in $sets.")
    end
    layout_bindings_vec = collect(values(sets))
    device = first(shaders).mod.device
    [DescriptorSetLayout(device, bindings) for bindings ∈ layout_bindings_vec]
end


text:
module KeyboardDriven

import Base.|

using Match
using TERMIOS

export
    CharKey,
    ModifiedKey,
    SpecialKey,
    alt,
    ctrl,
    getkeypress,
    keyboarddriven,
    prompt,
    shift

include("modes.jl")
include("keys.jl")
include("input.jl")

end # module


text:
using Conda

Conda.add("ase", channel = "conda-forge")
Conda.add("rdkit", channel = "conda-forge")
Conda.add("pymatgen", channel = "conda-forge")
Conda.rm("mkl")
Conda.add("mkl")


text:
#### Construct Normal spline (ReproducingKernel_0)

function _prepare(nodes::AbstractVecOfSVecs{n,T}, kernel::ReproducingKernel_0) where {n,T}
    # Normalize nodes out-of-place to avoid aliasing
    min_bound, max_bound, scale = _normalization_scaling(nodes)
    nodes = _normalize.(nodes, (min_bound,), (max_bound,), scale)

    if kernel.ε == 0
        kernel = _estimate_ε(kernel, nodes)
    end

    gram     = _gram(nodes, kernel)
    chol     = cholesky(gram)
    cond     = _estimate_cond(gram, chol)

    return NormalSpline{n,T,typeof(kernel)}(; _kernel = kernel, _nodes = nodes, _gram = gram, _chol = chol, _cond = cond, _min_bound = min_bound, _max_bound = max_bound, _scale = scale)
end

function _construct!(
        spl::NormalSpline{n,T,RK},
        values::AbstractVector{T},
    ) where {n, T, RK <: ReproducingKernel_0}
    n₁ = length(values)
    length(_get_nodes(spl)) != n₁ && error("Number of data values ($n₁) does not correspond to the number of nodes $(length(_get_nodes(spl))).")
    size(_get_chol(spl)) != (n₁, n₁) && error("Number of data values ($n₁) does not correspond to the size of the Gram matrix ($(size(_get_chol(spl)))).")

    # Resize buffers
    resize!(_get_values(spl), n₁)
    empty!(_get_d_nodes(spl))
    empty!(_get_d_dirs(spl))
    empty!(_get_d_values(spl))
    resize!(_get_mu(spl), n₁)
    resize!(_get_rhs(spl), n₁)

    # Copy values to avoid aliasing
    _get_values(spl) .= _get_rhs(spl) .= values

    # Compute spline coefficients
    ldiv!(_get_mu(spl), _get_chol(spl), _get_rhs(spl))

    return spl
end

#### Construct Normal spline (ReproducingKernel_1)

function _prepare(nodes::AbstractVecOfSVecs{n,T}, d_nodes::AbstractVecOfSVecs{n,T}, d_dirs::AbstractVecOfSVecs{n,T}, kernel::ReproducingKernel_1) where {n,T}
    # Normalize inputs out-of-place to avoid aliasing
    min_bound, max_bound, scale = _normalization_scaling(nodes, d_nodes)
    nodes   = _normalize.(nodes, (min_bound,), (max_bound,), scale)
    d_nodes = _normalize.(d_nodes, (min_bound,), (max_bound,), scale)
    d_dirs  = d_dirs ./ norm.(d_dirs)

    if kernel.ε == 0
        kernel = _estimate_ε(kernel, nodes, d_nodes)
    end

    gram     = _gram(nodes, d_nodes, d_dirs, kernel)
    chol     = cholesky(gram)
    cond     = _estimate_cond(gram, chol)

    NormalSpline{n,T,typeof(kernel)}(; _kernel = kernel, _nodes = nodes, _d_nodes = d_nodes, _d_dirs = d_dirs, _gram = gram, _chol = chol, _cond = cond, _min_bound = min_bound, _max_bound = max_bound, _scale = scale)
end

function _construct!(
        spl::NormalSpline{n,T,RK},
        values::AbstractVector{T},
        d_values::AbstractVector{T},
    ) where {n, T, RK <: ReproducingKernel_1}
    n₁ = length(values)
    n₂ = length(d_values)
    length(_get_nodes(spl)) != n₁ && error("Number of data values ($n₁) does not correspond to the number of nodes $(length(_get_nodes(spl))).")
    length(_get_d_nodes(spl)) != n₂ && error("Number of derivative values ($n₂) does not correspond to the number of derivative nodes.")
    size(_get_chol(spl)) != (n₁+n₂, n₁+n₂) && error("Number of data and derivative values ($(n₁+n₂)) do not correspond to the size of the Gram matrix ($(size(_get_chol(spl)))).")

    # Resize buffers
    resize!(_get_values(spl), n₁)
    resize!(_get_d_nodes(spl), n₂)
    resize!(_get_d_dirs(spl), n₂)
    resize!(_get_d_values(spl), n₂)
    resize!(_get_mu(spl), n₁+n₂)
    resize!(_get_rhs(spl), n₁+n₂)

    # Copy values to avoid aliasing
    _get_values(spl) .= view(_get_rhs(spl), 1:n₁) .= values

    # Nodes scaled down by `scale` -> directional derivative scaled up by `scale`; allocate new array to avoid aliasing
    _get_d_values(spl) .= view(_get_rhs(spl), n₁+1:n₁+n₂) .= _get_scale(spl) .* d_values

    # Compute spline coefficients and construct spline
    ldiv!(_get_mu(spl), _get_chol(spl), _get_rhs(spl))

    return spl
end

#### Evaluation

@inline function _evaluate!(
        vals::AbstractArray{<:Any,N},
        spl::AbstractNormalSpline{n,<:Any,<:ReproducingKernel_0},
        points::AbstractArrOfSVecs{n,<:Any,N},
    ) where {n, N}
    @inbounds for i in 1:length(points)
        vals[i] = _evaluate(spl, points[i])
    end
    return vals
end

@inline function _evaluate(
        spl::AbstractNormalSpline{n,<:Any,RK},
        x::SVector{n},
    ) where {n, RK <: ReproducingKernel_0}
    kernel, nodes, d_nodes, d_dirs, mu =
        _get_kernel(spl), _get_nodes(spl), _get_d_nodes(spl), _get_d_dirs(spl), _get_mu(spl)
    n₁ = length(nodes)
    n₂ = length(d_nodes)
    x  = _normalize(spl, x)
    v  = zero(promote_type(eltype(spl), eltype(kernel), eltype(x)))
    @inbounds for i in 1:n₁
        v += mu[i] * _rk(kernel, x, nodes[i])
    end
    if RK <: ReproducingKernel_1
        @inbounds for i in 1:n₂
            v += mu[i+n₁] * _∂rk_∂e(kernel, x, d_nodes[i], d_dirs[i])
        end
    end
    return v
end

@inline function _evaluate_gradient(
        spl::AbstractNormalSpline{n,<:Any,RK},
        x::SVector{n},
    ) where {n, RK <: ReproducingKernel_0}
    kernel, nodes, d_nodes, d_dirs, mu, scale =
        _get_kernel(spl), _get_nodes(spl), _get_d_nodes(spl), _get_d_dirs(spl), _get_mu(spl), _get_scale(spl)
    n₁ = length(nodes)
    n₂ = length(d_nodes)
    x  = _normalize(spl, x)
    ∇  = zero(SVector{n,promote_type(eltype(spl), eltype(kernel), eltype(x))})
    @inbounds for i in 1:n₁
        ∇ += mu[i] * _∂rk_∂η(kernel, x, nodes[i])
    end
    if RK <: ReproducingKernel_1
        @inbounds for i in 1:n₂
            ∇² = _∂²rk_∂η∂ξ(kernel, x, d_nodes[i])
            ∇ += mu[i+n₁] * (∇² * d_dirs[i])
        end
    end
    return ∇ ./ scale
end


text:
export DeepCFR

"""
    DeepCFR(;kwargs...)

Symbols used here follow the paper: [Deep Counterfactual Regret Minimization](https://arxiv.org/abs/1811.00164)

# Keyword arguments

- `K`, number of traverrsal.
- `t`, number of iteration.
- `Π`, the policy network.
- `V`, a dictionary of each player's advantage network.
- `MΠ`, a strategy memory.
- `MV`, a dictionary of each player's advantage memory.
- `reinitialize_freq=1`, the frequency of reinitializing the value networks.
"""
Base.@kwdef mutable struct DeepCFR{TP,TV,TMP,TMV,I,R,P} <: AbstractCFRPolicy
    Π::TP
    V::TV
    MΠ::TMP
    MV::TMV
    K::Int = 20
    t::Int = 1
    reinitialize_freq::Int = 1
    batch_size_V::Int = 32
    batch_size_Π::Int = 32
    n_training_steps_V::Int = 1
    n_training_steps_Π::Int = 1
    rng::R = Random.GLOBAL_RNG
    initializer::I = glorot_normal(rng)
    max_grad_norm::Float32 = 10.0f0
    # for logging
    Π_losses::Vector{Float32} = zeros(Float32, n_training_steps_Π)
    V_losses::Dict{P,Vector{Float32}} =
        Dict(k => zeros(Float32, n_training_steps_V) for (k, _) in MV)
    Π_norms::Vector{Float32} = zeros(Float32, n_training_steps_Π)
    V_norms::Dict{P,Vector{Float32}} =
        Dict(k => zeros(Float32, n_training_steps_V) for (k, _) in MV)
end

function RLBase.prob(π::DeepCFR, env::AbstractEnv)
    I = send_to_device(device(π.Π), state(env))
    m = send_to_device(device(π.Π), ifelse.(legal_action_space_mask(env), 0.0f0, -Inf32))
    logits = π.Π(Flux.unsqueeze(I, ndims(I) + 1)) |> vec
    σ = softmax(logits .+ m)
    send_to_host(σ)
end

(π::DeepCFR)(env::AbstractEnv) =
    sample(π.rng, action_space(env), Weights(prob(π, env), 1.0))

"Run one interation"
function RLBase.update!(π::DeepCFR, env::AbstractEnv)
    for p in players(env)
        if p != chance_player(env)
            for k in 1:π.K
                external_sampling!(π, copy(env), p)
            end
            update_advantage_networks(π, p)
        end
    end
    π.t += 1
end

"Update Π (policy network)"
function RLBase.update!(π::DeepCFR)
    Π = π.Π
    Π_losses = π.Π_losses
    Π_norms = π.Π_norms
    D = device(Π)
    MΠ = π.MΠ
    ps = Flux.params(Π)

    for x in ps
        x .= π.initializer(size(x)...)
    end

    for i in 1:π.n_training_steps_Π
        batch_inds = rand(π.rng, 1:length(MΠ), π.batch_size_Π)
        I = send_to_device(D, Flux.batch([MΠ[:I][i] for i in batch_inds]))
        σ = send_to_device(D, Flux.batch([MΠ[:σ][i] for i in batch_inds]))
        t = send_to_device(D, Flux.batch([MΠ[:t][i] / π.t for i in batch_inds]))
        m = send_to_device(
            D,
            Flux.batch([ifelse.(MΠ[:m][i], 0.0f0, -Inf32) for i in batch_inds]),
        )
        gs = gradient(ps) do
            logits = Π(I) .+ m
            loss = mean(reshape(t, 1, :) .* ((σ .- softmax(logits)) .^ 2))
            ignore() do
                # println(σ, "!!!",m, "===", Π(I))
                Π_losses[i] = loss
            end
            loss
        end
        Π_norms[i] = clip_by_global_norm!(gs, ps, π.max_grad_norm)
        update!(Π, gs)
    end
end

"Update advantage network"
function update_advantage_networks(π, p)
    V = π.V[p]
    V_losses = π.V_losses[p]
    V_norms = π.V_norms[p]
    MV = π.MV[p]
    if π.t % π.reinitialize_freq == 0
        for x in Flux.params(V)
            # TODO: inplace
            x .= π.initializer(size(x)...)
        end
    end
    if length(MV) >= π.batch_size_V
        for i in 1:π.n_training_steps_V
            batch_inds = rand(π.rng, 1:length(MV), π.batch_size_V)
            I = send_to_device(device(V), Flux.batch([MV[:I][i] for i in batch_inds]))
            r̃ = send_to_device(device(V), Flux.batch([MV[:r̃][i] for i in batch_inds]))
            t = send_to_device(device(V), Flux.batch([MV[:t][i] / π.t for i in batch_inds]))
            m = send_to_device(device(V), Flux.batch([MV[:m][i] for i in batch_inds]))
            ps = Flux.params(V)
            gs = gradient(ps) do
                loss = mean(reshape(t, 1, :) .* ((r̃ .- V(I) .* m) .^ 2))
                ignore() do
                    V_losses[i] = loss
                end
                loss
            end
            V_norms[i] = clip_by_global_norm!(gs, ps, π.max_grad_norm)
            update!(V, gs)
        end
    end
end

"CFR Traversal with External Sampling"
function external_sampling!(π::DeepCFR, env::AbstractEnv, p)
    if is_terminated(env)
        reward(env, p)
    elseif current_player(env) == chance_player(env)
        env(rand(π.rng, action_space(env)))
        external_sampling!(π, env, p)
    elseif current_player(env) == p
        V = π.V[p]
        s = state(env)
        I = send_to_device(device(V), Flux.unsqueeze(s, ndims(s) + 1))
        A = action_space(env)
        m = legal_action_space_mask(env)
        σ = masked_regret_matching(V(I) |> send_to_host |> vec, m)
        v = zeros(length(σ))
        v̄ = 0.0
        for i in 1:length(m)
            if m[i]
                v[i] = external_sampling!(π, child(env, A[i]), p)
                v̄ += σ[i] * v[i]
            end
        end
        push!(π.MV[p], I = s, t = π.t, r̃ = (v .- v̄) .* m, m = m)
        v̄
    else
        V = π.V[current_player(env)]
        s = state(env)
        I = send_to_device(device(V), Flux.unsqueeze(s, ndims(s) + 1))
        A = action_space(env)
        m = legal_action_space_mask(env)
        σ = masked_regret_matching(V(I) |> send_to_host |> vec, m)
        push!(π.MΠ, I = s, t = π.t, σ = σ, m = m)
        a = sample(π.rng, A, Weights(σ, 1.0))
        env(a)
        external_sampling!(π, env, p)
    end
end

"This is the specific regret matching method used in DeepCFR"
function masked_regret_matching(v, m)
    v⁺ = max.(v .* m, 0.0f0)
    s = sum(v⁺)
    if s > 0
        v⁺ ./= s
    else
        fill!(v⁺, 0.0f0)
        v⁺[findmax(v, m)[2]] = 1.0
    end
    v⁺
end


text:
# by lazarusA # HIDE
#using GLMakie, Random # HIDE
using GLMakie
GLMakie.activate!() #HIDE
let
    x = LinRange(-2, 0.5, 501)
    y = LinRange(-2, 2, 501)
    # objective function
    z = 100 .* (y' .- x .^ 2) .^ 2 .+ (1 .- x) .^ 2
    # build constraints
    zin = copy(z)
    zout1 = copy(z)
    zout2 = copy(z)
    for (i,x) in enumerate(x), (j,y) in enumerate(y)
        if x*y >=1
            zin[i,j] = NaN
            zout1[i,j] = z[i,j]
            zout2[i,j] = NaN
        elseif (x + y^2) >=1
            zin[i,j] = NaN
            zout1[i,j] = NaN
            zout2[i,j] = z[i,j]
        else
            zin[i,j] = z[i,j]
            zout1[i,j] = NaN
            zout2[i,j] = NaN
        end
    end
    fig = Figure(resolution = (1200,900), fontsize = 20)
    ax = Axis3(fig, aspect = (1,1,1), perspectiveness = 0.5, elevation = π/9,
        azimuth = 0.2π, zgridcolor = :grey,ygridcolor=:grey,xgridcolor = :grey)
    pltobj = surface!(ax, x, y, zin, colormap = Reverse(:inferno))
    surface!(ax, x, y, zout1, colormap = (:orangered,:red))
    surface!(ax, x, y, zout2, colormap = (:dodgerblue,:blue))
    cbar = Colorbar(fig, pltobj, label = "xy <=1 & x + y^2<=1",
        height = Relative(0.5), width = 20 )
    fig[1,1] = ax
    fig[1,2] = cbar
    fig
    save(joinpath(@__DIR__, "output", "surfaceConstraints.png"), fig, px_per_unit = 2.0) # HIDE
end
using Pkg # HIDE
Pkg.status(["GLMakie"]) # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie