text:
using Revise
using MathOptInterface
const MOI = MathOptInterface
using JuMP
using DataStructures

### ============================================================================
### Objective expression
### ============================================================================
mutable struct _Objective
    # constant value
    constant::Float64
    # terms
    terms::Vector{MOI.ScalarAffineTerm{Float64}}

    function _Objective(constant::Float64, terms::Vector{MOI.ScalarAffineTerm{Float64}})
        return new(constant, terms)
    end
end

### ============================================================================
### Variables
### ============================================================================
mutable struct _VariableInfo
    # Index of the variable
    index::MOI.VariableIndex
    # The block that the variable appears in.
    block::Int64
    # The quantifier of the variable.
    quantifier::String

    function _VariableInfo(index::MOI.VariableIndex, block::Int64, quantifier::String)
        return new(index, block, quantifier)
    end
end

### ============================================================================
### Variable constraints
### ============================================================================
struct _VariableConstraintInfo
    # Constraint index
    index::MOI.ConstraintIndex
    # Variable Index
    vindex::MOI.VariableIndex
    # Constraint set
    conSet::Any

    function _VariableConstraintInfo(ind::MOI.ConstraintIndex, vind::MOI.VariableIndex, set)
        return new(ind, vind, set)
    end
end

### ============================================================================
### Constraints
### ============================================================================
struct _ConstraintInfo
    # Constraint index
    index::MOI.ConstraintIndex
    # Constraint ScalarAffineFunction
    scalarAffineFunction::MOI.ScalarAffineFunction{Float64}
    # Constraint set
    conSet::Any
    # Constraint quantifier
    quantifier::String

    function _ConstraintInfo(ind::MOI.ConstraintIndex, scalarAffFunc::MOI.ScalarAffineFunction{Float64}, set, quantifier)
        return new(ind, scalarAffFunc, set, quantifier)
    end
end

### ============================================================================
### Results
### ============================================================================
struct _Results
    # status
    #raw_status_string::String
    #termination_status::MOI.TerminationStatusCode

    objective_value::Float64
    runtime::Float64
    #decisionNodes::Int64
    #propagationSteps::Int64
    #learntConstraints::Int64

    # quality
    solutionStatus::String
    gap::Float64

    # variable values
    values::Dict{String, Float64}

    function _Results(obj::Float64, runtime::Float64, solStatus::String, gap::Float64, values::Dict{String, Float64})
        return new(obj, runtime, solStatus, gap, values)
    end
end

"""
    AbstractSolverCommand
An abstract type that allows over-riding the call behavior of the solver.
"""
abstract type AbstractSolverCommand end

"""
    call_solver(
        solver::AbstractSolverCommand,
        qip_filename::String,
        options::Vector{String},
        stdin::IO,
        stdout::IO,
    )::String
Execute the `solver` given the QIP file at `qip_filename`, a vector of `options`,
and `stdin` and `stdout`. Return the filename of the resulting `.sol` file.
"""
function call_solver end

struct _DefaultSolverCommand{F} <: AbstractSolverCommand
    f::F
end

function call_solver(
    solver::_DefaultSolverCommand,
    solverPath::String,
    qip_filename::String,
    options::Vector{Int64},
    stdin::IO,
    stdout::IO,
    output::String,
)
    # parse optimizer attributes, value for time_limit is -1 if not supposed to be set
    cmd =  ``
    if options[2] == -1
        cmd = `$(solverPath) $(qip_filename) $(options[1])`
    else
        cmd = `$(solverPath) $(qip_filename) $(options[1]) $(options[2])`
    end

    solver.f() do solver_path
        ret = run(
            pipeline(
                cmd,
                stdin = stdin,
                stdout = output,
                append=true,
                stderr= output,
            ),
        )
        if ret.exitcode != 0
            error("Nonzero exit code: $(ret.exitcode)")
        end
    end
end

_solver_command(x::String) = _DefaultSolverCommand(f -> f(x))
_solver_command(x::Function) = _DefaultSolverCommand(x)
_solver_command(x::AbstractSolverCommand) = x

### ============================================================================
### Optimizer
### ============================================================================
mutable struct Optimizer <: MOI.AbstractOptimizer
    solver_command::AbstractSolverCommand
    stdin::Any
    stdout::Any
    # result information
    results::_Results
    # solve time
    solve_time::Float64
    # Store MOI.Name().
    name::String
    # The objective expression.
    o::_Objective
    sense::MOI.OptimizationSense
    #is_objective_set::Bool
    # A vector of variable constraints
    vc::Vector{_VariableConstraintInfo}
    # A vector of constraints
    c::Vector{_ConstraintInfo}
    cNumber::Int64
    # A dictionary of info for the variables.
    v::Dict{MOI.VariableIndex,_VariableInfo}
    # was optimizer called
    optimize_not_called::Bool
    # termination status
    termination_status::MOI.TerminationStatusCode

    # solver specific attributes
    # time limit in seconds
    time_limit::Int64
    # output information level
    output_info::Int64
    # problem file name
    problem_file::String
    # name of solver .exe
    solver_path::String
end

"""
    Optimizer(
        solver_command::Union{String,Function},
        stdin::Any = stdin,
        stdout:Any = stdout,
    )
Create a new Optimizer object.

`solver_command`:
* A `String` of the full path of an Yasol executable.

Redirect IO using `stdin` and `stdout`. These arguments are passed to
`Base.pipeline`. See the Julia documentation for more details.
"""
function Optimizer(
        solver_command::Union{AbstractSolverCommand,String,Function} = "",
        stdin::Any = stdin,
        stdout::Any = stdout,
)
    return Optimizer(
        _solver_command(solver_command),
        stdin,
        stdout,
        _Results(
            0.0,
            0.0,
            "",
            0.0,
            Dict{String,Float64}()
        ),
        NaN,
        "",
        _Objective(0.0, MOI.ScalarAffineTerm{Float64}[]),
        MOI.FEASIBILITY_SENSE,
        _VariableConstraintInfo[],
        _ConstraintInfo[],
        0,
        Dict{MOI.VariableIndex,_VariableInfo}(),
        true,
        MOI.OPTIMIZE_NOT_CALLED,
        -1,
        -1,
        "",
        ""
    )
end

Base.show(io::IO, ::Optimizer) = print(io, "A YASOL model")

MOI.get(model::Optimizer, ::MOI.SolverName) = "YASOL"
MOI.get(model::Optimizer, ::MOI.RawSolver) = model

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

MOI.get(model::Optimizer, ::MOI.NumberOfVariables) = length(model.v)

function MOI.empty!(model::Optimizer)
    #model.results = _QIPResults()
    model.solve_time = NaN
    model.o = _Objective(0.0, MOI.ScalarAffineTerm{Float64}[])
    model.sense = MOI.FEASIBILITY_SENSE
    model.vc = _VariableConstraintInfo[]
    model.c = _ConstraintInfo[]
    model.cNumber = 0
    model.v = Dict{MOI.VariableIndex,_VariableInfo}()
    model.time_limit = -1
    model.output_info = -1
    model.problem_file = ""
    model.optimize_not_called = true
    return
end

function MOI.is_empty(model::Optimizer)
    if isempty(model.vc) && isempty(model.c) && isempty(model.v)
        return true
    else
        return false
    end
end

# ========================================
#   Supported constraints and objectives
# ========================================
const _SCALAR_FUNCTIONS = Union{
    MOI.VariableIndex,
    MOI.ScalarAffineFunction{Float64},
}

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

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

MOI.supports(::Optimizer, ::MOI.ObjectiveSense) = true
MOI.supports(::Optimizer, ::MOI.ObjectiveFunction{<:MOI.ScalarAffineFunction}) = true

# ========================================
#   Copy_to functionality. No incremental modification supported.
# ========================================
#MOI.Utilities.supports_default_copy_to(::Optimizer, ::Bool) = false -> deprecated!!!
MOI.supports_incremental_interface(::Optimizer) = false

function MOI.copy_to(
    dest::Optimizer,
    model::MOI.ModelLike;
    copy_names::Bool = false,
)

    mapping = MOI.Utilities.IndexMap()

    # copy optimizer attributes
    try
        dest.time_limit = MOI.get(model, MOI.RawOptimizerAttribute("time limit"))
    catch
        dest.time_limit = -1
    end

    try
        dest.output_info = MOI.get(model, MOI.RawOptimizerAttribute("output info"))
    catch
        dest.output_info = 1
    end
    try
        dest.problem_file = MOI.get(model, MOI.RawOptimizerAttribute("problem file name"))
    catch
        @error string("Please provide a problem file name! No problem file could be written!")
        return
    end

    # copy objective sense
    dest.sense = MOI.get(model, MOI.ObjectiveSense())

    # copy variables

    # save the quantifier and block of the last variable
    last_block = 0
    last_quantifiers = []
    for v in MOI.get(model, MOI.ListOfVariableIndices())

        try
            quantifier = MOI.get(model, YasolSolver.VariableAttribute("quantifier"), v)
            block = MOI.get(model, YasolSolver.VariableAttribute("block"), v)

            if quantifier !== nothing && block !== nothing
                dest.v[v] = _VariableInfo(v, block, quantifier)

                if block > last_block
                    last_block = block
                end
            else
                @error string("You need to set a quantifier and a block for each variable when using Yasol solver!")
                return
            end

        catch err
            #println("You need to set a quantifier and a block for each variable when using Yasol solver!")
            #@warn string("You need to set a quantifier and a block for each variable when using Yasol solver!")
            @warn string(err)
        end

        mapping[v] = v
    end

    # show warning, if variable in last block is not existential
    for v in MOI.get(model, MOI.ListOfVariableIndices())
        try
            quantifier = MOI.get(model, YasolSolver.VariableAttribute("quantifier"), v)
            block = MOI.get(model, YasolSolver.VariableAttribute("block"), v)

            if block == last_block
                push!(last_quantifiers, quantifier)
            end
        catch err
            @warn string(err)
        end
    end

    for q in last_quantifiers
        if q != "exists"
            @error string("The variable in the last block needs to be existential! Please add a dummy variable!")
            return
        end
    end

    # copy objective function
    F = MOI.get(model, MOI.ObjectiveFunctionType())
    obj = MOI.get(model, MOI.ObjectiveFunction{F}())
    temp = _Objective(obj.constant, obj.terms)
    dest.o = temp

    # copy constraints
    for (F, S) in MOI.get(model, MOI.ListOfConstraintTypesPresent())
        for ci in MOI.get(model, MOI.ListOfConstraintIndices{F,S}())
            mapping[ci] = ci

            f = MOI.get(model, MOI.ConstraintFunction(), ci)
            s = MOI.get(model, MOI.ConstraintSet(), ci)

            # get constraint quantifier
            q = ""
            q_temp = MOI.get(model, YasolSolver.ConstraintAttribute("quantifier"), ci)
            if q_temp !== nothing
                q = q_temp
            else
                q = ""
            end

            if typeof(f) == MathOptInterface.VariableIndex
                vcon = _VariableConstraintInfo(ci, f, s)
                push!(dest.vc, vcon)
            else
                con = _ConstraintInfo(ci, f, s, q)
                push!(dest.c, con)
            end
        end
    end

    # count constraints
    values = []
    for con in dest.c
        push!(values, Int64(con.index.value))
    end
    dest.cNumber = length(values)

    return mapping
end

# ========================================
#   Write model to file.
# ========================================
function Base.write(io::IO, qipmodel::Optimizer)

    # print objective sense
    if(qipmodel.sense == MOI.MIN_SENSE)
        println(io, "MINIMIZE")
    elseif qipmodel.sense == MOI.MAX_SENSE
        println(io, "MAXIMIZE")
    end

    # number of variables
    numVar = length(qipmodel.v)
    exist = []
    all = []
    binaries = []
    generals = []

    # print objective function
    for term in qipmodel.o.terms
        # print coefficient and variables
        if term.coefficient < 0.0
            print(io, string(term.coefficient) * "x" * string(term.variable.value) * " ")
        elseif term.coefficient > 0.0
            print(io, "+" * string(term.coefficient) * "x" * string(term.variable.value) * " ")
        end
    end

    # print constant value if != 0
    if qipmodel.o.constant < 0.0
        print(io, string(qipmodel.o.constant))
    else qipmodel.o.constant > 0.0
        print(io, "+ " * string(qipmodel.o.constant))
    end

    println(io, "")

    # print constraints
    println(io, "SUBJECT TO")
    for con in qipmodel.c
        # print quantifier if set
        if con.quantifier !== ""
            if con.quantifier === "exists"
                print(io, "E_Constraint" * string(con.index.value) * ": ")
            elseif con.quantifier === "all"
                print(io, "U_Constraint" * string(con.index.value) * ": ")
            end
        end

        # print terms
        for term in con.scalarAffineFunction.terms
            if term.coefficient < 0.0
                print(io, string(term.coefficient) * "x" * string(term.variable.value) * " ")
            elseif term.coefficient > 0.0
                print(io, "+" * string(term.coefficient) * "x" * string(term.variable.value) * " ")
            end
        end

        temp = 0.0
        temp = temp + (con.scalarAffineFunction.constant)*-1

        """
        # print constant
        if con.scalarAffineFunction.constant < 0.0
            #print(io, string(con.scalarAffineFunction.constant))
            temp+= (con.scalarAffineFunction.constant)*-1
        else con.scalarAffineFunction.constant > 0.0
            #print(io, "+ " * string(con.scalarAffineFunction.constant))
            temp+= (con.scalarAffineFunction.constant)*-1
        end
        """

        # print set
        if typeof(con.conSet) == MathOptInterface.GreaterThan{Float64}
            if temp !== 0.0
                print(io, ">= " * string((con.conSet.lower) + temp))
            else
                print(io, ">= " * string((con.conSet.lower)))
            end
        elseif typeof(con.conSet) == MathOptInterface.LessThan{Float64}
            if temp !== 0.0
                print(io, "<= " * string((con.conSet.upper) + temp))
            else
                print(io, "<= " * string((con.conSet.upper)))
            end
        end
        println(io, "")
    end

    # print variable bounds
    println(io, "BOUNDS")
    for i in 1:numVar
        lower = -9999.9
        upper = 9999.9
        type = nothing
        for varCon in qipmodel.vc
            if varCon.vindex.value == i
                if typeof(varCon.conSet) == MathOptInterface.Integer
                    push!(generals, "x"*string(i))
                    type = "int"
                elseif typeof(varCon.conSet) == MathOptInterface.ZeroOne
                    push!(binaries, "x"*string(i))
                    type = "binary"
                elseif typeof(varCon.conSet) == MathOptInterface.GreaterThan{Float64}
                    lower = varCon.conSet.lower
                elseif typeof(varCon.conSet) == MathOptInterface.LessThan{Float64}
                    upper = varCon.conSet.upper
                end
            end
        end

        # show warning, if variable has no lower or upper bound
        if (lower == -9999.9 || upper == 9999.9)
            @error string("Every variable needs to be bounded from above and below (binary variables as well)!")
            return
        end

        # write bounds
        println(io, string(lower) * " <= " * "x" * string(i) * " <= " * string(upper))
    end

    # check, if all variables are integer or binary
    for a in all
        if !(a in binaries) && (!a in generals)
            @error string("All variables need to be binary or integer!")
            return
        end
    end

    # print binaries
    println(io, "BINARIES")
    bin = ""
    for b in binaries
        bin *= b * " "
    end
    println(io, bin)

    # print generals
    println(io, "GENERALS")
    gen = ""
    for g in generals
        gen *= g * " "
    end
    println(io, gen)

    # print exists
    println(io, "EXISTS")
    exists = ""
    for i in 1:MOI.get(qipmodel, MOI.NumberOfVariables())
        if qipmodel.v[MOI.VariableIndex(i)].quantifier === "exists"
            exists = exists * "x" * string(i) * " "
        end
    end
    println(io, exists)

    # print all
    println(io, "ALL")
    all = ""
    for i in 1:MOI.get(qipmodel, MOI.NumberOfVariables())
        if qipmodel.v[MOI.VariableIndex(i)].quantifier === "all"
            all = all * "x" * string(i) * " "
        end
    end
    println(io, all)

    # print order
    println(io, "ORDER")
    order = ""

    tempDict = OrderedDict{Int64,Int64}()
    for i in 1:MOI.get(qipmodel, MOI.NumberOfVariables())
        tempDict[i] = qipmodel.v[MathOptInterface.VariableIndex(i)].block
    end

    # sort temp dict
    last_block = 0
    last_block_variables = []
    sorted = sort!(tempDict, byvalue=true)

    # build order string
    for (key, value) in sorted
        order = order * "x" * string(key) * " "
        if value > last_block
            last_block = value
        end
    end

    # save all last block variables
    for (key, value) in sorted
        if value == last_block
            push!(last_block_variables, key)
        end
    end

    # make sure, that continuos variables are only allowed in the last block
    for (key, value) in sorted
        if value != last_block
            if !("x"*string(key) in generals) && !("x"*string(key) in binaries)
                @error string("Continuos variables are only allowed in the last block")
                return
            end
        end
    end

    #@warn string(last_block)
    #@warn string(last_block_variables)

    println(io, order)

    println(io, "END")
end

# ========================================
#   Call yasol solver with problem file and parameters.
# ========================================
function MOI.optimize!(model::Optimizer)

    # check if problem file name is set, show warning otherwise
    if(model.problem_file == "")
        @error string("Please provide a problem file name! No problem file could be written!")
        return
    end

    # check if Yasol.ini file is given in the solver folder
    path = joinpath(pwd(), "Yasol.ini")
    if !isfile(String(path))
        @warn string("No Yasol.ini file was found in the solver folder!")
    end

    # check if Yasol .exe is available under given path
    if !isfile(String(model.solver_path*".exe")) && !isfile(String(model.solver_path))
        @error string("No Yasol executable was found under the given path!")
        return
    end

    model.optimize_not_called = false

    options = [model.output_info, model.time_limit]

    start_time = time()
    # create problem file
    qlp_file = joinpath(pwd(), model.problem_file)

    output_file = joinpath(pwd(), (rsplit(model.problem_file, ".")[1]) * "_output.txt")
    touch(output_file)
    open(io -> write(io, model), qlp_file, "w")

    # call solver
    try
        call_solver(
            model.solver_command,
            model.solver_path,
            qlp_file,
            options,
            model.stdin,
            model.stdout,
            output_file,
        )
        # read solution & set results
        model.results = YasolSolver.importSolution(qlp_file * ".sol")

        if model.results.solutionStatus == "OPTIMAL"
            model.termination_status = MOI.OPTIMAL
        end

    catch err
        # TODO show error in results
    end

    model.solve_time = time() - start_time
    return
end

function MOI.write_to_file(model::Optimizer, filename::String)
    open(io -> write(io, model), filename, "w")
    return
end

# ========================================
#   Model attributes
# ========================================
MOI.supports(::Optimizer, ::MOI.RawOptimizerAttribute) = true

function MOI.set(model::Optimizer, param::MOI.RawOptimizerAttribute, value)
    if param == MOI.RawOptimizerAttribute("output info")
        model.output_info = Int64(value)
    elseif param == MOI.RawOptimizerAttribute("time limit")
        model.time_limit = Int64(value)
    elseif param == MOI.RawOptimizerAttribute("problem file name")
        model.problem_file = String(value)
        # check if problem file already exists
        if isfile(String(value))
            @warn "A file with the chosen name already exists. You are about to overwrite that file."
        end
        # check if solution file already exists
        if isfile(String(value) * ".sol")
            @warn "A solution file for the problem already exists. If you create another solution with the same name, you cannot import the new solution using JuMP."
        end
    elseif param == MOI.RawOptimizerAttribute("solver path")
        model.solver_path = String(value)
    end
    return
end

function MOI.get(model::Optimizer, param::MOI.RawOptimizerAttribute)
    if param == MOI.RawOptimizerAttribute("output info")
        return model.output_info
    elseif param == MOI.RawOptimizerAttribute("time limit")
        return model.time_limit
    elseif param == MOI.RawOptimizerAttribute("problem file name")
        return model.problem_file
    elseif param == MOI.RawOptimizerAttribute("solver path")
        return model.solver_path
    end
end

# ========================================
#   Variable attributes
# ========================================
struct VariableAttribute <: MOI.AbstractVariableAttribute
    name::String
end

function MOI.supports(
    model::Optimizer,
    attr::VariableAttribute,
    ::Type{MOI.VariableIndex},
)
    if attr.name === "quantifier" || attr.name === "block"
        return true
    else
        return false
    end
end

function MOI.get(
    model::Optimizer,
    attr::VariableAttribute,
    vi::MOI.VariableIndex,
)
    if attr === "quantifier"
        return model.v[vi].quantifier
    elseif attr === "block"
        return model.v[vi].block
    end
end

# variable attribute 'quantifier'
function MOI.set(
    model::Optimizer,
    attr::VariableAttribute,
    vi::MOI.VariableIndex,
    value::String,
)
    if attr === "quantifier" && value === "all"
        model.v[vi].quantifier = "all"
    elseif attr === "quantifier" && value === "exists"
        model.v[vi].quantifier = "exists"
    end
    return
end

# variable attribute 'block'
function MOI.set(
    model::Optimizer,
    attr::VariableAttribute,
    vi::MOI.VariableIndex,
    value::Int,
)
    if attr === "block"
        model.v[vi].block = Int64(value)
    end
    return
end

# ========================================
#   Constraint attributes
# ========================================
struct ConstraintAttribute <: MOI.AbstractConstraintAttribute
    name::String
end

function MOI.supports(
    model::Optimizer,
    attr::ConstraintAttribute,
    ::Type{<:MOI.ConstraintIndex},
)
    if attr.name === "quantifier"
        return true
    else
        return false
    end
end

function MOI.get(
    model::Optimizer,
    attr::ConstraintAttribute,
    ci::MOI.ConstraintIndex,
)
    if attr === "quantifier"
        for con in model.c
            if con.index == ci
                return con.quantifier
            end
        end
    end
end

function MOI.set(
    model::Optimizer,
    attr::ConstraintAttribute,
    ci::MOI.ConstraintIndex,
    value::String,
)
    if attr === "quantifier"
        for con in model.c
            if con.index == ci
                con.quantifier = value
            end
        end
    end
    return
end

# ========================================
#   Solution & TerminationStatus
# ========================================
MOI.supports(::Optimizer, ::MOI.TerminationStatus) = true
MOI.supports(::Optimizer, ::MOI.TimeLimitSec) = true

# return termination status
function MOI.get(model::Optimizer, attr::MOI.TerminationStatus)
    try
        if model.optimize_not_called
            return MOI.OPTIMIZE_NOT_CALLED
        elseif model.results.solutionStatus == "OPTIMAL"
            return MOI.OPTIMAL
        elseif model.results.solutionStatus == "INFEASIBLE"
            return MOI.INFEASIBLE
        elseif model.results.solutionStatus == "INCUMBENT"
            return MOI.TIME_LIMIT
        else
            return MOI.OTHER_ERROR
        end
    catch err
        println(err)
    end
end

# return solve time
function MOI.get(model::Optimizer, attr::MOI.SolveTimeSec)
    try
        return Float64(model.results.runtime)
    catch err
        println(err)
    end
end

# return objective value
function MOI.get(model::Optimizer, attr::MOI.ObjectiveValue)
    try
        return model.results.objective_value
    catch err
        println(err)
    end
end

# return variable value
function MOI.get(
    model::Optimizer,
    attr::MOI.VariablePrimal,
    x::MOI.VariableIndex,
)

    try
        return model.results.values["x" * string(x.value)]
    catch err
        println(err)
    end
end


text:
using Test
using MPI

MPI.Init()

function allgather(A)
    comm = MPI.COMM_WORLD
    C = MPI.Allgather(A, comm)
end

comm = MPI.COMM_WORLD

for typ in Base.uniontypes(MPI.MPIDatatype)
    A = typ[MPI.Comm_rank(comm) + 1]
    C = allgather(A)
    @test typeof(C) === Vector{typ}
    @test C == map(typ, collect(1:MPI.Comm_size(comm)))
    A = convert(typ,MPI.Comm_rank(comm) + 1)
    C = allgather(A)
    @test typeof(C) === Vector{typ}
    @test C == map(typ, collect(1:MPI.Comm_size(comm)))

    val = typ(MPI.Comm_rank(comm) + 1)

	# Test passing output buffer with set size
	A = typ[val]
	C = fill(val, MPI.Comm_size(comm))
	MPI.Allgather!(A, C, length(A), comm)
	@test C == map(typ, collect(1:MPI.Comm_size(comm)))

	# Test assertion error
	C = fill(val, MPI.Comm_size(comm)-1)
	@test_throws AssertionError MPI.Allgather!(A, C, length(A), comm)

	# Test explicit IN_PLACE
	A = typ[val]
	C = fill(val, MPI.Comm_size(comm))
	C[MPI.Comm_rank(comm)+1] = MPI.Comm_rank(comm) + 1
	MPI.Allgather!(MPI.IN_PLACE, C, 1, comm)
	@test typeof(C) === Vector{typ}
	@test C == map(typ, collect(1:MPI.Comm_size(comm)))

	# Test IN_PLACE
	A = typ[val]
	C = fill(val, MPI.Comm_size(comm))
	C[MPI.Comm_rank(comm)+1] = MPI.Comm_rank(comm) + 1
	MPI.Allgather!(C, 1, comm)
	@test typeof(C) === Vector{typ}
	@test C == map(typ, collect(1:MPI.Comm_size(comm)))

end

MPI.Finalize()
@test MPI.Finalized()

ext: jl


text:
"""
    Function from GridapPETSC to find PETSc lib location.
"""
function get_petsc_location()
    PETSC_DIR   = haskey(ENV,"PETSC_DIR") ? ENV["PETSC_DIR"] : "/usr/lib/petsc"
    PETSC_ARCH  = haskey(ENV,"PETSC_ARCH") ? ENV["PETSC_ARCH"] : ""
    PETSC_LIB = ""

    # Check PETSC_DIR exists
    if isdir(PETSC_DIR)

        # Define default paths
        PETSC_LIB_DIR = joinpath(PETSC_DIR,PETSC_ARCH,"lib")

        # Check PETSC_LIB (.../libpetsc.so or .../libpetsc_real.so file) exists
        if isfile(joinpath(PETSC_LIB_DIR,"libpetsc.so"))
            PETSC_LIB = joinpath(PETSC_LIB_DIR,"libpetsc.so")
        elseif isfile(joinpath(PETSC_LIB_DIR,"libpetsc_real.so"))
            PETSC_LIB = joinpath(PETSC_LIB_DIR,"libpetsc_real.so")
        end
    end

    # PETSc lib not found
    if(length(PETSC_LIB) == 0)
        # Workaround for automerging on RegistryCI or doc deployment
        if(haskey(ENV,"JULIA_REGISTRYCI_AUTOMERGE") || haskey(ENV, "DOC_DEPLOYMENT"))
            @warn "Setting fictive PETSc path because of RegistryCI or Doc deployment"
            PETSC_LIB = "JULIA_REGISTRYCI_AUTOMERGE"
        else
            throw(ErrorException("PETSc shared library (libpetsc.so) not found. Please check that PETSC_DIR and PETSC_ARCH env. variables are set."))
        end
    end

    return PETSC_LIB
end

# Absolute path to libpetsc.so
const libpetsc = get_petsc_location()


text:
const Region1_pt = (
n = [0.146_329_712_131_67,
-0.845_481_871_691_14,
-0.375_636_036_720_40E1,
 0.338_551_691_683_85E1,
-0.957_919_633_878_72,
 0.157_720_385_132_28,
-0.166_164_171_995_01E-1,
 0.812_146_299_835_68E-3,
 0.283_190_801_238_04E-3,
-0.607_063_015_658_74E-3,
-0.189_900_682_184_19E-1,
-0.325_297_487_705_05E-1,
-0.218_417_171_754_14E-1,
-0.528_383_579_699_30E-4,
-0.471_843_210_732_67E-3,
-0.300_017_807_930_26E-3,
 0.476_613_939_069_87E-4,
-0.441_418_453_308_46E-5,
-0.726_949_962_975_94E-15,
-0.316_796_448_450_54E-4,
-0.282_707_979_853_12E-5,
-0.852_051_281_201_03E-9,
-0.224_252_819_080_00E-5,
-0.651_712_228_956_01E-6,
-0.143_417_299_379_24E-12,
-0.405_169_968_601_17E-6,
-0.127_343_017_416_41E-8,
-0.174_248_712_306_34E-9,
-0.687_621_312_955_31E-18,
 0.144_783_078_285_21E-19,
 0.263_357_816_627_95E-22,
-0.119_476_226_400_71E-22,
 0.182_280_945_814_04E-23,
-0.935_370_872_924_58E-25],

I = [0,
 0,
 0,
 0,
 0,
 0,
 0,
 0,
 1,
 1,
 1,
 1,
 1,
 1,
 2,
 2,
 2,
 2,
 2,
 3,
 3,
 3,
 4,
 4,
 4,
 5,
 8,
 8,
 21,
 23,
 29,
 30,
 31,
 32],

J = [-2,
 -1,
 0,
 1,
 2,
 3,
 4,
 5,
 -9,
 -7,
 -1,
 0,
 1,
 3,
 -3,
 0,
 1,
 3,
 17,
 -4,
 0,
 6,
 -5,
 -2,
 10,
 -8,
 -11,
 -6,
 -29,
 -31,
 -38,
 -39,
 -40,
 -41])
"""
    IF97Region{:r1} <: ThermoModel

Returns all the property values in region 1.
273.15K ≤ T ≤ 623.15K  Psat(T) ≤ P ≤ 100MPa
Pressures in MPa and temperature in [K]
"""
const IF97Region1 =  IF97Region{:r1}

function mass_gibbs_impl(mt::SinglePT,model::IF97Region{:r1},p,t)
    P = convert_unit(u"Pa",u"MPa",p)
    T = normalize_units(t)
    Pstar = 16.53   #MPa
    Tstar = 1386.0  #K
    n = Region1_pt.n
    I = Region1_pt.I
    J = Region1_pt.J
    π = P / Pstar
    τ = Tstar / T
    ππ = 7.1 - π
    ττ = τ - 1.222
    γ    = sum(n[i]*(ππ^I[i])*(ττ^J[i]) for i=1:34)
    res =  IF97_R*T*γ
    return convert_unit(u"kJ/kg",u"J/kg",res)
end

function mass_helmholtz_impl(mt::SinglePT,model::IF97Region{:r1},p,t)
    P = convert_unit(u"Pa",u"MPa",p)
    T = normalize_units(t)
    Pstar = 16.53   #MPa
    Tstar = 1386.0  #K
    n = Region1_pt.n
    I = Region1_pt.I
    J = Region1_pt.J
    π = P / Pstar
    τ = Tstar / T
    ππ = 7.1 - π
    ττ = τ - 1.222
    γ    = sum(n[i]*(ππ^I[i])*(ττ^J[i]) for i=1:34)
    γ_π  = sum(-n[i]*I[i]*(ππ^(I[i]-1))*(ττ^J[i]) for i=1:34)

    res =  IF97_R*T*(γ - π*γ_π/1000)
    return convert_unit(u"kJ/kg",u"J/kg",res)
end

function mass_volume_impl(mt::SinglePT,model::IF97Region{:r1},p,t)
    P = convert_unit(u"Pa",u"MPa",p)
    T = normalize_units(t)
    Pstar = 16.53   #MPa
    Tstar = 1386.0  #K
    n = Region1_pt.n
    I = Region1_pt.I
    J = Region1_pt.J
    π = P / Pstar
    τ = Tstar / T
    ππ = 7.1 - π
    ττ = τ - 1.222
    γ_π  = sum(-n[i]*I[i]*(ππ^(I[i]-1))*(ττ^J[i]) for i=1:34)

    res = IF97_R*T*π*γ_π/P/1000
    return res
end

function mass_internal_energy_impl(mt::SinglePT,model::IF97Region{:r1},p,t)
    P = convert_unit(u"Pa",u"MPa",p)
    T = normalize_units(t)
    Pstar = 16.53   #MPa
    Tstar = 1386.0  #K
    n = Region1_pt.n
    I = Region1_pt.I
    J = Region1_pt.J
    π = P / Pstar
    τ = Tstar / T
    ππ = 7.1 - π
    ττ = τ - 1.222
    γ_π  = sum(-n[i]*I[i]*(ππ^(I[i]-1))*(ττ^J[i]) for i=1:34)
    γ_τ  = sum(n[i]*(ππ^I[i])*J[i]*(ττ^(J[i]-1)) for i=1:34)
    res =  IF97_R*T*(τ*γ_τ - π*γ_π)
    return convert_unit(u"kJ/kg",u"J/kg",res)
end

function mass_entropy_impl(mt::SinglePT,model::IF97Region{:r1},p,t)
    P = convert_unit(u"Pa",u"MPa",p)
    T = normalize_units(t)
    Pstar = 16.53   #MPa
    Tstar = 1386.0  #K
    n = Region1_pt.n
    I = Region1_pt.I
    J = Region1_pt.J
    π = P / Pstar
    τ = Tstar / T
    ππ = 7.1 - π
    ττ = τ - 1.222
    γ    = sum(n[i]*(ππ^I[i])*(ττ^J[i]) for i=1:34)
    γ_τ  = sum(n[i]*(ππ^I[i])*J[i]*(ττ^(J[i]-1)) for i=1:34)
    res =  IF97_R*(τ*γ_τ - γ)
    return convert_unit(u"kJ/(kg*K)",u"J/(kg*K)",res)
end

function mass_enthalpy_impl(mt::SinglePT,model::IF97Region{:r1},p,t)
    P = convert_unit(u"Pa",u"MPa",p)
    T = normalize_units(t)
    Pstar = 16.53   #MPa
    Tstar = 1386.0  #K
    n = Region1_pt.n
    I = Region1_pt.I
    J = Region1_pt.J
    π = P / Pstar
    τ = Tstar / T
    ππ = 7.1 - π
    ττ = τ - 1.222
    γ    = sum(n[i]*(ππ^I[i])*(ττ^J[i]) for i=1:34)
    γ_τ  = sum(n[i]*(ππ^I[i])*J[i]*(ττ^(J[i]-1)) for i=1:34)

    res = IF97_R*T*τ*γ_τ
    return convert_unit(u"kJ/kg",u"J/kg",res)
end

function mass_cp_impl(mt::SinglePT,model::IF97Region{:r1},p,t)
    P = convert_unit(u"Pa",u"MPa",p)
    T = normalize_units(t)
    Pstar = 16.53   #MPa
    Tstar = 1386.0  #K
    n = Region1_pt.n
    I = Region1_pt.I
    J = Region1_pt.J
    π = P / Pstar
    τ = Tstar / T
    ππ = 7.1 - π
    ττ = τ - 1.222
    γ_ττ = sum(n[i]*(ππ^I[i])*J[i]*(J[i]-1)*(ττ^(J[i]-2)) for i=1:34)

    res = return -IF97_R*τ^2*γ_ττ
    return convert_unit(u"kJ/(kg*K)",u"J/(kg*K)",res)

end

function mass_cv_impl(mt::SinglePT,model::IF97Region{:r1},p,t)
    P = convert_unit(u"Pa",u"MPa",p)
    T = normalize_units(t)
    Pstar = 16.53   #MPa
    Tstar = 1386.0  #K
    n = Region1_pt.n
    I = Region1_pt.I
    J = Region1_pt.J
    π = P / Pstar
    τ = Tstar / T
    ππ = 7.1 - π
    ττ = τ - 1.222
    γ_ππ = sum(n[i]*I[i]*(I[i]-1)*(ππ^(I[i]-2))*(ττ^J[i]) for i=1:34)
    γ_πτ = sum(-n[i]*I[i]*(ππ^(I[i]-1))*J[i]*(ττ^(J[i]-1)) for i=1:34)
    γ_π  = sum(-n[i]*I[i]*(ππ^(I[i]-1))*(ττ^J[i]) for i=1:34)
    γ_ττ = sum(n[i]*(ππ^I[i])*J[i]*(J[i]-1)*(ττ^(J[i]-2)) for i=1:34)

    res =  IF97_R*((-τ^2)*γ_ττ + (γ_π - τ*γ_πτ)^2/γ_ππ)
    return convert_unit(u"kJ/(kg*K)",u"J/(kg*K)",res)

end

function sound_speed_impl(mt::SinglePT,model::IF97Region{:r1},p,t)
    P = convert_unit(u"Pa",u"MPa",p)
    T = normalize_units(t)
    Pstar = 16.53   #MPa
    Tstar = 1386.0  #K
    n = Region1_pt.n
    I = Region1_pt.I
    J = Region1_pt.J
    π = P / Pstar
    τ = Tstar / T
    ππ = 7.1 - π
    ττ = τ - 1.222
    γ_ππ = sum(n[i]*I[i]*(I[i]-1)*(ππ^(I[i]-2))*(ττ^J[i]) for i=1:34)
    γ_ττ = sum(n[i]*(ππ^I[i])*J[i]*(J[i]-1)*(ττ^(J[i]-2)) for i=1:34)
    γ_πτ = sum(-n[i]*I[i]*(ππ^(I[i]-1))*J[i]*(ττ^(J[i]-1)) for i=1:34)
    γ_π  = sum(-n[i]*I[i]*(ππ^(I[i]-1))*(ττ^J[i]) for i=1:34)

    res =  return √(1000*IF97_R*T*(γ_π)^2/((γ_π - τ*γ_πτ)^2/τ^2/γ_ττ - γ_ππ))
    return res
end

const Region1_t_ph =(
n = [-0.238_724_899_245_21E3,
0.404_211_886_379_45E3,
0.113_497_468_817_18E3,
-0.584_576_160_480_39E1,
-0.152_854_824_131_40E-3,
-0.108_667_076_953_77E-5,
-0.133_917_448_726_02E2,
0.432_110_391_835_59E2,
-0.540_100_671_705_06E2,
0.305_358_922_039_16E2,
-0.659_647_494_236_38E1,
0.939_654_008_783_63E-2,
0.115_736_475_053_40E-6,
-0.258_586_412_820_73E-4,
-0.406_443_630_847_99E-8,
0.664_561_861_916_35E-7,
0.806_707_341_030_27E-10,
-0.934_777_712_139_47E-12,
0.582_654_420_206_01E-14,
-0.150_201_859_535_03E-16],

I = [0,
0,
0,
0,
0,
0,
1,
1,
1,
1,
1,
1,
1,
2,
2,
3,
3,
4,
5,
6],

J = [0,
1,
2,
6,
22,
32,
0,
1,
2,
3,
4,
10,
32,
10,
32,
10,
32,
32,
32,
32])

function temperature_impl(mt::SinglePH,model::IF97Region{:r1},p,_h)
    P = convert_unit(u"Pa",u"MPa",p)
    h = convert_unit(u"J/kg",u"kJ/kg",normalize_units(_h))
    hstar = 2500 #kJ/kg

    n = Region1_t_ph.n
    I = Region1_t_ph.I
    J = Region1_t_ph.J
    η = h / hstar
    return sum(n[i]*P^I[i]*(η+1)^J[i] for i=1:20)
end

const Region1_t_ps =(
    n = [0.174_782_680_583_07E3,
    0.348_069_308_928_73E2,
    0.652_925_849_784_55E1,
    0.330_399_817_754_89,
   -0.192_813_829_231_96E-6,
   -0.249_091_972_445_73E-22,
   -0.261_076_364_893_32,
    0.225_929_659_815_86,
   -0.642_564_633_952_26E-1,
    0.788_762_892_705_26E-2,
    0.356_721_106_073_66E-9,
    0.173_324_969_948_95E-23,
    0.566_089_006_548_37E-3,
   -0.326_354_831_397_17E-3,
    0.447_782_866_906_32E-4,
   -0.513_221_569_085_07E-9,
   -0.425_226_570_422_07E-25,
    0.264_004_413_606_89E-12,
    0.781_246_004_597_23E-28,
   -0.307_321_999_036_68E-30],

I = [0,
    0,
    0,
    0,
    0,
    0,
    1,
    1,
    1,
    1,
    1,
    1,
    2,
    2,
    2,
    2,
    2,
    3,
    3,
    4],

J = [0,
    1,
    2,
    3,
    11,
    31,
    0,
    1,
    2,
    3,
    12,
    31,
    0,
    1,
    2,
    9,
    31,
    10,
    32,
    32]

)

function temperature_impl(mt::SinglePS,model::IF97Region{:r1},p,_s)
    P = convert_unit(u"Pa",u"MPa",p)
    s = convert_unit(u"J/(kg*K)",u"kJ/(kg*K)",normalize_units(_s))
    n = Region1_t_ps.n
    I = Region1_t_ps.I
    J = Region1_t_ps.J

    return sum(n[i]*P^I[i]*(s+2)^J[i] for i=1:20)
end


text:
using CompileBot

botconfig = BotConfig(
  "MatLang";
  yml_path = "SnoopCompile.yml",
)

snoop_bot(
  botconfig,
  "$(@__DIR__)/example_script.jl",
)

max_stars_repo_name: juliamatlab/MatLang

max_issues_repo_name: juliamatlab/MatLang

max_forks_repo_name: juliamatlab/MatLang


text:
module Benchmarks

using BenchmarkTools, MathOptInterface

const MOI = MathOptInterface
const BENCHMARKS = Dict{String,Function}()

"""
    suite(
        new_model::Function;
        exclude::Vector{Regex} = Regex[]
    )

Create a suite of benchmarks. `new_model` should be a function that takes no
arguments, and returns a new instance of the optimizer you wish to benchmark.

Use `exclude` to exclude a subset of benchmarks.

### Examples

```julia
suite() do
    GLPK.Optimizer()
end
suite(exclude = [r"delete"]) do
    Gurobi.Optimizer(OutputFlag=0)
end
```
"""
function suite(new_model::Function; exclude::Vector{Regex} = Regex[])
    group = BenchmarkGroup()
    for (name, func) in BENCHMARKS
        any(occursin.(exclude, Ref(name))) && continue
        group[name] = @benchmarkable $func($new_model)
    end
    return group
end

"""
    create_baseline(suite, name::String; directory::String = ""; kwargs...)

Run all benchmarks in `suite` and save to files called `name` in `directory`.

Extra `kwargs` are based to `BenchmarkTools.run`.

### Examples

```julia
my_suite = suite(() -> GLPK.Optimizer())
create_baseline(my_suite, "glpk_master"; directory = "/tmp", verbose = true)
```
"""
function create_baseline(
    suite::BenchmarkTools.BenchmarkGroup,
    name::String;
    directory::String = "",
    kwargs...,
)
    tune!(suite)
    BenchmarkTools.save(
        joinpath(directory, name * "_params.json"),
        params(suite),
    )
    results = run(suite; kwargs...)
    BenchmarkTools.save(joinpath(directory, name * "_baseline.json"), results)
    return
end

"""
    compare_against_baseline(
        suite, name::String; directory::String = "",
        report_filename::String = "report.txt"
    )

Run all benchmarks in `suite` and compare against files called `name` in
`directory` that were created by a call to `create_baseline`.

A report summarizing the comparison is written to `report_filename` in
`directory`.

Extra `kwargs` are based to `BenchmarkTools.run`.

### Examples

```julia
my_suite = suite(() -> GLPK.Optimizer())
compare_against_baseline(
    my_suite, "glpk_master"; directory = "/tmp", verbose = true
)
```
"""
function compare_against_baseline(
    suite::BenchmarkTools.BenchmarkGroup,
    name::String;
    directory::String = "",
    report_filename::String = "report.txt",
    kwargs...,
)
    params_filename = joinpath(directory, name * "_params.json")
    baseline_filename = joinpath(directory, name * "_baseline.json")
    if !isfile(params_filename) || !isfile(baseline_filename)
        error("You create a baseline with `create_baseline` first.")
    end
    loadparams!(
        suite,
        BenchmarkTools.load(params_filename)[1],
        :evals,
        :samples,
    )
    new_results = run(suite; kwargs...)
    old_results = BenchmarkTools.load(baseline_filename)[1]
    open(joinpath(directory, report_filename), "w") do io
        println(stdout, "\n========== Results ==========")
        println(io, "\n========== Results ==========")
        for key in keys(new_results)
            judgement = judge(
                BenchmarkTools.median(new_results[key]),
                BenchmarkTools.median(old_results[key]),
            )
            println(stdout, "\n", key)
            println(io, "\n", key)
            show(stdout, MIME"text/plain"(), judgement)
            show(io, MIME"text/plain"(), judgement)
        end
    end
    return
end

###
### Benchmarks
###

macro add_benchmark(f)
    name = f.args[1].args[1]
    return quote
        $(esc(f))
        BENCHMARKS[String($(Base.Meta.quot(name)))] = $(esc(name))
    end
end

@add_benchmark function add_variable(new_model)
    model = new_model()
    for i in 1:10_000
        MOI.add_variable(model)
    end
    return model
end

@add_benchmark function add_variables(new_model)
    model = new_model()
    MOI.add_variables(model, 10_000)
    return model
end

@add_benchmark function add_variable_constraint(new_model)
    model = new_model()
    x = MOI.add_variables(model, 10_000)
    for (i, xi) in enumerate(x)
        MOI.add_constraint(model, MOI.SingleVariable(xi), MOI.LessThan(1.0 * i))
    end
    return model
end

@add_benchmark function add_variable_constraints(new_model)
    model = new_model()
    x = MOI.add_variables(model, 10_000)
    MOI.add_constraints(
        model,
        MOI.SingleVariable.(x),
        MOI.LessThan.(1.0:10_000.0),
    )
    return model
end

@add_benchmark function delete_variable(new_model)
    model = new_model()
    x = MOI.add_variables(model, 1_000)
    MOI.add_constraint.(model, MOI.SingleVariable.(x), Ref(MOI.LessThan(1.0)))
    MOI.delete.(model, x)
    return model
end

@add_benchmark function delete_variable_constraint(new_model)
    model = new_model()
    x = MOI.add_variables(model, 1_000)
    cons =
        MOI.add_constraint.(
            model,
            MOI.SingleVariable.(x),
            Ref(MOI.LessThan(1.0)),
        )
    for con in cons
        MOI.delete(model, con)
    end
    cons =
        MOI.add_constraint.(
            model,
            MOI.SingleVariable.(x),
            Ref(MOI.LessThan(1.0)),
        )
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, x), 0.0),
    )
    MOI.optimize!(model)
    for con in cons
        MOI.delete(model, con)
    end
    return model
end

@add_benchmark function add_constraint(new_model)
    model = new_model()
    index = MOI.add_variables(model, 10_000)
    for (i, x) in enumerate(index)
        MOI.add_constraint(
            model,
            MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x)], 0.0),
            MOI.LessThan(1.0 * i),
        )
    end
    return model
end

@add_benchmark function add_constraints(new_model)
    model = new_model()
    x = MOI.add_variables(model, 10_000)
    MOI.add_constraints(
        model,
        [
            MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, xi)], 0.0) for
            xi in x
        ],
        MOI.LessThan.(1:1.0:10_000),
    )
    return model
end

@add_benchmark function delete_constraint(new_model)
    model = new_model()
    index = MOI.add_variables(model, 1_000)
    cons = Vector{
        MOI.ConstraintIndex{
            MOI.ScalarAffineFunction{Float64},
            MOI.LessThan{Float64},
        },
    }(
        undef,
        1_000,
    )
    for (i, x) in enumerate(index)
        cons[i] = MOI.add_constraint(
            model,
            MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x)], 0.0),
            MOI.LessThan(1.0 * i),
        )
    end
    for con in cons
        MOI.delete(model, con)
    end
    return model
end

@add_benchmark function copy_model(new_model)
    model = new_model()
    index = MOI.add_variables(model, 1_000)
    cons = Vector{
        MOI.ConstraintIndex{
            MOI.ScalarAffineFunction{Float64},
            MOI.LessThan{Float64},
        },
    }(
        undef,
        1_000,
    )
    for (i, x) in enumerate(index)
        cons[i] = MOI.add_constraint(
            model,
            MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x)], 0.0),
            MOI.LessThan(1.0 * i),
        )
    end

    model2 = new_model()
    MOI.copy_to(model2, model)
    # MOI.copy_to(model2, model, filter_constraints=(x) -> x in cons[1:500])

    return model2
end

end


text:
module MetNets

    import MAT
    import SparseArrays
    import SparseArrays: AbstractSparseArray, sparse
    import Statistics: mean

    using ProgressMeter
    using Serialization

    # Types
    include("Types/Rxns.jl")
    include("Types/Mets.jl")
    include("Types/MetNet.jl")
    include("Types/MetStates.jl")

    # Utils
    include("Utils/iders.jl") # Need to be included first
    include("Utils/balance_str.jl")
    include("Utils/base.jl")
    include("Utils/check_dims.jl")
    include("Utils/clampfields.jl")
    include("Utils/compressed_model.jl")
    include("Utils/defaults.jl")
    include("Utils/del_blocked.jl")
    include("Utils/del_met.jl")
    include("Utils/del_rxn.jl")
    include("Utils/del_void_rxns_mets.jl")
    include("Utils/extended_model.jl")
    include("Utils/fixxing.jl")
    include("Utils/force_dims.jl")
    include("Utils/getter.jl")
    include("Utils/interchange_rxn.jl")
    include("Utils/invert_bkwds.jl")
    include("Utils/invert_rxn.jl")
    include("Utils/make_intake_info_compat.jl")
    include("Utils/queries.jl")
    include("Utils/read_mat.jl")
    include("Utils/rxn_str.jl")
    include("Utils/scale_polytope.jl")
    include("Utils/search.jl")
    include("Utils/set_met.jl")
    include("Utils/set_rxn.jl")
    include("Utils/setter.jl")
    include("Utils/similar_rxns.jl")
    include("Utils/sparse.jl")
    include("Utils/split_revs.jl")
    include("Utils/stoi_err.jl")
    include("Utils/summary.jl")
    include("Utils/TNmarginal.jl")
    include("Utils/update.jl")
    include("Utils/utils.jl")
    include("Utils/well_scaled_model.jl")

    # Test
    include("Test/ecoli_core.jl")
    include("Test/lineal_model.jl")
    include("Test/toy_model.jl")

end


text:
using Documenter, ReactionMechanismSimulator

makedocs(
    format = Documenter.HTML(),
    sitename = "ReactionMechanismSimulator.jl",
    modules = [ReactionMechanismSimulator],
    pages = [
        "Home" => "index.md",
        "Input.md",
        "Simulating.md",
        "Analysis.md",
    ]
)

deploydocs(
    repo = "github.com/ReactionMechanismGenerator/ReactionMechanismSimulator.jl.git",
    target = "build",
)


text:
@doc raw"""
    DiffusionStateModel(f::Function, g::Function, init)

Returns a diffusion process hidden state model ``dX_t = f(X_t)dt + g(X_t)dW_t``, where ``f`` is the `drift_function`, ``g`` is the `observation_function`, ``X_t`` is the ``n``-dimensional hidden state at time ``t``, and ``W_t`` is an ``m``-dimensional Brownian motion process.

Argument `init` stands for the initial condition of the process, which is either
* A vector of length `n` for a fixed (deterministic) initial condition
* A `Distributions.Sampleable` type for a random initial condition
"""
struct DiffusionStateModel{S, F1, F2, TI} <: HiddenStateModel{Vector{S}, ContinuousTime}
    n::Int
    m::Int
    drift_function::F1
    diffusion_function::F2
    init::TI
    function DiffusionStateModel(f::Function, g::Function, init::TI) where TI<:Union{Distributions.Sampleable, Any}
        if init isa Distributions.Sampleable
            x = rand(init)
            x isa Vector || (x = [x])
        else
            x = init
        end

        F = f(x)
        G = g(x)

        if length(F) == length(x)
            n = length(F)
        else
            error("Error: output of drift_function and initial condition have incompatible lengths.")
        end

        if length(F) != size(G, 1)
            error("Error: drift_function and diffusion_function have incompatible output sizes.")
        end

        m = size(G, 2)

        if eltype(x) == eltype(F) == eltype(G)
            S = eltype(x)
            if x isa Vector
                return new{S, typeof(f), typeof(g), typeof(init)}(n, m, f, g, init)
            else
                error("Error: initial condition has incorrect type.")
            end
        else
            error("Error: initial condition and outputs of drift_function and diffusion_function have incompatible element types.")
        end
    end
end

#####################
### BASIC METHODS ###
#####################

function drift(model::DiffusionStateModel)
    if hasmethod(model.drift_function, Tuple{AbstractMatrix})
        return model.drift_function
    else
        f(x) = model.drift_function(x)
        f(x::AbstractMatrix) = mapslices(f, x, dims=1)
        return f
    end
end

function diffusion(model::DiffusionStateModel)
    if hasmethod(model.diffusion_function, Tuple{AbstractMatrix})
        return model.diffusion_function
    else
        g(x) = model.diffusion_function(x)
        g(x::AbstractMatrix) = cat([model.diffusion_function(col) for col in eachcol(x)]..., dims=3)
        return g
    end
end

initial_condition(model::DiffusionStateModel) = model.init

state_dim(model::DiffusionStateModel) = model.n

"""
    noise_dim(model)

Returns the dimension of the Brownian motion ``W_t`` in the diffusion model ``dX_t = f(X_t)dt + g(X_t)dW_t``.
"""
noise_dim(model::DiffusionStateModel) = model.m
function initialize(model::DiffusionStateModel{T, F1, F2, TI}) where {T, F1, F2, TI<:Distributions.Sampleable}
    x = rand(model.init)
    if x isa Vector
        return x
    else
        return [x]
    end
end

function initialize(model::DiffusionStateModel{T, F1, F2, TI}) where
    TI<:AbstractVector{T} where {T, F1, F2}
    return model.init
end

"""
    drift_function(model)

Returns the drift function ``f`` of the diffusion model ``dX_t = f(X_t)dt + g(X_t)dW_t``.
"""
drift_function(model::DiffusionStateModel) = drift(model)

"""
    diffusion_function(model)

Returns the diffusion function ``g`` of the diffusion model ``dX_t = f(X_t)dt + g(X_t)dW_t``.
"""
diffusion_function(model::DiffusionStateModel) = diffusion(model)

function Base.show(io::IO, ::MIME"text/plain", model::DiffusionStateModel)
    print(io, "Diffusion process model for the hidden state
    type of hidden state:                   ", state_dim(model),"-dimensional vector
    number of independent Brownian motions: ", noise_dim(model),"
    initial condition:                      ", model.init isa Distributions.Sampleable ? "random" : "fixed")
end

function Base.show(io::IO, model::DiffusionStateModel)
    print(io, state_dim(model),"-dimensional diffusion with ", noise_dim(model), "-dimensional Brownian motion")
end

function Base.show(io::IO, ::MIME"text/plain", model::Type{DiffusionStateModel{S1, S2, S3, S4}}) where {S1, S2, S3, S4}
    print(io, "DiffusionStateModel{", S1, ",", S2, ",", S3, ",", S4, "}")
end

function Base.show(io::IO, model::Type{DiffusionStateModel{S1, S2, S3, S4}}) where {S1, S2, S3, S4}
    print(io, "DiffusionStateModel")
end

######################
### TIME EVOLUTION ###
######################

(model::DiffusionStateModel)(x::AbstractVector, dt) = x + model.drift_function(x)*dt + sqrt(dt)*model.diffusion_function(x)*randn(eltype(x), noise_dim(model))

function (model::DiffusionStateModel)(x::AbstractMatrix, dt)
    N   = size(x, 2)
    F   = drift_function(model)(x)
    G   = diffusion_function(model)(x)
    if ndims(G) == 2
        G = repeat(G, 1, 1, N)
    end
    out = zeros(eltype(x), state_dim(model), N)
    sqr = sqrt(dt)
    @inbounds for b in 1:N
        for c in 1:noise_dim(model), a in 1:state_dim(model)
            out[a,b]  +=  x[a,b]  +  dt * F[a,b]  +  sqr * G[a,c,b] * randn(eltype(x))
        end
    end
    return out
end

function propagate!(x::AbstractVector, model::DiffusionStateModel, dt)
    N   = size(x, 2)
    F   = drift_function(model)(x)
    G   = diffusion_function(model)(x)
    sqr = sqrt(dt)
    @inbounds for c in 1:noise_dim(model), a in 1:state_dim(model)
            x[a]  +=  dt * F[a]  +  sqr * G[a,c] * randn(eltype(x))
    end
    return x
end

function propagate!(x::AbstractMatrix, model::DiffusionStateModel, dt)
    N   = size(x, 2)
    F   = drift_function(model)(x)
    G   = diffusion_function(model)(x)
    if ndims(G) == 2
        G = repeat(G, 1, 1, N)
    end
    sqr = sqrt(dt)
    @inbounds for b in 1:N
        for c in 1:noise_dim(model)
            for a in 1:state_dim(model)
                x[a,b]  +=  dt * F[a,b]  +  sqr * G[a,c,b] * randn(eltype(x))
            end
        end
    end
    return x
end

################################
### CONVENIENCE CONSTRUCTORS ###
################################

@doc raw"""
    ScalarDiffusionStateModel(f::Function, g::Function, init)

Returns a scalar diffusion process hidden state model ``dX_t = f(X_t)dt + g(X_t)dW_t``.
"""
function ScalarDiffusionStateModel(f::Function, g::Function, init::ContinuousUnivariateDistribution)
    F(x::AbstractVector) = [f(x[1])]
    F(x::AbstractMatrix) = mapslices(F, x, dims=1)
    G(x::AbstractVector) = g(x[1])*ones(1,1)
    G(x::AbstractMatrix) = permutedims(repeat(mapslices(y->g(y[1]), x, dims=1), 1, 1, 1), [1,3,2])
    return DiffusionStateModel(F, G, init)
end

max_stars_repo_name: simsurace/FeedbackParticleFilters

max_issues_repo_name: simsurace/FeedbackParticleFilters

max_forks_repo_name: simsurace/FeedbackParticleFilters


text:
using CUDAdrv
using LLVM

# TODO: import from CUDAapi.jl
import CUDAdrv: debug

## auxiliary routines

# helper methods for querying version DBs
# (tool::VersionNumber => devices::Vector{VersionNumber})
search(db, predicate) =
    Set(Base.Iterators.flatten(valvec for (key,valvec) in db if predicate(key)))

# database of LLVM versions and the supported devices
const llvm_db = [
    v"3.2" => [v"2.0", v"2.1", v"3.0", v"3.5"],
    v"3.5" => [v"5.0"],
    v"3.7" => [v"3.2", v"3.7", v"5.2", v"5.3"],
    v"3.9" => [v"6.0", v"6.1", v"6.2"]
]

# check support for the LLVM version
function check_llvm()
    llvm_version = LLVM.version()
    info("Using LLVM $llvm_version")

    InitializeAllTargets()
    "nvptx" in LLVM.name.(collect(targets())) ||
        error("Your LLVM does not support the NVPTX back-end. Fix this, and rebuild LLVM.jl and CUDAnative.jl.")

    llvm_support = search(llvm_db, ver -> ver <= llvm_version)
    isempty(llvm_support) && error("LLVM $llvm_version does not support any compatible device")

    return llvm_version, llvm_support
end

# database of CUDA versions and the supported devices
const cuda_db = [
    v"4.0" => [v"2.0", v"2.1"],
    v"4.2" => [v"3.0"],
    v"5.0" => [v"3.5"],
    v"6.0" => [v"3.2", v"5.0"],
    v"6.5" => [v"3.7"],
    v"7.0" => [v"5.2"],
    v"7.5" => [v"5.3"],
    v"8.0" => [v"6.0", v"6.1", v"6.2"]
]

# check support for the CUDA version
function check_cuda()
    cuda_version = CUDAdrv.version()
    info("Using CUDA $cuda_version")

    cuda_support = search(cuda_db, ver -> ver <= cuda_version)
    isempty(cuda_support) && error("CUDA $cuda_version does not support any compatible device")

    return cuda_version, cuda_support
end

# check support for the Julia version
function check_julia(llvm_version)
    julia_llvm_version = VersionNumber(Base.libllvm_version)
    info("Using Julia's LLVM $julia_llvm_version")

    if julia_llvm_version != llvm_version
        error("LLVM $llvm_version incompatible with Julia's LLVM $julia_llvm_version")
    end

    return julia_llvm_version
end

## discovery routines

# find CUDA toolkit
function find_cuda()
    cuda_envvars = ["CUDA_PATH", "CUDA_HOME", "CUDA_ROOT"]
    cuda_envvars_set = filter(var -> haskey(ENV, var), cuda_envvars)
    if length(cuda_envvars_set) > 0
        cuda_paths = unique(map(var->ENV[var], cuda_envvars_set))
        if length(unique(cuda_paths)) > 1
            warn("Multiple CUDA path environment variables set: $(join(cuda_envvars_set, ", ", " and ")). ",
                 "Arbitrarily selecting CUDA at $(first(cuda_paths)). ",
                 "To ensure a consistent path, ensure only a single unique CUDA path is set.")
        end
        cuda_path = Nullable(first(cuda_paths))
    else
        cuda_path = Nullable{String}()
    end

    return cuda_path
end

# find device library bitcode files
function find_libdevice(cuda_path, supported_capabilities)
    libraries = Dict{VersionNumber,String}()

    # find the root directory
    if haskey(ENV, "NVVMIR_LIBRARY_DIR")
        dirs = [ENV["NVVMIR_LIBRARY_DIR"]]
    elseif !isnull(cuda_path)
        dirs = ["$(get(cuda_path))/libdevice",
                "$(get(cuda_path))/nvvm/libdevice"]
    else
        dirs = ["/usr/lib/nvidia-cuda-toolkit/libdevice",
                "/usr/local/cuda/nvvm/libdevice",
                "/opt/cuda/nvvm/libdevice"]
    end
    dirs = unique(filter(isdir, dirs))
    if isempty(dirs)
        error("CUDA device library path not found. ",
              "Specify by setting a CUDA path variable, or by setting NVVMIR_LIBRARY_DIR.")
    elseif length(dirs) > 1
        warn("Multiple locations found with device code: $(join(dirs, ", ", " and ")). ",
             "Arbitrarily selecting those at $(first(dirs)).")
    end
    dir = first(dirs)

    # discover individual device library files
    for cap in supported_capabilities
        path = joinpath(dir, "libdevice.compute_$(cap.major)$(cap.minor).10.bc")
        if isfile(path)
            libraries[cap] = path
        end
    end
    isempty(libraries) && error("No device libraries found in $dir for your hardware.")
    info("Found libdevice for $(join(sort(map(ver->"sm_$(ver.major)$(ver.minor)", keys(libraries))), ", ", " and "))")

    return libraries
end

## main

const ext = joinpath(@__DIR__, "ext.jl")

function main()
    # check version support
    llvm_version, llvm_support = check_llvm()
    cuda_version, cuda_support = check_cuda()
    julia_llvm_version = check_julia(llvm_version)

    # figure out supported capabilities
    supported_capabilities = Vector{VersionNumber}()
    append!(supported_capabilities, llvm_support ∩ cuda_support)
    debug("Supported capabilities: $(join(supported_capabilities, ", "))")

    # discover stuff
    cuda_path = find_cuda()
    libdevice_libraries = find_libdevice(cuda_path, supported_capabilities)

    # check if we need to rebuild
    if isfile(ext)
        debug("Checking validity of existing ext.jl...")
        @eval module Previous; include($ext); end
        if  isdefined(Previous, :cuda_version)           && Previous.cuda_version == cuda_version &&
            isdefined(Previous, :llvm_version)           && Previous.llvm_version == llvm_version &&
            isdefined(Previous, :julia_llvm_version)     && Previous.julia_llvm_version == julia_llvm_version &&
            isdefined(Previous, :supported_capabilities) && Previous.supported_capabilities == supported_capabilities &&
            isdefined(Previous, :libdevice_libraries)    && Previous.libdevice_libraries == libdevice_libraries
            info("CUDAnative.jl has already been built for this set-up, no need to rebuild")
            return
        end
    end

    # write ext.jl
    open(ext, "w") do fh
        write(fh, """
            # Toolchain properties
            const cuda_version = v"$cuda_version"
            const llvm_version = v"$llvm_version"
            const julia_llvm_version = v"$julia_llvm_version"

            const supported_capabilities = $supported_capabilities
            const libdevice_libraries = $libdevice_libraries
            """)
    end
    nothing
end

try
    main()
catch ex
    # if anything goes wrong, wipe the existing ext.jl to prevent the package from loading
    rm(ext; force=true)
    rethrow(ex)
end