text:
using FreeSteam

const eps = 1e-8

@static is_apple() ? println("no osx support") : begin

#freesteam

libpath = abspath(joinpath(@__FILE__,"..","..","lib"))
@static if is_linux()
  gsl = Libdl.dlopen(joinpath(libpath, "libgsl.so"));
  gslcblas = Libdl.dlopen(joinpath(libpath, "libgslcblas.so"));
end

function test_region1_point(T::Float64 ,p::Float64 ,v::Float64 ,h::Float64 ,u::Float64 ,s::Float64 ,cp::Float64 ,w::Float64)

  #=
    Note: inputs to this function need units conversion!
      Temperature temp;    ///< K
      Pressure pres;       ///< MPa
      SpecificVolume v;    ///< m�/kg
      SpecificEnergy h;    ///< enthalpy / kJ/kg
      SpecificEnergy u;    ///< internal energy / kJ/kg
      SpecificEntropy s;   ///< kJ/kg.K
      SpecHeatCap cp;      ///< specific heat capacity at constant pressure
      Velocity w;          ///< speed of sound
  =#

	S = freesteam_region1_set_pT(p*1e6, T);
	@test_approx_eq_eps freesteam_p(S)  p*1e6     eps*p*1e6
	@test_approx_eq_eps freesteam_T(S)  T         eps*T
	@test_approx_eq_eps freesteam_v(S)  v         eps*v
	@test_approx_eq_eps freesteam_h(S)  h*1e3     eps*h*1e3
	@test_approx_eq_eps freesteam_u(S)  u*1e3     eps*u*1e3
	@test_approx_eq_eps freesteam_s(S)  s*1e3     eps*s*1e3
	@test_approx_eq_eps freesteam_cp(S) cp*1e3    eps*cp*1e3
	@test_approx_eq_eps freesteam_w(S)  w         eps*w

end

function testregion1()
	println("REGION 1 TESTS")
	test_region1_point(300., 3., 0.100215168E-2, 0.115331273E3, 0.112324818E3, 0.392294792, 0.417301218E1, 0.150773921E4)
	test_region1_point(300., 80., 0.971180894E-3, 0.184142828E3, 0.106448356E3, 0.368563852, 0.401008987E1, 0.163469054E4)
	test_region1_point(500., 3., 0.120241800E-2, 0.975542239E3, 0.971934985E3, 0.258041912E1, 0.465580682E1, 0.124071337E4)
end

function test_region2_point(T::Float64 ,p::Float64 ,v::Float64 ,h::Float64 ,u::Float64 ,s::Float64 ,cp::Float64 ,w::Float64)

	# units of measurement as for region1 test

	S = freesteam_region2_set_pT(p*1e6, T);
	@test_approx_eq_eps freesteam_p(S)  p*1e6     eps*p*1e6
	@test_approx_eq_eps freesteam_T(S)  T         eps*T
	@test_approx_eq_eps freesteam_v(S)  v         eps*v
	@test_approx_eq_eps freesteam_h(S)  h*1e3     eps*h*1e3
	@test_approx_eq_eps freesteam_u(S)  u*1e3     eps*u*1e3
	@test_approx_eq_eps freesteam_s(S)  s*1e3     eps*s*1e3
	@test_approx_eq_eps freesteam_cp(S) cp*1e3    eps*cp*1e3
	@test_approx_eq_eps freesteam_w(S)  w         eps*w

end

function testregion2()
	println("REGION 2 TESTS")
	test_region2_point(300., 0.0035, 0.394913866E2, 0.254991145E4, 0.241169160E4, 0.852238967E1, 0.191300162E1, 0.427920172E3)
	test_region2_point(700., 0.0035, 0.923015898E2, 0.333568375E4, 0.301262819E4, 0.101749996E2, 0.208141274E1, 0.644289068E3)
	test_region2_point(700., 30., 0.542946619E-02, 0.263149474E+4, 0.246861076E+4, 0.517540298E+1, 0.103505092E+2, 0.480386523E+3)
end

function test_region3_point(T::Float64 ,rho::Float64 ,p::Float64 ,h::Float64 ,u::Float64 ,s::Float64 ,cp::Float64 ,w::Float64)

	# units of measurement as for region1 test

	S = freesteam_region3_set_rhoT(rho, T)
	@test_approx_eq_eps freesteam_p(S)  p*1e6     eps*p*1e6
	@test_approx_eq_eps freesteam_T(S)  T         eps*T
	@test_approx_eq_eps freesteam_v(S)  1./rho    eps*1./rho
	@test_approx_eq_eps freesteam_h(S)  h*1e3     eps*h*1e3
	@test_approx_eq_eps freesteam_u(S)  u*1e3     eps*u*1e3
	@test_approx_eq_eps freesteam_s(S)  s*1e3     eps*s*1e3
	@test_approx_eq_eps freesteam_cp(S) cp*1e3    eps*cp*1e3
	@test_approx_eq_eps freesteam_w(S)  w         eps*w

end

function testregion3()
	println("REGION 3 TESTS")
	test_region3_point(650., 500., 0.255837018E2, 0.186343019E4, 0.181226279E4, 0.405427273E1,0.138935717E2, 0.502005554E3)
	test_region3_point(650., 200., 0.222930643E2, 0.237512401E4, 0.226365868E4, 0.485438792E1, 0.446579342E2, 0.383444594E3)
	test_region3_point(750., 500., 0.783095639E2, 0.225868845E4, 0.210206932E4, 0.446971906E1, 0.634165359E1, 0.760696041E3);
end

# --- REGION 4 SATURATION LINE

function test_region4_point(T::Float64 ,p::Float64)
	S = freesteam_region4_set_Tx(T,0.)
	p1 = freesteam_p(S)
	@test_approx_eq_eps p1  p*1e6 eps*p*1e6
	T1 = freesteam_region4_Tsat_p(p1)
	@test_approx_eq_eps T1  T eps*T
end

function testregion4()
	println("REGION 4 TESTS")
	test_region4_point(300.,	0.353658941E-2)
	test_region4_point(500.,	0.263889776E1)
	test_region4_point(600.,	0.123443146E2)
end

# --- REGION 1 BACKWARDS (P,H)

function test_region1_ph_point(p::Float64, h::Float64, T::Float64)
	T1 = freesteam_region1_T_ph(p*1e6,h*1e3)
	@test_approx_eq_eps T1 T T/eps
end

function testregion1ph()
	println("REGION 1 (P,H)")
	test_region1_ph_point(3.,	500.,	0.391798509e3)
	test_region1_ph_point(80.,	500.,	0.378108626e3)
	test_region1_ph_point(80.,	1500.,	0.611041229e3)
end

testregion1()
testregion2()
testregion3()
testregion4()
testregion1ph()

T = 400. # in Kelvin!
p = 1e5  # 1 bar
ss = freesteam_set_pT(p, T)
s = freesteam_s(ss);
@test_approx_eq 7502.40089208754 s # J/kgK
ss2 = freesteam_set_ps(p*10, s)
T2 = freesteam_T(ss2)
@test_approx_eq 684.5051229099 T2
p2 = freesteam_p(ss2)
@test_approx_eq 7502.40089208754 freesteam_s(ss2) # J/kgK

@static if is_linux()
  Libdl.dlclose(gsl);
  Libdl.dlclose(gslcblas);
end

end


text:
" Generate datasets "
# using JLD
using Random, Distributions

# N = 50
# d = 5
# n = 10
# k = 2

# zero error
# A = zeros(d, N)
# U = rand(d, k)
# V = rand(k, n)
# A[:, 1:n] = U*V
# A[:, n+1:N] = rand(d, N-n)
# save("datasets/N_$(N)_d_$(d).jld", "data", collect(eachcol(A)), "n_$(n)_k_$(k)", (0.0, collect(1:n)))

"""
generate N×d matrix A where the upper left p×n submatrix of A is of rank k, the rest Uniform[0,1] noise
"""
function matrix_with_subspace_noise(N, d, n, k, p)
    @assert p <= d
    @assert k < p
    A = rand(Uniform(-1.0, 1.0), d, N)
    U = rand(Uniform(-1.0, 1.0), p, k)
    V = rand(Uniform(-1.0, 1.0), k, n)
    A[1:p, 1:n] = U*V
    return A
end
# save("datasets/N_$(N)_d_$(d)_nonZeroError.jld", "data", collect(eachcol(A)), "n_$(n)_k_$(k)", (err, collect(1:n)))

"""
generate N×d matrix (A + M) where the the left d×n submatrix of A is of rank k,
the rest Uniform[0,1] noise, and M being a noise matrix with entries in Normal(0, std)
"""
function matrix_with_normal_noise(N, d, n, k, std)
    A = rand(Uniform(-1.0, 1.0), d, N)
    U = rand(Uniform(-1.0, 1.0), d, k)
    V = rand(Uniform(-1.0, 1.0), k, n)
    A[:, 1:n] = U*V
    noise = rand(Normal(0, std), d, N)
    return A + noise
end

ext: jl


text:
# Produces n copies of a density matrix

export copies

""" `rhoOut = copies(rho,n)`

Returns *rho*^(tensor *n*), i.e., *n* copies of the input *rho*
"""

function copies(rho::AbstractMatrix, n::Int)
	return reduce(kron, repeated(rho, n))
end

function copies(rho::AbstractMatrix, n::Float64)
	return copies(rho, convert(Int, n))
end


text:
"""
# module Cards

# Examples

```jldoctest
julia>
```
"""
module Cards

using Parsso.CoreDataStructure.Elements

export Card,
    OptionedCard

abstract type AbstractCard end

struct Card <: AbstractCard
    name::String
    values
end

struct OptionedCard <: AbstractCard
    name::String
    option::String
    values
end

end


text:
using DiffPrecTest
using Statistics, StatsBase, LinearAlgebra
using SparseArrays
using ProximalBase, CoordinateDescent, CovSel
using Random, Distributions
using JLD

@show gethostname()

rep   = parse(Int,ARGS[1])
ip    = parse(Int,ARGS[2])
dir = ARGS[3]

pArr = [50, 100, 150]
n = 300

p = pArr[ip]
Random.seed!(54298)

# generate model
Ωx = Matrix{Float64}(I, p, p)
for l=1:p-1
  Ωx[l  , l+1] = 0.6
  Ωx[l+1, l  ] = 0.6
end
for l=1:p-2
  Ωx[l  , l+2] = 0.3
  Ωx[l+2, l  ] = 0.3
end

# generate Delta
d = Vector{Float64}(undef, p)
rand!(Uniform(0.5, 2.5), d)
d .= sqrt.(d)
D = Diagonal(d)
Σ = inv(Symmetric(D * Ωx * D))

dist_X = MvNormal(convert(Matrix, Σ))

# generate data
Random.seed!(12346 + rep)
X = rand(dist_X, n)'
Y = rand(dist_X, n)'

Sx = Symmetric( X'X / n )
Sy = Symmetric( Y'Y / n )

eS_init = Array{BitArray}(undef, div((p + 1)*p, 2))
for j=1:div((p + 1)*p, 2)
    eS_init[j] = falses(p, p)
end

@time boot_res, eS = bootstrap(X, Y; estimSupport=eS_init)

@save "$(dir)/res_$(ip)_$(rep).jld" boot_res


text:
function dgpnorm(n::Int, x::Array{Float64}, xl::Array{Float64},
                 xu::Array{Float64}, g::Array{Float64})
    inf_norm = 0.0
    for i=1:n
        if xl[i] != xu[i]
            if x[i] == xl[i]
                v = (min(g[i], 0.0))^2
            elseif x[i] == xu[i]
                v = (max(g[i], 0.0))^2
            else
                v = g[i]^2
            end

            v = sqrt(v)
            inf_norm = (inf_norm > v) ? inf_norm : v
        end
    end

    return inf_norm
end

@inline function dgpnorm(n::Int, x::CuDeviceArray{Float64,1}, xl::CuDeviceArray{Float64,1},
                         xu::CuDeviceArray{Float64,1}, g::CuDeviceArray{Float64,1})
    tx = threadIdx().x

    v = 0.0
    if tx <= n
        @inbounds begin
            if xl[tx] != xu[tx]
                if x[tx] == xl[tx]
                    v = min(g[tx], 0.0)
                elseif x[tx] == xu[tx]
                    v = max(g[tx], 0.0)
                else
                    v = g[tx]
                end

                v = abs(v)
            end
        end
    end

    # shfl_down_sync() will automatically sync threads in a warp.

    offset = 16
    while offset > 0
        v = max(v, CUDA.shfl_down_sync(0xffffffff, v, offset))
        offset >>= 1
    end
    v = CUDA.shfl_sync(0xffffffff, v, 1)

    return v
end


text:
#= Run as:
cd(".../Modia/test")
include("runsymbolictests.jl")
=#

# Desired:
#   using Test
#
# In order that Test need not to be defined in the user environment, it is included via Modia:
using Modia

@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Modia.Test
end

@testset "SymbolicTests" begin

    include("symbolic/BLTandPantelides/setup.jl")
    include("symbolic/BLTandPantelides/testBLTandPantelides.jl")

    #  include("symbolic/DAEquations/setup.jl")
    include("symbolic/DAEquations/testSymbolicTransform.jl")

    include("symbolic/TestExactlyRemoveSingularities.jl")

    include("symbolic/TestStateSelectionAlgorithm.jl")

    include("symbolic/TestTearingAlgorithm.jl")

end
nothing


text:
module TestLQNash

using QuantEcon
using Base.Test
using FactCheck

# set up for when agents don't interact with each other
a = [.95 0.
     0 .95]
b1 = [.95; 0.]
b2 = [0.; .95]
r1 = [-.25 0.
      0. 0.]
r2 = [0. 0.
      0. -.25]
q1 = [-.15]
q2 = [-.15]

f1, f2, p1, p2 = nnash(a, b1, b2, r1, r2, q1, q2, 0, 0, 0, 0, 0, 0,
                       tol=1e-8, max_iter=10000)

alq = .95
blq = .95
rlq = -.25
qlq = -.15

lq_obj = LQ(qlq, rlq, alq, blq, bet=1)

p, f, d = stationary_values(lq_obj)

facts("Testing lqnash.jl") do

    context("Checking the policies") do

        @fact sum(f1) --> roughly(sum(f2))
        @fact sum(f1) --> roughly(sum(f))
    end

    context("Checking the Value Function") do

        @fact p1[1, 1] --> roughly(p2[2, 2])
        @fact p1[1, 1] --> roughly(p[1])
    end

    context("Judd test case")  do
        # Define Parameters
        delta = 0.02
        d = [-1 0.5
           0.5 -1]
        B = [25; 25]
        c1 = [1; -2; 1]
        c2 = [1; -2; 1]
        e1 = [10; 10; 3]
        e2 = [10; 10; 3]
        delta_1 = 1 - delta

        ## Define matrices
        a = [delta_1 0 -delta_1*B[1]
           0 delta_1 -delta_1*B[2]
           0 0 1]

        b1 = delta_1 * [1 -d[1, 1]
                      0 -d[2, 1]
                      0 0]
        b2 = delta_1 * [0 -d[1, 2]
                      1 -d[2, 2]
                      0 0]

        r1 = -[0.5*c1[3] 0 0.5*c1[2]
             0 0 0
             0.5*c1[2] 0 c1[1]]
        r2 = -[0 0 0
             0 0.5*c2[3] 0.5*c2[2]
             0 0.5*c2[2] c2[1]]

        q1 = [-0.5*e1[3] 0
            0 d[1, 1]]
        q2 = [-0.5*e2[3] 0
            0 d[2, 2]]

        s1 = zeros(2, 2)
        s2 = copy(s1)

        w1 = [0 0
            0 0
           -0.5*e1[2] B[1]/2.]
        w2 = [0 0
            0 0
           -0.5*e2[2] B[2]/2.]

        m1 = [0 0
            0 d[1, 2]/2.]
        m2 = copy(m1)

        # build model and solve it
        f1, f2, p1, p2 = nnash(a, b1, b2, r1, r2, q1, q2, s1, s2, w1, w2, m1, m2)

        aaa = a - b1*f1 - b2*f2
        aa = aaa[1:2, 1:2]
        tf = eye(2) - aa
        tfi = inv(tf)
        xbar = tfi*aaa[1:2, 3]

        # Define answers from matlab. TODO: this is ghetto
        f1_ml = [0.243666582208565   0.027236062661951 -6.827882928738190
                 0.392370733875639   0.139696450885998 -37.734107291009138]

        f2_ml = [0.027236062661951   0.243666582208565  -6.827882928738186
                 0.139696450885998   0.392370733875639 -37.734107291009131]

        xbar_ml = [1.246871007582702, 1.246871007582685]

        @fact f1 --> roughly(f1_ml)
        @fact f2 --> roughly(f2_ml)
        @fact xbar --> roughly(xbar_ml)

    end
end
end  # module


text:
"""
    trig()
``trig'' is an MINLPLib test instance with trigonometric functions
(see http://www.minlplib.org/trig.html). The best known primal feasible
objective value for this problem is -3.76250036. This example illustrates
how to obtain lower bounds to the optimal objective value using the MILP and
LP relaxations obtained using the PolyhedralRelaxations package.
trig is a non-linear non-convex minimization problem with one decision variable
and one inequality constraint.
Variables: -2 <= x <= 5; start = 1
Constraint: 5 * sin(x) - x <= 0.0
Objective: Minimize sin(11 * x) - cos(13 * x) + sin(17 * x) + cos(19 * x)
We reformulate this problem using additional variables as
Variables:
    -2 <= x <= 5; start = 1
    -1 <= y[1:5] <= 1
Constraints:
    5 * y[1] - x <= 0.0
    y[1] = sin(x)
    y[2] = sin(11 * x)
    y[3] = cos(13 * x)
    y[4] = sin(17 * x)
    y[5] = cos(19 * x)
Objective: Minimize y[2] + y[3] - y[4] - y[5]
"""

function trig()
    best_known_objective = -3.76250036
    cplex_optimizer =
        JuMP.optimizer_with_attributes(CPLEX.Optimizer, "CPX_PARAM_SCRIND" => 0)
    error_tolerances = [NaN64, 1e-1, 1e-2]
    # base partition holds the inflection points and the end-points for the trigonometric functions
    milp_statistics = zeros(length(error_tolerances), 5)
    lp_statistics = zeros(length(error_tolerances), 5)
    base_partition = Dict{Int,Vector{Float64}}()
    base_partition[1] = [-2.0, 0.0, π, 5.0]
    base_partition[2] = [-2.0, collect(-7*π/11:π/11:17*π/11)..., 5.0]
    base_partition[3] = [-2.0, collect(-15*π/26:π/13:41*π/26)..., 5.0]
    base_partition[4] = [-2.0, collect(-10*π/17:π/17:27*π/17)..., 5.0]
    base_partition[5] = [-2.0, collect(-23*π/38:π/19:59*π/38)..., 5.0]
    # functions
    functions = Dict{Int,Function}()
    functions[1] = x -> sin(x)
    functions[2] = x -> sin(11 * x)
    functions[3] = x -> cos(13 * x)
    functions[4] = x -> sin(17 * x)
    functions[5] = x -> cos(19 * x)
    # solve MILP for various error tolerances
    for i = 1:length(error_tolerances)
        error_tolerance = error_tolerances[i]
        milp_var = Dict{Int,Any}()
        x_indexes = Dict{Int,Int}()
        y_indexes = Dict{Int,Int}()
        # construct MILP relaxations
        milp = Model(cplex_optimizer)
        for j = 1:5
            var, x_index, y_index = add_milp_relaxation(
                functions[j],
                base_partition[j],
                error_tolerance,
                milp,
                name = string(j),
            )
            milp_var[j] = var
            x_indexes[j] = x_index
            y_indexes[j] = y_index
            (j == 1) && (continue)
            @constraint(milp, var[x_index] == milp_var[1][x_indexes[1]])
        end
        @constraint(milp, 5 * milp_var[1][y_indexes[1]] - milp_var[1][x_indexes[1]] <= 0.0)
        @objective(
            milp,
            Min,
            milp_var[2][y_indexes[2]] + milp_var[3][y_indexes[3]] -
            milp_var[4][y_indexes[4]] - milp_var[5][y_indexes[5]]
        )
        _, solve_time, solve_bytes_alloc, sec_in_gc = @timed optimize!(milp)
        relaxation_objective = objective_value(milp)
        @assert relaxation_objective ≤ best_known_objective
        relative_gap =
            abs(best_known_objective - relaxation_objective) / abs(relaxation_objective)
        milp_statistics[i, :] = [
            error_tolerance,
            best_known_objective,
            relaxation_objective,
            relative_gap,
            round(solve_time, digits = 2),
        ]
    end

    for i = 1:length(error_tolerances)
        error_tolerance = error_tolerances[i]
        lp_var = Dict{Int,Any}()
        x_indexes = Dict{Int,Int}()
        y_indexes = Dict{Int,Int}()
        # construct LP relaxations
        lp = Model(cplex_optimizer)
        for j = 1:5
            var, x_index, y_index = add_lp_relaxation(
                functions[j],
                base_partition[j],
                error_tolerance,
                lp,
                name = string(j),
            )
            lp_var[j] = var
            x_indexes[j] = x_index
            y_indexes[j] = y_index
            (j == 1) && (continue)
            @constraint(lp, var[x_index] == lp_var[1][x_indexes[1]])
        end
        @constraint(lp, 5 * lp_var[1][y_indexes[1]] - lp_var[1][x_indexes[1]] <= 0.0)
        @objective(
            lp,
            Min,
            lp_var[2][y_indexes[2]] + lp_var[3][y_indexes[3]] - lp_var[4][y_indexes[4]] -
            lp_var[5][y_indexes[5]]
        )
        _, solve_time, solve_bytes_alloc, sec_in_gc = @timed optimize!(lp)
        relaxation_objective = objective_value(lp)
        @assert relaxation_objective ≤ best_known_objective
        relative_gap =
            abs(best_known_objective - relaxation_objective) / abs(relaxation_objective)
        lp_statistics[i, :] = [
            error_tolerance,
            best_known_objective,
            relaxation_objective,
            relative_gap,
            round(solve_time, digits = 2),
        ]
    end

    open("csv/trig_milp.csv", "w") do io
        writedlm(io, ["error_tol" "g_opt" "relaxation" "relative_gap" "time"], ',')
        writedlm(io, milp_statistics, ",")
    end

    open("csv/trig_lp.csv", "w") do io
        writedlm(io, ["error_tol" "g_opt" "relaxation" "relative_gap" "time"], ',')
        writedlm(io, lp_statistics, ",")
    end

end

max_stars_repo_name: sujeevraja/univariate-experiments

max_issues_repo_name: sujeevraja/univariate-experiments

max_forks_repo_name: sujeevraja/univariate-experiments


text:
using Dates

@testset "SO6.jl" begin
    filename = "data\\test.so6"
    df = DDR2import.SO6.read(filename)
    @test df[2,:SEGMENT_ID] == "*AM80_!AAEW"
    @test df[2,:ADEP] == "EHAM"
    @test df[2,:ADES] == "UUEE"
    @test df[2,:ACTYPE] == "A321"
    @test df[2,:FLBEGINSEGMENT] == 10
    @test df[2,:FLENDSEGMENT] == 25
    @test df[2,:STATUS] == 0
    @test df[2,:CALLSIGN] == "AFL2193"
    @test df[2,:LATBEGINSEGMENT_DEG] == 3137.483333 / 60.0
    @test df[2,:LONBEGINSEGMENT_DEG] == 284.183333 / 60.0
    @test df[2,:LATENDSEGMENT_DEG] == 3135.833333 / 60.0
    @test df[2,:LONENDSEGMENT_DEG] == 284.033333 / 60.0
    @test df[2,:FLIGHT_ID] == 213765625
    @test df[2,:SEQUENCE] == 2
    @test df[2,:SEGMENT_LENGTH_M] == 1.652550 * 1852.0
    @test df[2,:SEGMENT_PARITY] == 0
    @test Dates.year(df[2,:DATETIMEBEGINSEGMENT]) == 2017
    @test Dates.hour(df[2,:DATETIMEBEGINSEGMENT]) == 21
    @test Dates.second(df[2,:DATETIMEBEGINSEGMENT]) == 26
    @test Dates.day(df[2,:DATETIMEENDSEGMENT]) == 14
    @test Dates.minute(df[2,:DATETIMEENDSEGMENT]) == 50
    #LAST
    @test df[3,:SEGMENT_ID] == "!AAEW_!AAEX"
end


text:
# -----------------------------------------------------------------------------
# Integration with AtomsBase
# -----------------------------------------------------------------------------

# Internal type used to represent a MassBody and retain the element metadata
struct ElementMassBody{cType<:Real,mType<:Real} <: Body
    r::SVector{3,cType} # in LENGTH_UNIT
    v::SVector{3,cType} # in VELOCITY_UNIT
    m::mType            # in MASS_UNIT
    symbol::Symbol
    number::Int
    data::Dict{Symbol,Any}
end
function ElementMassBody(r::SVector{3,<:Unitful.Length}, v::SVector{3,<:Unitful.Velocity}, e::Element; data...)
    ElementMassBody{Float64,Float64}(austrip.(r), austrip.(v), austrip(e.atomic_mass), Symbol(e.symbol), e.number, Dict(data...))
end
function ElementMassBody(body::ElementMassBody{cType,mType}, r::SVector{3,cType}, v::SVector{3,cType}) where {cType<:Real,mType<:Real}
    ElementMassBody{cType,mType}(r, v, body.m, body.symbol, body.number, body.data)
end

# Convert AtomsBase Atom to NBodySimulator body
function ElementMassBody(atom::Atom)
    ElementMassBody{Float64,Float64}(
        austrip.(position(atom)),
        austrip.(velocity(atom)),
        austrip(atomic_mass(atom)),
        AtomsBase.atomic_symbol(atom),
        atomic_number(atom),
        atom.data
    )
end
# Convert NBodySimulator body to AtomsBase Atom
# TODO: support more boundary conditions
function AtomsBase.Atom(b::ElementMassBody, boundary_conditions::CubicPeriodicBoundaryConditions)
    Atom(b.symbol, mod.(b.r, boundary_conditions.L) .* LENGTH_UNIT, b.v .* VELOCITY_UNIT; b.data...)
end

# Convert AtomsBase boundary conditions to NBodySimulator boundary conditions
function nbs_boundary_conditions(system::AbstractSystem{3})
    # TODO: support more boundary conditions
    box = bounding_box(system)
    @assert all(ustrip(box[i][j]) == 0 for i ∈ 1:3 for j ∈ 1:3 if i != j)
    @assert box[1][1] == box[2][2] == box[3][3]
    @assert boundary_conditions(system) == [Periodic(), Periodic(), Periodic()]
    CubicPeriodicBoundaryConditions(austrip(box[1][1]))
end
# Convert NBodySimulator boundary conditions to AtomsBase boundary conditions
# TODO: support more boundary conditions
ab_boundary_conditions(::CubicPeriodicBoundaryConditions) = [Periodic(), Periodic(), Periodic()]

# Convert NBodySimulator boundary conditions to AtomsBase bounding box
# TODO: support more boundary conditions
function ab_bounding_box(boundary_conditions::CubicPeriodicBoundaryConditions)
    box_size = boundary_conditions.L * LENGTH_UNIT
    z = zero(typeof(box_size))
    [[box_size, z, z], [z, box_size, z], [z, z, box_size]]
end

# Convert AtomsBase AbstractSystem to Vector of NBodySimulator bodies
bodies(system::AbstractSystem{3}) = ElementMassBody.(system)
# Convert Vector of NBodySimulator bodies to AtomsBase FlexibleSystem
function AtomsBase.FlexibleSystem(bodies::AbstractVector{<:ElementMassBody}, boundary_conditions::BoundaryConditions)
    particles = Fix2(Atom, boundary_conditions).(bodies)
    FlexibleSystem(particles, ab_bounding_box(boundary_conditions), ab_boundary_conditions(boundary_conditions))
end
# Convert Vector of NBodySimulator bodies to AtomsBase FastSystem
function AtomsBase.FastSystem(bodies::AbstractVector{<:ElementMassBody}, boundary_conditions::BoundaryConditions)
    particles = Fix2(Atom, boundary_conditions).(bodies)
    FastSystem(particles, ab_bounding_box(boundary_conditions), ab_boundary_conditions(boundary_conditions))
end

# -----------------------------------------------------------------------------
# Convenience methods for generating starting configurations with element data
# -----------------------------------------------------------------------------

function NBodySimulator.generate_bodies_in_cell_nodes(n::Integer, symbol::Union{Integer,AbstractString,Symbol}, L::Unitful.Length, reference_temp::Unitful.Temperature; rng = MersenneTwister(n))
    e = elements[symbol]
    average_velocity = √(u"k" * reference_temp / e.atomic_mass)
    generate_bodies_in_cell_nodes(n, symbol, average_velocity, L, rng = rng)
end
function NBodySimulator.generate_bodies_in_cell_nodes(n::Integer, symbol::Union{Integer,AbstractString,Symbol}, average_velocity::Unitful.Velocity, L::Unitful.Length; rng = MersenneTwister(n))
    velocities = average_velocity * randn(rng, Float64, (3, n))
    e = elements[symbol]
    bodies = ElementMassBody[]

    count = 1
    dL = L / (ceil(n^(1 / 3)))
    for x ∈ dL/2:dL:L, y ∈ dL/2:dL:L, z ∈ dL/2:dL:L
        if count > n
            break
        end
        push!(bodies, ElementMassBody(SVector(x, y, z), SVector{3}(velocities[:, count]), e))
        count += 1
    end
    return bodies
end


text:
using Distributions: rand
using BenchmarkTools: length, Benchmark, Trial
using BenchmarkTools
using Distributions
using LinearAlgebra
using BenchmarkPlots, StatsPlots

include("mydct.jl")
global res
global suite
global sort_res
global sort_dct2
global sort_slow

function runTest(start=10, stop=100)
    square_dimensions = collect(range(start, stop, step=5))
    n_dimensions = length(square_dimensions)

    suite = BenchmarkGroup()
    rangearr = collect(range(1, n_dimensions, step=1))
    suite["dct2"] = BenchmarkGroup(sort(rangearr))
    suite["fast"] = BenchmarkGroup(sort(rangearr))
    suite["slow"] = BenchmarkGroup(sort(rangearr))
    # suite["slow"] = BenchmarkGroup(sort(rangearr))

    param = Exponential(100 * rand())
    for i = 1:n_dimensions
        m = rand(param, square_dimensions[i], square_dimensions[i])
        suite["dct2"][i] = @benchmarkable mDCT2n3($m)  evals=1 samples=5
        suite["fast"][i] = @benchmarkable dct($m) evals=1 samples=5
        if (mod(square_dimensions[i], 15) == 0)
            suite["slow"][i] = @benchmarkable mDCT2n4($m) evals=1 samples=1
        end
    end

    # tune!(suite)
    global suite

    global res = run(suite; verbose=true)
    global sort_res = sort(collect(mean(res["fast"])), by = x->x[1])
    global sort_dct2 = sort(collect(mean(res["dct2"])), by = x->x[1])
    global sort_slow = sort(collect(mean(res["slow"])), by = x->x[1])

end

# i=1
# val = zeros(Int, n_dimensions)
# time_s = zeros(n_dimensions)
# for p in sort_res
#     val[i], time_s[i] = p
#     i = i+1
# end

max_stars_repo_name: illoxian/bitmap-compressor-tool

max_issues_repo_name: illoxian/bitmap-compressor-tool

max_forks_repo_name: illoxian/bitmap-compressor-tool


text:
mutable struct Haploid <: AbstractAgent
    id::Int
    trait::Float64
end
function sampleModelInitialize(;n = 100)
  sampleModel = ABM(Haploid)
  for i in 1:n
      add_agent!(sampleModel, rand())
  end
  return sampleModel
end

mutable struct HaploidS <: AbstractAgent
    id::Int
    pos::Tuple{Int, Int}
    trait::Float64
end
function sampleModelInitializeS(;n = 100, dims=(2,3))
  sampleModel = ABM(HaploidS, Space(dims))
  for i in 1:n
      add_agent!(sampleModel, rand())
  end
  return sampleModel
end


text:
abstract type AbstractNLsolveSolver end
abstract type AbstractNLsolveCache end
mutable struct NLSolverCache{rateType,uType,W,uToltype,cType,gType} <: AbstractNLsolveCache
  κ::uToltype
  tol::uToltype
  min_iter::Int
  max_iter::Int
  nl_iters::Int
  new_W::Bool
  z::uType
  W::W # NLNewton -> `W` operator; NLAnderson -> Vectors; NLFunctional -> Nothing
  γ::gType
  c::cType
  ηold::uToltype
  # The following fields will alias for immutable cache
  z₊::uType # Only used in `NLAnderson` and `NLFunctional`
  dz::uType
  tmp::uType
  b::uType # can be aliased with `k` if no unit
  k::rateType
end

struct NLFunctional{iip} <: AbstractNLsolveSolver
  cache::NLSolverCache
end
struct NLAnderson{iip} <: AbstractNLsolveSolver
  cache::NLSolverCache
  n::Int
end
struct NLNewton{iip} <: AbstractNLsolveSolver
  cache::NLSolverCache
end

NLSolverCache(;κ=nothing, tol=nothing, min_iter=1, max_iter=10) =
NLSolverCache(κ, tol, min_iter, max_iter, 0, true,
              (nothing for i in 1:10)...)

# Default `iip` to `true`, but the whole type will be reinitialized in `alg_cache`
NLFunctional(;kwargs...) = NLFunctional{true}(NLSolverCache(;kwargs...))
NLAnderson(n=5; kwargs...) = NLAnderson{true}(NLSolverCache(;kwargs...), n)
NLNewton(;kwargs...) = NLNewton{true}(NLSolverCache(;kwargs...))

oop_nlsolver(s::NLFunctional{true}) = NLFunctional{false}(s.cache)
oop_nlsolver(s::NLAnderson{true}) = NLAnderson{false}(s.cache, s.n)
oop_nlsolver(s::NLNewton{true}) = NLNewton{false}(s.cache)


text:
using CUDA
using GemmKernels

M = parse(Int, ARGS[1])
N = parse(Int, ARGS[2])
K = parse(Int, ARGS[3])

transpose_a = (ARGS[4] == "n" ? false : ARGS[4] == "t" ? true : error("Invalid memory layout for A: $(ARGS[4])"))
transpose_b = (ARGS[5] == "n" ? false : ARGS[5] == "t" ? true : error("Invalid memory layout for B: $(ARGS[5])"))

function benchmark_matmul(a, b, c, d)
    CUDA.@sync begin
        conf = GemmKernels.get_config(
                gemm_shape = (M = M, N = N, K = K),
                operator = Operator.WMMADualOp{16, 16, 16},

                global_a_layout = transpose_a ? Layout.InterleavedRowMajor{Float16} : Layout.InterleavedColMajor{Float16},
                global_b_layout = transpose_b ? Layout.InterleavedRowMajor{Float16} : Layout.InterleavedColMajor{Float16},
                global_c_layout = Layout.InterleavedColMajor{Float32},
                global_d_layout = Layout.InterleavedColMajor{Float32},

                shared_a_layout = transpose_a ? Layout.Padded{Layout.SplitRowMajor{Float16}, 8} : Layout.Padded{Layout.SplitColMajor{Float16}, 8},
                shared_b_layout = transpose_b ? Layout.Padded{Layout.SplitRowMajor{Float16}, 8} : Layout.Padded{Layout.SplitColMajor{Float16}, 8},
                shared_c_layout = Layout.SplitColMajor{Float32},
                shared_d_layout = Layout.SplitColMajor{Float32},

                warps_per_block = 8,

                compute_warp = (M = 16, N = 32, K = 16),

                block_shape = (M = 64, N = 64, K = 32),

                mem_a_warp = transpose_a ? (M = 4, K = 32) : (M = 64, K = 2),
                mem_b_warp = transpose_b ? (K = 2, N = 64) : (K = 32, N = 4),
                mem_cd_warp = (M = 64, N = 1),

                mem_a_thread = transpose_a ? (M = 1, K = 4) : (M = 4, K = 1),
                mem_b_thread = transpose_b ? (K = 1, N = 4) : (K = 4, N = 1),
                mem_cd_thread = (M = 2, N = 1),

                is_a_col_major = !transpose_a,
                is_b_col_major = !transpose_b
            )

        GemmKernels.matmul(a, b, c, d, conf;
                           kernel = Kernel.matmul_pipelined)
    end
end

a_h = rand(Complex{Float16}, (M, K)) / sqrt(Float16(K))
b_h = rand(Complex{Float16}, (K, N)) / sqrt(Float16(K))
c_h = rand(Complex{Float32}, (M, N))

# Transpose input if necessary
a_h = transpose_a ? transpose(a_h) : a_h
b_h = transpose_b ? transpose(b_h) : b_h

a   = CuArray(a_h)
b   = CuArray(b_h)
c   = CuArray(c_h)
d   = similar(c)

# warmup
benchmark_matmul(a, b, c, d)

# profile
for i = 1 : 10
    CUDA.@profile benchmark_matmul(a, b, c, d)
end


text:
module LP

import ..FileFormats
import MathOptInterface
const MOI = MathOptInterface

# Julia 1.6 removes Grisu from Base. Previously, we went
#   print_shortest(io, x) = Base.Grisu.print_shortest(io, x)
# To avoid adding Grisu as a dependency, use the following printing heuristic.
# TODO(odow): consider printing 1.0 as 1.0 instead of 1, i.e., without the
# rounding branch.
function print_shortest(io::IO, x::Real)
    x_int = round(Int, x)
    if isapprox(x, x_int)
        print(io, x_int)
    else
        print(io, x)
    end
    return
end

MOI.Utilities.@model(
    Model,
    (MOI.ZeroOne, MOI.Integer),
    (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan, MOI.Interval),
    (),
    (MOI.SOS1, MOI.SOS2),
    (),
    (MOI.ScalarAffineFunction,),
    (MOI.VectorOfVariables,),
    ()
)
function MOI.supports(
    ::Model{T},
    ::MOI.ObjectiveFunction{<:MOI.ScalarQuadraticFunction{T}},
) where {T}
    return false
end

struct Options
    maximum_length::Int
    warn::Bool
end

function get_options(m::Model)
    default_options = Options(255, false)
    return get(m.ext, :LP_OPTIONS, default_options)
end

"""
    Model(; kwargs...)

Create an empty instance of FileFormats.LP.Model.

Keyword arguments are:

 - `maximum_length::Int=255`: the maximum length for the name of a variable.
   lp_solve 5.0 allows only 16 characters, while CPLEX 12.5+ allow 255.

 - `warn::Bool=false`: print a warning when variables or constraints are renamed.
"""
function Model(; maximum_length::Int = 255, warn::Bool = false)
    model = Model{Float64}()
    options = Options(maximum_length, warn)
    model.ext[:LP_OPTIONS] = options
    return model
end

function Base.show(io::IO, ::Model)
    print(io, "A .LP-file model")
    return
end

# ==============================================================================
#
#   Base.write
#
# ==============================================================================

const START_REG = r"^([\.0-9eE])"
const NAME_REG = r"([^a-zA-Z0-9\!\"\#\$\%\&\(\)\/\,\.\;\?\@\_\`\'\{\}\|\~])"

function write_function(
    io::IO,
    model::Model,
    func::MOI.VariableIndex,
    variable_names::Dict{MOI.VariableIndex,String},
)
    print(io, variable_names[func])
    return
end

function write_function(
    io::IO,
    model::Model,
    func::MOI.ScalarAffineFunction{Float64},
    variable_names::Dict{MOI.VariableIndex,String},
)
    is_first_item = true
    if !(func.constant ≈ 0.0)
        print_shortest(io, func.constant)
        is_first_item = false
    end
    for term in func.terms
        if !(term.coefficient ≈ 0.0)
            if is_first_item
                print_shortest(io, term.coefficient)
                is_first_item = false
            else
                print(io, term.coefficient < 0 ? " - " : " + ")
                print_shortest(io, abs(term.coefficient))
            end

            print(io, " ", variable_names[term.variable])
        end
    end
    return
end

function write_constraint_suffix(io::IO, set::MOI.LessThan)
    print(io, " <= ")
    print_shortest(io, set.upper)
    println(io)
    return
end

function write_constraint_suffix(io::IO, set::MOI.GreaterThan)
    print(io, " >= ")
    print_shortest(io, set.lower)
    println(io)
    return
end

function write_constraint_suffix(io::IO, set::MOI.EqualTo)
    print(io, " = ")
    print_shortest(io, set.value)
    println(io)
    return
end

function write_constraint_suffix(io::IO, set::MOI.Interval)
    print(io, " <= ")
    print_shortest(io, set.upper)
    println(io)
    return
end

function write_constraint_prefix(io::IO, set::MOI.Interval)
    print_shortest(io, set.lower)
    print(io, " <= ")
    return
end

write_constraint_prefix(io::IO, set) = nothing

function write_constraint(
    io::IO,
    model::Model,
    index::MOI.ConstraintIndex,
    variable_names::Dict{MOI.VariableIndex,String};
    write_name::Bool = true,
)
    func = MOI.get(model, MOI.ConstraintFunction(), index)
    set = MOI.get(model, MOI.ConstraintSet(), index)
    if write_name
        print(io, MOI.get(model, MOI.ConstraintName(), index), ": ")
    end
    write_constraint_prefix(io, set)
    write_function(io, model, func, variable_names)
    return write_constraint_suffix(io, set)
end

const SCALAR_SETS = (
    MOI.LessThan{Float64},
    MOI.GreaterThan{Float64},
    MOI.EqualTo{Float64},
    MOI.Interval{Float64},
)

function write_sense(io::IO, model::Model)
    if MOI.get(model, MOI.ObjectiveSense()) == MOI.MAX_SENSE
        println(io, "maximize")
    else
        println(io, "minimize")
    end
    return
end

function write_objective(
    io::IO,
    model::Model,
    variable_names::Dict{MOI.VariableIndex,String},
)
    print(io, "obj: ")
    obj_func_type = MOI.get(model, MOI.ObjectiveFunctionType())
    obj_func = MOI.get(model, MOI.ObjectiveFunction{obj_func_type}())
    write_function(io, model, obj_func, variable_names)
    println(io)
    return
end

"""
    Base.write(io::IO, model::FileFormats.LP.Model)

Write `model` to `io` in the LP file format.
"""
function Base.write(io::IO, model::Model)
    options = get_options(model)
    FileFormats.create_unique_names(
        model,
        warn = options.warn,
        replacements = [
            s -> match(START_REG, s) !== nothing ? "_" * s : s,
            s -> replace(s, NAME_REG => "_"),
            s -> s[1:min(length(s), options.maximum_length)],
        ],
    )
    variable_names = Dict{MOI.VariableIndex,String}(
        index => MOI.get(model, MOI.VariableName(), index) for
        index in MOI.get(model, MOI.ListOfVariableIndices())
    )
    write_sense(io, model)
    write_objective(io, model, variable_names)
    println(io, "subject to")
    for S in SCALAR_SETS
        for index in MOI.get(
            model,
            MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{Float64},S}(),
        )
            write_constraint(
                io,
                model,
                index,
                variable_names;
                write_name = true,
            )
        end
    end
    println(io, "Bounds")
    free_variables = Set(keys(variable_names))
    for S in SCALAR_SETS
        for index in
            MOI.get(model, MOI.ListOfConstraintIndices{MOI.VariableIndex,S}())
            delete!(free_variables, MOI.VariableIndex(index.value))
            write_constraint(
                io,
                model,
                index,
                variable_names;
                write_name = false,
            )
        end
    end
    for index in MOI.get(
        model,
        MOI.ListOfConstraintIndices{MOI.VariableIndex,MOI.ZeroOne}(),
    )
        delete!(free_variables, MOI.VariableIndex(index.value))
    end
    for variable in sort(collect(free_variables), by = x -> x.value)
        println(io, variable_names[variable], " free")
    end
    for (S, str_S) in [(MOI.Integer, "General"), (MOI.ZeroOne, "Binary")]
        indices =
            MOI.get(model, MOI.ListOfConstraintIndices{MOI.VariableIndex,S}())
        if length(indices) > 0
            println(io, str_S)
            for index in indices
                write_function(
                    io,
                    model,
                    MOI.get(model, MOI.ConstraintFunction(), index),
                    variable_names,
                )
                println(io)
            end
        end
    end
    println(io, "End")
    return
end

# ==============================================================================
#
#   Base.read!
#
# ==============================================================================
const _COMMENT_REG = r"(.*?)\\(.*)"
const _READ_START_REG = r"^([\.0-9])"
function _strip_comment(line::String)
    if occursin("\\", line)
        m = match(_COMMENT_REG, line)
        return strip(String(m[1]))
    else
        return strip(line)
    end
end

# a list of section keywords in lower-case
const _KEYWORDS = Dict(
    "max" => Val{:obj},
    "maximize" => Val{:obj},
    "maximise" => Val{:obj},
    "maximum" => Val{:obj},
    "min" => Val{:obj},
    "minimize" => Val{:obj},
    "minimise" => Val{:obj},
    "minimum" => Val{:obj},
    "subject to" => Val{:constraints},
    "such that" => Val{:constraints},
    "st" => Val{:constraints},
    "s.t." => Val{:constraints},
    "bounds" => Val{:bounds},
    "bound" => Val{:bounds},
    "gen" => Val{:integer},
    "general" => Val{:integer},
    "generals" => Val{:integer},
    "bin" => Val{:binary},
    "binary" => Val{:binary},
    "binaries" => Val{:binary},
    "end" => Val{:quit},
)

const _SENSE_ALIAS = Dict(
    "max" => MOI.MAX_SENSE,
    "maximize" => MOI.MAX_SENSE,
    "maximise" => MOI.MAX_SENSE,
    "maximum" => MOI.MAX_SENSE,
    "min" => MOI.MIN_SENSE,
    "minimize" => MOI.MIN_SENSE,
    "minimise" => MOI.MIN_SENSE,
    "minimum" => MOI.MIN_SENSE,
)

const _SUBJECT_TO_ALIAS = ["subject to", "such that", "st", "s.t."]

const _CONSTRAINT_SENSE = Dict(
    "<" => :le,
    "<=" => :le,
    "=" => :eq,
    "==" => :eq,
    ">" => :ge,
    ">=" => :ge,
)

function _verify_name(variable::String, maximum_length::Int)
    if length(variable) > maximum_length
        return false
    end
    m = match(_READ_START_REG, variable)
    if m !== nothing
        return false
    end
    m = match(NAME_REG, variable)
    if m !== nothing
        return false
    end
    return true
end

mutable struct CacheLPModel
    objective_function::MOI.ScalarAffineFunction
    linear_constraint_function::MOI.ScalarAffineFunction
    linear_constraint_set::MOI.AbstractScalarSet
    linear_constraint_open::Bool
    linear_constraint_name::String
    num_linear_constraints::Int
    variables_in_model::Dict{String,MOI.VariableIndex}
    function CacheLPModel()
        return new(
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm{Float64}[], 0.0),
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm{Float64}[], 0.0),
            MOI.EqualTo(0.0),
            false,
            "",
            0,
            Dict{String,MOI.VariableIndex}(),
        )
    end
end

_set_sense!(T, model::Model, line) = nothing
function _set_sense!(::Type{Val{:obj}}, model::Model, line)
    return MOI.set(model, MOI.ObjectiveSense(), _SENSE_ALIAS[lowercase(line)])
end

function _add_new_variable!(model::Model, data_cache::CacheLPModel, name)
    var = MOI.add_variable(model)
    MOI.set(model, MOI.VariableName(), var, name)
    data_cache.variables_in_model[name] = var
    return var
end

function _get_variable_from_name(
    model::Model,
    data_cache::CacheLPModel,
    variable_name,
)
    var_inside_model = get(data_cache.variables_in_model, variable_name, "")
    if var_inside_model != ""
        return var_inside_model
    end
    options = get_options(model)
    if !_verify_name(variable_name, options.maximum_length)
        error("Invalid variable name $variable_name")
    end
    return _add_new_variable!(model, data_cache, variable_name)
end

function _tokenize(line::AbstractString)
    return String.(split(line, " "; keepempty = false))
end

function _parse_float_from_bound(val::String)
    lower_case_val = lowercase(val)
    if lower_case_val == "-inf" || lower_case_val == "-infinity"
        return -Inf
    elseif lower_case_val == "+inf" || lower_case_val == "+infinity"
        return Inf
    else
        return parse(Float64, lower_case_val)
    end
end

function _parse_affine_terms!(
    model::Model,
    data_cache::CacheLPModel,
    tokens::Vector{String},
    section::String,
    line::AbstractString,
)
    affine_terms = MOI.ScalarAffineTerm{Float64}[]
    while length(tokens) > 0
        variable = String(pop!(tokens))
        # In the case of objective functions this can be an objective constant
        if section == "objective"
            try
                obj_constant = parse(Float64, variable)
                if length(tokens) > 0
                    _sign = pop!(tokens)
                    if _sign == "-"
                        obj_constant *= -1
                    elseif _sign == "+"
                    else
                        error(
                            "Unable to parse $section due to bad operator: $(_sign) $(line)",
                        )
                    end
                end
                data_cache.objective_function.constant += obj_constant
                continue
            catch
            end
        end
        var = _get_variable_from_name(model, data_cache, variable)
        if length(tokens) > 0
            coef_token = pop!(tokens)
        else
            coeff = 1.0
            push!(affine_terms, MOI.ScalarAffineTerm(coeff, var))
            continue
        end
        try
            if coef_token == "+"
                coeff = 1.0
                push!(affine_terms, MOI.ScalarAffineTerm(coeff, var))
                continue
            elseif coef_token == "-"
                coeff = -1.0
                push!(affine_terms, MOI.ScalarAffineTerm(coeff, var))
                continue
            end
            coeff = parse(Float64, coef_token)
        catch
            error(
                "Unable to parse $section due to bad operator: $(_sign) $(line)",
            )
        end
        if length(tokens) > 0
            _sign = pop!(tokens)
            if _sign == "-"
                coeff *= -1
            elseif _sign == "+"
            else
                error(
                    "Unable to parse $section due to bad operator: $(_sign) $(line)",
                )
            end
        end
        push!(affine_terms, MOI.ScalarAffineTerm(coeff, var))
    end
    return affine_terms
end

function _parse_sos!(
    model::Model,
    data_cache::CacheLPModel,
    line::AbstractString,
)
    tokens = _tokenize(line)
    if length(tokens) < 3
        error(string("Malformed SOS constraint: ", line))
    end
    sos_con_name = String.(split(tokens[1], ":"))[1]
    if tokens[2] == "S1::"
        order = 1
    elseif tokens[2] == "S2::"
        order = 2
    else
        error("SOS of type $(tokens[2]) not recognised")
    end
    variables = MOI.VariableIndex[]
    weights = Float64[]
    for token in tokens[3:end]
        items = String.(split(token, ":"))
        if length(items) != 2
            error(string("Invalid sequence: ", token))
        end
        push!(variables, _get_variable_from_name(model, data_cache, items[1]))
        push!(weights, parse(Float64, items[2]))
    end
    sos_con = MOI.add_constraint(
        model,
        variables,
        order == 1 ? MOI.SOS1(weights) : MOI.SOS2(weights),
    )
    MOI.set(model, MOI.ConstraintName(), sos_con, sos_con_name)
    # TODO I think this only works for SOS of one line
    return
end

function _parse_variable_type!(
    model::Model,
    data_cache::CacheLPModel,
    line::AbstractString,
    set::MOI.AbstractSet,
)
    items = _tokenize(line)
    for v in items
        var = _get_variable_from_name(model, data_cache, v)
        MOI.add_constraint(model, var, set)
    end
    return nothing
end

function _parse_section!(
    ::Type{Val{:none}},
    model::Model,
    data_cache::CacheLPModel,
    line::AbstractString,
)
    return nothing
end
function _parse_section!(
    ::Type{Val{:quit}},
    model::Model,
    data_cache::CacheLPModel,
    line::AbstractString,
)
    return error("Corrupted LP File. You have the lne $(line) after an end.")
end
function _parse_section!(::Type{Val{:integer}}, model, data_cache, line)
    return _parse_variable_type!(model, data_cache, line, MOI.Integer())
end
function _parse_section!(::Type{Val{:binary}}, model, data_cache, line)
    return _parse_variable_type!(model, data_cache, line, MOI.ZeroOne())
end

function _parse_section!(
    ::Type{Val{:obj}},
    model::Model,
    data_cache::CacheLPModel,
    line::AbstractString,
)
    # okay so line should be the start of the objective
    if occursin(":", line)
        # throw away name
        m = match(r"(.*?)\:(.*)", line)
        line = String(m[2])
    end
    tokens = _tokenize(line)
    if length(tokens) == 0 # no objective
        return MOI.set(model, MOI.ObjectiveSense(), MOI.FEASIBILITY_SENSE)
    end
    affine_terms =
        _parse_affine_terms!(model, data_cache, tokens, "objective", line)
    push!(data_cache.objective_function.terms, affine_terms...)
    return
end

function _parse_section!(
    ::Type{Val{:constraints}},
    model::Model,
    data_cache::CacheLPModel,
    line::AbstractString,
)
    if match(r" S([0-9]):: ", line) !== nothing
        # it's an SOS constraint
        _parse_sos!(model, data_cache, line)
        return
    end
    if data_cache.linear_constraint_open == false
        # parse the number of rows and add this name
        data_cache.linear_constraint_name = "R$(data_cache.num_linear_constraints)"
    end
    if occursin(":", line)
        if data_cache.linear_constraint_open == true
            error("Malformed constraint $(line). Is the previous one valid?")
        end
        # throw away name
        m = match(r"(.*?)\:(.*)", line)
        data_cache.linear_constraint_name = String(m[1])
        line = String(m[2])
    end
    data_cache.linear_constraint_open = true

    tokens = _tokenize(line)
    if length(tokens) == 0 # no entries
        return
    elseif length(tokens) >= 2 && haskey(_CONSTRAINT_SENSE, tokens[end-1])# test if constraint ends this line
        rhs = parse(Float64, pop!(tokens))
        sym = pop!(tokens)
        if _CONSTRAINT_SENSE[sym] == :le
            data_cache.linear_constraint_set = MOI.LessThan(rhs)
        elseif _CONSTRAINT_SENSE[sym] == :ge
            data_cache.linear_constraint_set = MOI.GreaterThan(rhs)
        elseif _CONSTRAINT_SENSE[sym] == :eq
            data_cache.linear_constraint_set = MOI.EqualTo(rhs)
        end
        # Finished
        # Add constraint
        c = MOI.add_constraint(
            model,
            data_cache.linear_constraint_function,
            data_cache.linear_constraint_set,
        )
        MOI.set(
            model,
            MOI.ConstraintName(),
            c,
            data_cache.linear_constraint_name,
        )
        data_cache.num_linear_constraints += 1
        # Clear the constraint part of data_cache
        data_cache.linear_constraint_set = MOI.EqualTo(0.0)
        data_cache.linear_constraint_function =
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm{Float64}[], 0.0)
        data_cache.linear_constraint_name = ""
        data_cache.linear_constraint_open = false
    end
    affine_terms =
        _parse_affine_terms!(model, data_cache, tokens, "constraint", line)
    push!(data_cache.linear_constraint_function.terms, affine_terms...)
    return
end

_bound_error(line::AbstractString) = error("Unable to parse bound: $(line)")
function _parse_section!(
    ::Type{Val{:bounds}},
    model::Model,
    data_cache::CacheLPModel,
    line::AbstractString,
)
    items = _tokenize(line)
    v = ""
    lb = -Inf
    ub = Inf
    if length(items) == 5 # ranged bound
        v = items[3]
        if (items[2] == "<=" || items[2] == "<") &&
           (items[4] == "<=" || items[4] == "<") # le
            lb = _parse_float_from_bound(items[1])
            ub = _parse_float_from_bound(items[5])
        elseif (items[2] == ">=" || items[2] == ">") &&
               (items[4] == ">=" || items[4] == ">") # ge
            lb = _parse_float_from_bound(items[5])
            ub = _parse_float_from_bound(items[1])
        else
            _bound_error(line)
        end
    elseif length(items) == 3 # one sided
        v = items[1]
        if items[2] == "<=" || items[2] == "<" # le
            ub = _parse_float_from_bound(items[3])
            if ub > 0.0
                lb = 0.0
            else
                lb = -Inf
            end
        elseif items[2] == ">=" || items[2] == ">" # ge
            lb = _parse_float_from_bound(items[3])
            ub = +Inf
        elseif items[2] == "==" || items[2] == "=" # eq
            lb = ub = _parse_float_from_bound(items[3])
        else
            _bound_error(line)
        end
    elseif length(items) == 2 # free
        if items[2] != "free"
            _bound_error(line)
        end
        v = items[1]
    else
        _bound_error(line)
    end
    var = _get_variable_from_name(model, data_cache, v)
    set = bounds_to_set(lb, ub)
    if set !== nothing
        MOI.add_constraint(model, var, set)
    end
    return
end

function bounds_to_set(lower::Float64, upper::Float64)
    if -Inf < lower < upper < Inf
        return MOI.Interval(lower, upper)
    elseif -Inf < lower && upper == Inf
        return MOI.GreaterThan(lower)
    elseif -Inf == lower && upper < Inf
        return MOI.LessThan(upper)
    elseif lower == upper
        return MOI.EqualTo(upper)
    end
    return  # free variable
end

function _add_objective!(model::Model, data_cache::CacheLPModel)
    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        data_cache.objective_function,
    )
    return
end

"""
    Base.read!(io::IO, model::FileFormats.LP.Model)

Read `io` in the LP file format and store the result in `model`.
"""
function Base.read!(io::IO, model::Model)
    if !MOI.is_empty(model)
        error("Cannot read in file because model is not empty.")
    end
    data_cache = CacheLPModel()
    section = Val{:none}
    while !eof(io)
        line = string(strip(readline(io)))
        line = _strip_comment(line)
        if line == "" # skip blank lines
            continue
        end
        if haskey(_KEYWORDS, lowercase(line)) # section has changed
            section = _KEYWORDS[lowercase(line)]
            _set_sense!(section, model, line)
            continue
        end
        _parse_section!(section, model, data_cache, line)
    end
    _add_objective!(model, data_cache)
    return
end

end


text:
#
# This file is part of the DiscreteEvents.jl Julia package, MIT license
#
# Paul Bayer, 2020
#
# This is a Julia package for discrete event simulation
#

abstract type AbstractResource end

"""
    Resource{T}(capacity)

A Resource implements a [`Deque`](https://juliacollections.github.io/DataStructures.jl/latest/deque/)
with a limited capacity. If used in multithreading applications, the user must avoid race
conditions by explicitly wrapping modifying calls with `lock-unlock`.

# Fields
- `items::Deque{T}`: resource buffer
- `capacity::Number=Inf`: the capacity is limited to the given integer,
- `lock::ReentrantLock`: a lock for coordinating resource access by tasks.

# Example
```jldoctest
julia>
```
!!! note
    In order to use the full interface to `Resource` you have to load `DataStructures`.
"""
mutable struct Resource{T} <: AbstractResource
    items::Deque{T}
    capacity::Int
    lock::ReentrantLock

    function Resource{T}(capacity::U=Inf) where {T,U<:Number}
        if isinf(capacity)
            new{T}(Deque{T}(), typemax(Int), ReentrantLock())
        else
            new{T}(Deque{T}(), capacity, ReentrantLock())
        end
    end
end

"""
    capacity(r::Resource)

Get the capacity of a resource
"""
DataStructures.capacity(r::Resource) = r.capacity

"""
    isfull(r::Resource)

Test whether the resource is full
"""
DataStructures.isfull(r::Resource) = length(r) ≥ capacity(r)

"""
    isready(r::Resource)

Test whether an item is available.
"""
Base.isready(r::Resource) = !isempty(r.items)

"""
    isempty(r::Resource)

Test whether the resource is empty.
"""
Base.isempty(r::Resource) = isempty(r.items)

"""
    empty!(r::Resource)

Reset the resource buffer (deque).
"""
Base.empty!(r::Resource) = empty!(r.items)

"""
    length(r::Resource)

Get the number of elements available.
"""
Base.length(r::Resource) = length(r.items)

"""
    push!(r::Resource, x)

Add an element to the back of a resource deque.
"""
Base.push!(r::Resource, x) =
    !isfull(r) ? push!(r.items, x) : throw(ArgumentError("Resource must be non-full"))

"""
    pop!(r::Resource)

Remove an element from the back of a resource deque.
"""
Base.pop!(r::Resource) = pop!(r.items)

"""
    pushfirst!(r::Resource, x)

Add an element to the front of a resource deque.
"""
Base.pushfirst!(r::Resource, x) =
    !isfull(r) ? pushfirst!(r.items, x) : throw(ArgumentError("Resource must be non-full"))

"""
    popfirst!(r::Resource)

Remove an element from the front of a resource deque.
"""
Base.popfirst!(r::Resource) = popfirst!(r.items)

"""
    first(r::Resource)

Get the element at the front of a resource deque.
"""
Base.first(r::Resource) = first(r.items)

"""
    last(r::Resource)

Get the element at the back of a resource deque.
"""
Base.last(r::Resource) = last(r.items)

"""
    lock(r::Resource)

Acquire the resource lock when it becomes available. If the lock is already locked by a
different task/thread, wait for it to become available.

Each lock must be matched by an unlock.
"""
Base.lock(r::Resource) = lock(r.lock)

"""
    unlock(r::Resource)

Releases ownership of the resource lock.

If this is a recursive lock which has been acquired before, decrement an internal counter
and return immediately.
"""
Base.unlock(r::Resource) = unlock(r.lock)

"""
    islocked(r::Resource)

Check whether the lock is held by any task/thread. This should not be used for
synchronization (see instead trylock).
"""
Base.islocked(r::Resource) = islocked(r.lock)

"""
    trylock(r::Resource)

Acquire the resource lock if it is available, and return true if successful. If the lock
is already locked by a different task/thread, return false.

Each successful trylock must be matched by an unlock.
"""
Base.trylock(r::Resource) = trylock(r.lock)

# -----------------------------------------------------
# extended interface to Julia channels
# -----------------------------------------------------

"""
    capacity(ch::Channel)

Get the capacity of a channel.
"""
DataStructures.capacity(ch::Channel) = ch.sz_max

"""
    length(ch::Channel)

Get the number of items in a channel.
"""
Base.length(ch::Channel) = length(ch.data)

"""
    isfull(ch::Channel)

Test whether a channel is full.
"""
DataStructures.isfull(ch::Channel) = length(ch) ≥ capacity(ch)

"""
    isempty(ch::Channel)

Test whether a channel is empty.
"""
Base.isempty(ch::Channel) = isempty(ch.data)

"""
    first(ch::Channel)

Get the first item in a channel without removing it. Throws an error if the channel is empty.
If you want a blocking operation, use [`fetch`](https://docs.julialang.org/en/v1/base/parallel/#Base.fetch-Tuple{Channel}).
"""
Base.first(ch::Channel) = first(ch.data)

"""
    empty!(ch::Channel)

Reset a channel, throw away the elements stored to it.
"""
Base.empty!(ch::Channel) = foreach((x)->take!(ch), 1:length(ch))