text:
# Defining new Python types from Julia (ugly simulation of C headers)

################################################################
# Python expects the PyMethodDef and similar strings to be constants,
# so we define anonymous globals to hold them, returning the pointer
const permanent_strings = String[]
function gstring_ptr(s::AbstractString)
    g = String(s)
    push!(permanent_strings, g)
    unsafe_convert(Ptr{UInt8}, g)
end
gstring_ptr(s::Ptr) = convert(Ptr{UInt8}, s)

gstring_ptr_ornull(s::AbstractString) =
    isempty(s) ? NULL_UInt8_Ptr : gstring_ptr(s)
gstring_ptr_ornull(s::Ptr) = gstring_ptr(s)

################################################################
# mirror of Python API types and constants from methodobject.h

pymethoddef(name=C_NULL, meth=C_NULL, flags=0)

struct PyMethodDef
    ml_name::Ptr{UInt8}
    ml_meth::Ptr{Cvoid}
    ml_flags::Cint
    ml_doc::Ptr{UInt8}
    function PyMethodDef(name=C_NULL, meth=C_NULL, flags=0, doc=C_NULL)
        new(gstring_ptr(name), convert(Ptr{Cvoid}, meth), convert(Cint, flags), gstring_ptr_ornull(doc))
    end
end

const NULL_UInt8_Ptr = convert(Ptr{UInt8}, C_NULL)

################################################################
# mirror of Python API types and constants from descrobject.h

struct PyGetSetDef
    name::Ptr{UInt8}
    get::Ptr{Cvoid}
    set::Ptr{Cvoid} # may be NULL for read-only members
    doc::Ptr{UInt8} # may be NULL
    closure::Ptr{Cvoid} # pass-through thunk, may be NULL
    function PyGetSetDef(_name=C_NULL, _get=C_NULL, _set=C_NULL, _doc=C_NULL,_closure=C_NULL; name=_name, get=_get, set=_set, doc=_doc, closure=_closure)
        new(gstring_ptr(name), convert(Ptr{Cvoid}, get), convert(Ptr{Cvoid}, set), gstring_ptr_ornull(doc), convert(Ptr{Cvoid}, closure))
    end
end

################################################################
# from Python structmember.h:

# declare immutable because we need a C-like array of these
struct PyMemberDef
    name::Ptr{UInt8}
    typ::Cint
    offset::Int # warning: was Cint for Python <= 2.4
    flags::Cint
    doc::Ptr{UInt8}
    function PyMemberDef(name=C_NULL,typ=0,offset=0,flags=0,doc=C_NULL)
        new(gstring_ptr(name),
            convert(Cint,typ),
            convert(Int,offset),
            convert(Cint,flags),
            gstring_ptr_ornull(doc))
    end
end

################################################################
# Mirror of PyNumberMethods in Python object.h

const PyNumberMethods_fields = [
     (:nb_add, Ptr{Cvoid}, C_NULL),
     (:nb_subtract, Ptr{Cvoid}, C_NULL),
     (:nb_multiply, Ptr{Cvoid}, C_NULL),
     (:nb_remainder, Ptr{Cvoid}, C_NULL),
     (:nb_divmod, Ptr{Cvoid}, C_NULL),
     (:nb_power, Ptr{Cvoid}, C_NULL),
     (:nb_negative, Ptr{Cvoid}, C_NULL),
     (:nb_positive, Ptr{Cvoid}, C_NULL),
     (:nb_absolute, Ptr{Cvoid}, C_NULL),
     (:nb_bool, Ptr{Cvoid}, C_NULL),
     (:nb_invert, Ptr{Cvoid}, C_NULL),
     (:nb_lshift, Ptr{Cvoid}, C_NULL),
     (:nb_rshift, Ptr{Cvoid}, C_NULL),
     (:nb_and, Ptr{Cvoid}, C_NULL),
     (:nb_xor, Ptr{Cvoid}, C_NULL),
     (:nb_or, Ptr{Cvoid}, C_NULL),
     (:nb_int, Ptr{Cvoid}, C_NULL),
     (:nb_reserved, Ptr{Cvoid}, C_NULL),
     (:nb_float, Ptr{Cvoid}, C_NULL),
     (:nb_inplace_add, Ptr{Cvoid}, C_NULL),
     (:nb_inplace_subtract, Ptr{Cvoid}, C_NULL),
     (:nb_inplace_multiply, Ptr{Cvoid}, C_NULL),
     (:nb_inplace_remainder, Ptr{Cvoid}, C_NULL),
     (:nb_inplace_power, Ptr{Cvoid}, C_NULL),
     (:nb_inplace_lshift, Ptr{Cvoid}, C_NULL),
     (:nb_inplace_rshift, Ptr{Cvoid}, C_NULL),
     (:nb_inplace_and, Ptr{Cvoid}, C_NULL),
     (:nb_inplace_xor, Ptr{Cvoid}, C_NULL),
     (:nb_inplace_or, Ptr{Cvoid}, C_NULL),
     (:nb_floordivide, Ptr{Cvoid}, C_NULL),
     (:nb_truedivide, Ptr{Cvoid}, C_NULL),
     (:nb_inplace_floordivide, Ptr{Cvoid}, C_NULL),
     (:nb_inplace_truedivide, Ptr{Cvoid}, C_NULL),
     (:nb_index, Ptr{Cvoid}, C_NULL),
     (:nb_matrixmultiply, Ptr{Cvoid}, C_NULL),
     (:nb_imatrixmultiply, Ptr{Cvoid}, C_NULL),
]

@eval struct PyNumberMethods
    $([:($n::$t) for (n,t,d) in PyNumberMethods_fields]...)
    PyNumberMethods(; $([Expr(:kw, n, d) for (n,t,d) in PyNumberMethods_fields]...)) =
        new($([:(convert($t, $n)) for (n,t,d) in PyNumberMethods_fields]...))
end

################################################################
# Mirror of PySequenceMethods in Python object.h

const PySequenceMethods_fields = [
    (:sq_length, Ptr{Cvoid}, C_NULL),
    (:sq_concat, Ptr{Cvoid}, C_NULL),
    (:sq_repeat, Ptr{Cvoid}, C_NULL),
    (:sq_item, Ptr{Cvoid}, C_NULL),
    (:was_sq_item, Ptr{Cvoid}, C_NULL),
    (:sq_ass_item, Ptr{Cvoid}, C_NULL),
    (:was_sq_ass_slice, Ptr{Cvoid}, C_NULL),
    (:sq_contains, Ptr{Cvoid}, C_NULL),
    (:sq_inplace_concat, Ptr{Cvoid}, C_NULL),
    (:sq_inplace_repeat, Ptr{Cvoid}, C_NULL),
]

@eval struct PySequenceMethods
    $([:($n :: $t) for (n,t,d) in PySequenceMethods_fields]...)
    PySequenceMethods(; $([Expr(:kw, n, d) for (n,t,d) in PySequenceMethods_fields]...)) =
        new($([:(convert($t, $n)) for (n,t,d) in PySequenceMethods_fields]...))
end

################################################################
# Mirror of PyMappingMethods in Python object.h

const PyMappingMethods_fields = [
    (:mp_length, Ptr{Cvoid}, C_NULL),
    (:mp_subscript, Ptr{Cvoid}, C_NULL),
    (:mp_ass_subscript, Ptr{Cvoid}, C_NULL),
]

@eval struct PyMappingMethods
    $([:($n :: $t) for (n,t,d) in PyMappingMethods_fields]...)
    PyMappingMethods(; $([Expr(:kw, n, d) for (n,t,d) in PyMappingMethods_fields]...)) =
        new($([:(convert($t, $n)) for (n,t,d) in PyMappingMethods_fields]...))
end

################################################################
# Mirror of PyTypeObject in Python object.h
#  -- assumes non-debugging Python build (no Py_TRACE_REFS)
#  -- most fields can default to 0 except where noted

PyTypeObject_defaultflags() =
    pyversion.major >= 3 ?
        (Py_TPFLAGS_HAVE_STACKLESS_EXTENSION[] |
         Py_TPFLAGS_HAVE_VERSION_TAG) :
        (Py_TPFLAGS_HAVE_GETCHARBUFFER |
         Py_TPFLAGS_HAVE_SEQUENCE_IN |
         Py_TPFLAGS_HAVE_INPLACEOPS |
         Py_TPFLAGS_HAVE_RICHCOMPARE |
         Py_TPFLAGS_HAVE_WEAKREFS |
         Py_TPFLAGS_HAVE_ITER |
         Py_TPFLAGS_HAVE_CLASS |
         Py_TPFLAGS_HAVE_STACKLESS_EXTENSION[] |
         Py_TPFLAGS_HAVE_INDEX)

const PyTypeObject_fields = [
    # PyObject_HEAD (for non-Py_TRACE_REFS build):
    (:ob_refcnt, Int, 0),
    (:ob_type, PyPtr, C_NULL),
    (:ob_size, Int, 0), # PyObject_VAR_HEAD, C_NULL

    # PyTypeObject fields:
    (:tp_name, Ptr{UInt8}, C_NULL), # required, should be in format "<module>.<name>"

    # warning: these two were Cint for Python <= 2.4
    (:tp_basicsize, Int, 0), # required, = sizeof(instance)
    (:tp_itemsize, Int, 0),

    (:tp_dealloc, Ptr{Cvoid}, C_NULL),
    (:tp_print, Ptr{Cvoid}, C_NULL),
    (:tp_getattr, Ptr{Cvoid}, C_NULL),
    (:tp_setattr, Ptr{Cvoid}, C_NULL),
    (:tp_compare, Ptr{Cvoid}, C_NULL),
    (:tp_repr, Ptr{Cvoid}, C_NULL),

    (:tp_as_number, Ptr{PyNumberMethods}, C_NULL),
    (:tp_as_sequence, Ptr{PySequenceMethods}, C_NULL),
    (:tp_as_mapping, Ptr{PyMappingMethods}, C_NULL),

    (:tp_hash, Ptr{Cvoid}, C_NULL),
    (:tp_call, Ptr{Cvoid}, C_NULL),
    (:tp_str, Ptr{Cvoid}, C_NULL),
    (:tp_getattro, Ptr{Cvoid}, C_NULL),
    (:tp_setattro, Ptr{Cvoid}, C_NULL),

    (:tp_as_buffer, Ptr{Cvoid}, C_NULL),

    (:tp_flags, Clong, 0), # Required, should default to Py_TPFLAGS_DEFAULT

    (:tp_doc, Ptr{UInt8}, C_NULL), # normally set in example code, but may be NULL

    (:tp_traverse, Ptr{Cvoid}, C_NULL),

    (:tp_clear, Ptr{Cvoid}, C_NULL),

    (:tp_richcompare, Ptr{Cvoid}, C_NULL),

    (:tp_weaklistoffset, Int, 0),

    # added in Python 2.2:
    (:tp_iter, Ptr{Cvoid}, C_NULL),
    (:tp_iternext, Ptr{Cvoid}, C_NULL),

    (:tp_methods, Ptr{PyMethodDef}, C_NULL),
    (:tp_members, Ptr{PyMemberDef}, C_NULL),
    (:tp_getset, Ptr{PyGetSetDef}, C_NULL),
    (:tp_base, Ptr{Cvoid}, C_NULL),

    (:tp_dict, PyPtr, C_NULL),
    (:tp_descr_get, Ptr{Cvoid}, C_NULL),
    (:tp_descr_set, Ptr{Cvoid}, C_NULL),
    (:tp_dictoffset, Int, 0),

    (:tp_init, Ptr{Cvoid}, C_NULL),
    (:tp_alloc, Ptr{Cvoid}, C_NULL),
    (:tp_new, Ptr{Cvoid}, C_NULL),
    (:tp_free, Ptr{Cvoid}, C_NULL),
    (:tp_is_gc, Ptr{Cvoid}, C_NULL),

    (:tp_bases, PyPtr, C_NULL),
    (:tp_mro, PyPtr, C_NULL),
    (:tp_cache, PyPtr, C_NULL),
    (:tp_subclasses, PyPtr, C_NULL),
    (:tp_weaklist, PyPtr, C_NULL),
    (:tp_del, Ptr{Cvoid}, C_NULL),

    # added in Python 2.6:
    (:tp_version_tag, Cuint, 0),

    # only used for COUNT_ALLOCS builds of Python
    (:tp_allocs, Int, 0),
    (:tp_frees, Int, 0),
    (:tp_maxalloc, Int, 0),
    (:tp_prev, Ptr{Cvoid}, C_NULL),
    (:tp_next, Ptr{Cvoid}, C_NULL),
]

@eval mutable struct PyTypeObject
    $([:($n :: $t) for (n,t,d) in PyTypeObject_fields]...)
    # cache of julia objects referenced by this type, to prevent them being garbage-collected
    jl_cache::Dict{Symbol,Any}

    function PyTypeObject(; unsafe_null=false, opts...)
        t = new($([:(convert($t, $d)) for (n,t,d) in PyTypeObject_fields]...), Dict{Symbol,Any}())
        unsafe_null ? t : PyTypeObject_init!(t; opts...)
    end
end

function PyTypeObject_init!(t::PyTypeObject; opts...)
    for (k, x) in pairs(opts)
        setproperty!(t, k, x)
    end
    t.tp_name == C_NULL && error("required: tp_name")
    t.tp_basicsize == 0 && !haskey(opts, :tp_basicsize) && error("required: tp_basicsize")
    t.tp_flags == 0 && !haskey(opts, :tp_flags) && (t.tp_flags = PyTypeObject_defaultflags())
    t.tp_new == C_NULL && !haskey(opts, :tp_new) && (t.tp_new = @pyglobal(:PyType_GenericNew))
    @pycheckz ccall((@pysym :PyType_Ready), Cint, (Ref{PyTypeObject},), t)
    ccall((@pysym :Py_IncRef), Cvoid, (Any,), t)
    return t
end

function Base.setproperty!(t::PyTypeObject, k::Symbol, x)
    if k == :tp_name && x isa AbstractString
        z = t.jl_cache[k] = Base.cconvert(Ptr{UInt8}, x)
        setfield!(t, k, unsafe_convert(Ptr{UInt8}, z))
    elseif k == :tp_as_number && x isa PyNumberMethods
        z = t.jl_cache[k] = Ref(x)
        setfield!(t, k, unsafe_convert(Ptr{PyNumberMethods}, z))
    elseif k == :tp_as_sequence && x isa PySequenceMethods
        z = t.jl_cache[k] = Ref(x)
        setfield!(t, k, unsafe_convert(Ptr{PySequenceMethods}, z))
    elseif k == :tp_as_mapping && x isa PyMappingMethods
        z = t.jl_cache[k] = Ref(x)
        setfield!(t, k, unsafe_convert(Ptr{PyMappingMethods}, z))
    elseif k == :tp_members && x isa AbstractVector{PyMemberDef}
        z = t.jl_cache[k] = push!(copy(convert(Vector{PyMemberDef}, x)), PyMemberDef())
        setfield!(t, k, pointer(z))
    elseif k == :tp_methods && x isa AbstractVector{PyMethodDef}
        z = t.jl_cache[k] = push!(copy(convert(Vector{PyMethodDef}, x)), PyMethodDef())
        setfield!(t, k, pointer(z))
    elseif k == :tp_getset && x isa AbstractVector{PyGetSetDef}
        z = t.jl_cache[k] = push!(copy(convert(Vector{PyGetSetDef}, x)), PyGetSetDef())
        setfield!(t, k, pointer(z))
    elseif k == :tp_dict && x isa NamedTuple
        d = t.jl_cache[k] = PyObject(@pycheckn ccall(@pysym(:PyDict_New), PyPtr, ()))
        for (k, v) in pairs(x)
            @pycheckz ccall(@pysym(:PyDict_SetItemString), Cint, (PyPtr, Cstring, PyPtr), d, string(k), PyObject(v))
        end
        setfield!(t, k, PyPtr(d))
    elseif k == :tp_base && x isa PyTypeObject
        t.jl_cache[k] = x
        setfield!(t, k, pointer_from_objref(x))
    else
        setfield!(t, k, convert(fieldtype(PyTypeObject, k), x))
    end
end

unsafe_pytype(o::PyPtr) =
    convert(Ptr{PyTypeObject}, unsafe_load(o).ob_type)

PyObject(t::PyTypeObject) = pyincref(convert(PyPtr, pointer_from_objref(t)))


text:
cd(@__DIR__)
allowed = ["case", "space", "flux", "collision"]
D = KitBase.read_dict("config.txt", allowed)

D = KitBase.read_dict("config.txt")
ks, ctr, face, simTime = KitBase.initialize(D)

ks, ctr, face, simTime = KitBase.initialize("config.txt")
KitBase.write_jld(ks, ctr)

using KitBase.Plots
plot_line(ks, ctr)
#plot(ks, ctr, legend=:none)

ext: jl


text:
function generate_object_defs(domain::Domain, state::State,
                              domain_type::Symbol, state_type::Symbol)
    objects = sort(collect(get_objects(state)), by=x->x.name)
    object_ids = (; ((o.name, i) for (i, o) in enumerate(objects))...)
    get_objects_def =
        :(get_objects(::$state_type) = $(QuoteNode(Tuple(objects))))
    get_objtypes_def =
        :(get_objtypes(::$state_type) = $(QuoteNode(get_objtypes(state))))
    objectindices_def =
        :(objectindices(::$state_type) = $(QuoteNode(object_ids)))
    objectindex_def =
        :(objectindex(state::$state_type, o::Symbol) =
            getfield(objectindices(state), o))
    typed_defs = !get_requirements(domain)[:typing] ? Expr(:block) :
        generate_object_typed_defs(domain, state, domain_type, state_type)
    return Expr(:block, get_objects_def, get_objtypes_def,
                objectindices_def, objectindex_def, typed_defs)
end

function generate_object_typed_defs(domain::Domain, state::State,
                                    domain_type::Symbol, state_type::Symbol)
    object_ids = Dict()
    for type in get_types(domain)
        objs = sort(get_objects(domain, state, type), by=x->x.name)
        object_ids[type] = (; ((o.name, i) for (i, o) in enumerate(objs))...)
    end
    object_ids = (; (ty => ids for (ty, ids) in object_ids)...)
    get_objects_def =
        :(get_objects(state::$state_type, type::Symbol) =
            Const.(keys(objectindices(state, type))))
    objectindices_def =
        :(objectindices(state::$state_type, type::Symbol) =
            getfield($(QuoteNode(object_ids)), type))
    objectindex_def =
        :(objectindex(state::$state_type, type::Symbol, obj::Symbol) =
            getfield(objectindices(state, type), obj))
    return Expr(:block, get_objects_def, objectindices_def, objectindex_def)
end


text:
module HidApi

import hidapi_jll
const hidapi = hidapi_jll.hidapi

include(joinpath(@__DIR__, "hidapi_common.jl"))
include(joinpath(@__DIR__, "hidapi_api.jl"))
include(joinpath(@__DIR__, "high_level_api.jl"))

export HidDevice
export init, shutdown, open, close, read, write
export find_device, enumerate_devices

end # module


text:
"""
    NeighborSampler{S<:AbstractGraph, T} <: Sampleable{Univariate, Discrete}

A sampler that makes it easy to sample `vertex`'s neighbors in `graph`.

##### Fields
- `graph::S` : a graph.
- `vertex::T` : a vertex that is in `graph`'s vertex set.
"""
struct NeighborSampler{S<:AbstractGraph, T} <: Sampleable{Univariate, Discrete}
    graph::S
    vertex::T
end

"""
    NeighborSampler(graph::AbstractGraph, vertex::Int; checkargs::Bool = true) -> (
                    NeighborSampler)

Construct a `NeighborSampler`, optionally checking to ensure that the arguments are valid.

##### Complexity
- Time complexity: `O(n)`, when `checkargs == true`; otherwise, `O(1)`.
"""
function NeighborSampler(graph::AbstractGraph, vertex::Int; checkargs::Bool = true)
    # don't check arguments
    !checkargs && return NeighborSampler{typeof(graph)}(graph, vertex)
    # check arguments
    if !(vertex in vertices(graph))
        throw(ArgumentError("The vertex must be in the graph's vertex set."))
    else
        return NeighborSampler{typeof(graph), typeof(vertex)}(graph, vertex)
    end
end

"""
    Base.show(io::IO, neighborsampler::NeighborSampler) -> nothing

Show an `NeighborSampler` in a human-friendly manner.
"""
function Base.show(io::IO, neighborsampler::NeighborSampler)
    println(io, typeof(neighborsampler))
    println(io, " graph: ", neighborsampler.graph)
    println(io, " vertex: ", neighborsampler.vertex)
    return nothing
end

"""
    Base.rand(rng::AbstractRNG, neighborsampler::NeighborSampler) -> Int

Select and return one of `vertex`'s neighbors uniformly-at-random.
"""
function Base.rand(rng::AbstractRNG, neighborsampler::NeighborSampler)
    vertneighbors = neighbors(neighborsampler.graph, neighborsampler.vertex)
    if !isempty(vertneighbors)
        return rand(vertneighbors)
    # if vertex has no neighbors, then throw an error
    else
        error("vertex has no neighbors.")
    end
end

max_stars_repo_name: LoganDGraham/StochasticProcesses

max_issues_repo_name: LoganDGraham/StochasticProcesses

max_forks_repo_name: LoganDGraham/StochasticProcesses


text:
include("exactprimal.jl")
include("subgradient.jl")
include("cuttingdual.jl")

abstract type RegressionMethod end
abstract type WarmStartMethod <: RegressionMethod end
struct NoWarmStart <: WarmStartMethod end

"""
    ExactPrimalCuttingPlane

Use MIP solver to solve for optimal s.
"""
mutable struct ExactPrimalCuttingPlane <: RegressionMethod
    gamma::Float64
    time_limit::Float64
    solver::MathProgBase.AbstractMathProgSolver
    warm_start::WarmStartMethod
end
ExactPrimalCuttingPlane() = ExactPrimalCuttingPlane(0.0, 30.0, UnsetSolver, NoWarmStart())
"""
    ExactPrimalCuttingPlane(gamma::Float64, time_limit::Float64,
            stype::Type{S}, warm_start::WarmStartMethod) where {S <: MathProgBase.AbstractMathProgSolver}
"""
function ExactPrimalCuttingPlane(gamma::Float64, time_limit::Float64, stype::Type{S}, warm_start::WarmStartMethod) where {S <: MathProgBase.AbstractMathProgSolver}
      solver = getsolver(stype, time_limit)
      ExactPrimalCuttingPlane(gamma, time_limit, solver, warm_start)
end
function ExactPrimalCuttingPlane(stype::Type{S}) where {S <: MathProgBase.AbstractMathProgSolver}
      ExactPrimalCuttingPlane(0.0, 30.0, stype, NoWarmStart())
end

"""
    PrimalWithHeuristics

Use MIP solver and supply it with node heuristics.
"""
mutable struct PrimalWithHeuristics <: RegressionMethod
    gamma::Float64
    time_limit::Float64
    solver::MathProgBase.AbstractMathProgSolver
    warm_start::WarmStartMethod
end
PrimalWithHeuristics() = PrimalWithHeuristics(0.0, 30.0, UnsetSolver, NoWarmStart())
"""
    PrimalWithHeuristics(gamma::Float64, time_limit::Float64,
            stype::Type{S}, warm_start::WarmStartMethod) where {S <: MathProgBase.AbstractMathProgSolver}
"""
function PrimalWithHeuristics(gamma::Float64, time_limit::Float64, stype::Type{S}, warm_start::WarmStartMethod) where {S <: MathProgBase.AbstractMathProgSolver}
      solver = getsolver(stype, time_limit)
      ExactPrimalCuttingPlane(gamma, time_limit, solver, warm_start)
end
function PrimalWithHeuristics(stype::Type{S}) where {S <: MathProgBase.AbstractMathProgSolver}
      ExactPrimalCuttingPlane(0.0, 30.0, stype,  NoWarmStart())
end

"""
    RelaxDualSubgradient

Use subgradient descent on the dual of the convex relaxation.
"""
mutable struct RelaxDualSubgradient{SR <: SteppingRule} <: WarmStartMethod
    gamma::Float64
    sr::SR
    maxiter::Int
end
RelaxDualSubgradient() = RelaxDualSubgradient(0.0, ConstantStepping(1e-3), 100)

"""
    RelaxDualCutting

Use Kelley's cuts on dual of convex relaxation.
"""
mutable struct RelaxDualCutting <: WarmStartMethod
    gamma::Float64
    solver::MathProgBase.AbstractMathProgSolver
end
RelaxDualCutting() = RelaxDualCutting(0.0, UnsetSolver())

function solve_problem(::RegressionMethod, ::Array{Float64,2}, ::YVector, ::Int)
    error("You need to define `solve_problem` for $m.")
end

function solve_problem(m::ExactPrimalCuttingPlane, X::Array{Float64,2}, Y::YVector, sparsity::Int)
    if m.warm_start != NoWarmStart()
        inds, _ = solve_problem(m.warm_start, X, Y, sparsity)
        indices0, w0, Δt, status, Gap, cutCount = oa_formulation_bm(SubsetSelection.OLS(), Y, X, sparsity, m.gamma, m.solver, indices0 = inds)
    else
        indices0, w0, Δt, status, Gap, cutCount = oa_formulation_bm(SubsetSelection.OLS(), Y, X, sparsity, m.gamma, m.solver)
    end
    indices0, w0
end

function solve_problem(m::PrimalWithHeuristics, X::Array{Float64,2}, Y::YVector, sparsity::Int)
    if m.warm_start != NoWarmStart()
        inds, _ = solve_problem(m.warm_start, X, Y, sparsity)
        indices0, w0, Δt, status, Gap, cutCount = oa_formulation_bm(SubsetSelection.OLS(), Y, X, sparsity, m.gamma, m.solver, node_heuristics=true, indices0=inds)
    else
        indices0, w0, Δt, status, Gap, cutCount = oa_formulation_bm(SubsetSelection.OLS(), Y, X, sparsity, m.gamma, m.solver, node_heuristics=true)
    end
    indices0, w0
end

function solve_problem(m::RelaxDualSubgradient, X::Array{Float64,2}, Y::YVector, k::Int)
    sparsity = SubsetSelection.Constraint(k)
    sp = subsetSelection_bm(SubsetSelection.OLS(), sparsity, Y, X, γ = m.gamma, sr = m.sr, maxIter = m.maxiter)
    sp.indices, sp.w
end

function solve_problem(m::RelaxDualCutting, X::Array{Float64,2}, Y::YVector, k::Int)
    solve_dualcutting(X, Y, k, m.gamma, m.solver)
end


text:
mutable struct Tariff
    provider::String
    category::String
    e_cost::Array{Any}
    p_cost::Array{Any}
    increase::Float64
    access::Float64
    street_light::Float64
    access_increase::Float64
    behaviour::Array{Float64}
	min_cost::Float64
end

max_stars_repo_name: arnavgautam/PVCR

max_issues_repo_name: arnavgautam/PVCR

max_forks_repo_name: arnavgautam/PVCR


text:
# utilitaries
using LinearAlgebra

include("quadrature.jl")
# include("connectique.jl")

@inline function mod_wrap(x, a)
    #= return the index x modulo a =#
    return mod(x-1, a) + 1
end

function polygonArea(vtx)
    area = 0.
    N = size(vtx, 2)
    for i=1:N
        j = mod_wrap(i+1, N)
        area += vtx[1, i]*vtx[2, j] - vtx[2, i]*vtx[1, j]
    end
    return 0.5*abs(area)
end

function polygonDiameter(vtx)
    diameter = 0.
    N = size(vtx, 2)
    for i=1:(N-1)
        for j=(i+1):N
            diameter = max(diameter, norm(vtx[:, i] - vtx[:, j]))
        end
    end
    return diameter
end

function polygonCentroid(vtx)
    N        = size(vtx, 2)
    centroid = zeros(Float64, 2)
    #
    den = 0.
    for i=1:N
        j = mod_wrap(i+1, N)
        w = (vtx[1, i]*vtx[2, j] - vtx[1, j]*vtx[2, i])
        centroid[1] += (vtx[1, i] + vtx[1, j])*w
        centroid[2] += (vtx[2, i] + vtx[2, j])*w
        den += w
    end
    #
    return centroid/(3 * den)
end

# INTEGRATION OF MULTIPLE FUNCTIONS OVER A POLYGON
function integratePolygon(pverts::Array{Float64,2}, xE::Array{Float64,1},
    ref::Quadrature2d, fs::Array{Function, 1})
    #=
    Integration of all functions in fs over a polygon by subdividing it
    into triangles
    =#
    n   = size(pverts, 2)
    nf  = length(fs)
    res = zeros(Float64, nf)
    for i=1:n
        ip = mod_wrap(i+1, n)
        @inbounds quadT = quadGL2d(ref, [pverts[1,i] pverts[1,ip] xE[1];
        pverts[2,i] pverts[2,ip] xE[2]])
        for f = 1:nf
            @inbounds res[f] += integrate(quadT, fs[f])
        end
    end
    return  res
end

function integratePolygon(pverts::Array{Float64,2}, xE::Array{Float64,1},
    ref::Quadrature2d, fs::Function)
    #= Same as before, but for a single function fs =#
    n   = size(pverts, 2)
    res = 0.
    for i=1:n
        ip = mod_wrap(i+1, n)
        @inbounds quadT = quadGL2d(ref, [pverts[1,i] pverts[1,ip] xE[1];
        pverts[2,i] pverts[2,ip] xE[2]])
        res += integrate(quadT, fs)
    end
    return  res
end


text:
import NLPModelsIpopt: ipopt

"""
ipopt(nlp) DOESN'T CHECK THE WRONG KWARGS, AND RETURN AN ERROR.

ipopt(::NLPStopping)

"""
function NLPModelsIpopt.ipopt(stp :: NLPStopping; kwargs...)

  #xk = solveIpopt(stop.pb, stop.current_state.x)
  nlp = stp.pb
  stats = ipopt(nlp, print_level     = 0,
                      tol             = stp.meta.rtol,
                      x0              = stp.current_state.x,
                      max_iter        = stp.meta.max_iter,
                      max_cpu_time    = stp.meta.max_time,
                      dual_inf_tol    = stp.meta.atol,
                      constr_viol_tol = stp.meta.atol,
                      compl_inf_tol   = stp.meta.atol,
                      kwargs...)

  #Update the meta boolean with the output message
  #=
  if stats.status == :first_order stp.meta.suboptimal      = true end
  if stats.status == :acceptable  stp.meta.suboptimal      = true end
  if stats.status == :infeasible  stp.meta.infeasible      = true end
  if stats.status == :small_step  stp.meta.stalled         = true end
  if stats.status == :max_eval    stp.meta.max_eval        = true end
  if stats.status == :max_iter    stp.meta.iteration_limit = true end
  if stats.status == :max_time    stp.meta.tired           = true end
  if stats.status ∈ [:neg_pred,
                     :not_desc]   stp.meta.fail_sub_pb     = true end
  if stats.status == :unbounded   stp.meta.unbounded       = true end
  if stats.status == :user        stp.meta.stopbyuser      = true end
  if stats.status ∈ [:stalled,
                     :small_residual,
                     :small_step]   stp.meta.stalled       = true end
  #if stats.status == :exception   stp.meta.exception       = true end #available ≥ 0.2.6
  =#
  stp = stats_status_to_meta!(stp, stats)

  if status(stp) == :Unknown
    @warn "Error in StoppingInterface statuses: return status is $(stats.status)"
  end

  stp.meta.nb_of_stop = stats.iter
  #stats.elapsed_time

  x = stats.solution

  #Not mandatory, but in case some entries of the State are used to stop
  fill_in!(stp, x) #too slow

  stop!(stp)

  return stp
end


text:
### A Pluto.jl notebook ###
# v0.14.4

using Markdown
using InteractiveUtils

# ╔═╡ c9cf879e-0b01-4e0a-b31c-2ca5c6f03e63
begin
	using PlutoUI
	PlutoUI.TableOfContents(aside=true)
end

# ╔═╡ 94dde840-54f6-4fe4-846b-bc15f8fca9fa
begin
	# Packages
	using LightGraphs
	using GraphPlot
	using Clustering
	using SparseArrays
	using LinearAlgebra
	using Plots
	using Distances
	using Random
end

# ╔═╡ 81efb137-7e14-4091-a7f0-c6b72e2d77f7
md"""
# Spectral Graph K-partitioning

Instead of using recursive spectral bipartitioning, the graph $k$-partitioning problem can be solved using $k$ eigenvectors which correspond to $k$ smallest eigenvalues of Laplace matrix (Laplacian) or normalized Laplace matrix (normalized Laplacian), respectively.

Suggested reading is [U. von Luxburg, A Tutorial on Spectral Clustering](http://www.kyb.mpg.de/fileadmin/user_upload/files/publications/attachments/luxburg06_TR_v2_4139%5b1%5d.pdf), which includes the quote

_"Spectral clustering cannot serve as a “black box algorithm” which automatically detects the correct clusters in any given data set. But it can be considered as a powerful tool which can produce good results if applied with care."_

__Prerequisites__

The reader should be familiar with k-means algorithm, spectral graph bipartitioning and recursive bipartitioning.

__Competences__

The reader should be able to apply graph spectral k-partitioning to data clustering problems.

__Credits.__ The notebook was initially derived from M.Sc. Thesis of Ivančica Mirošević.
"""

# ╔═╡ 75cd6187-4c51-426c-8a3c-d7b6c922ae0b
md"""

## Definitions

Let $G=(V,E)$ be a weighted graph with weights $\omega$, with weights matrix $W$, Laplacian matrix $L=D-W$, and normalized Laplacian matrix $L_n=D^{-1/2}(D-W)D^{-1/2}$.

Let the $k$-partition $\pi_k =\{V_{1},V_{2},...,V_{k}\}$, the cut $\mathop{\mathrm{cut}}(\pi_k)$, the proportional cut $\mathop{\mathrm{pcut}}(\pi_k)$ and the normalized cut $\mathop{\mathrm{ncut}}(\pi_k)$ be defined as in the
[Spectral Graph Bipartitioning](https://ivanslapnicar.github.io/NumericalMathematics/C2%20Spectral%20Graph%20Bipartitioning.jl.html) notebook.

__Partition vectors__ of a $k$-partition $\pi_k$ are

$$
\begin{aligned}
h_{1} &=[\overset{\displaystyle |V_{1}|}{\overbrace{1,\cdots
,1}},0,\cdots
,0,\cdots ,0,\cdots ,0]^{T}  \\
h_{2} &=[0,\cdots
,0,\overset{\displaystyle |V_{2}|}{\overbrace{1,\cdots ,1}} ,\cdots ,0,\cdots
,0]^{T} \\
&\vdots \\
h_{k} &=[0,\cdots ,0,0,\cdots ,0,\cdots
,\overset{\displaystyle |V_{k}|}{ \overbrace{1,\cdots ,1}}]^{T}.
\end{aligned}$$

## Facts

1. Set

$$
\begin{aligned}
X&=\begin{bmatrix} x_1 & x_2 &\cdots & x_k \end{bmatrix}, \quad
x_i=\displaystyle \frac{h_i}{\|h_i\|_2}, \\
Y&=\begin{bmatrix} y_1 & y_2 &\cdots & y_k \end{bmatrix}, \quad
y_i=\displaystyle \frac{D^{1/2}h_i}{\|D^{1/2}h_i\|_2}.
\end{aligned}$$

It holds

$$
\begin{aligned}
& \mathop{\mathrm{cut}}(V_{i},V\backslash V_{i})=h_{i}^{T}(D-W)h_{i}=h_{i}^{T}L h_{i},\quad
\omega( C_{i})=h_{i}^{T}D h_{i},\quad |C_{i}| =h_{i}^{T}h_{i},\\
& \mathop{\mathrm{pcut}}(\pi_k) =\frac{h_{1}^{T}L h_{1}}{h_{1}^{T} h_{1}}+
\cdots + \frac{h_{k}^{T}L h_{k}}{h_{k}^{T}h_{k}}
=x_{1}^{T}L x_{1}+\cdots +x_{k}^{T}Lx_{k}=\mathop{\mathrm{trace}}(X^{T}LX),\\
& \mathop{\mathrm{ncut}}(\pi_k)=\frac{h_{1}^{T}L h_{1}}{h_{1}^{T}D h_{1}}+\cdots
+\frac{h_{k}^{T}L h_{k}}{h_{k}^{T}D h_{k}}
=\mathop{\mathrm{trace}}(Y^{T}L_{n}Y).
\end{aligned}$$

2. The __relaxed__ $k$-partitioning problems are trace-minimization problems,

$$
\begin{aligned}
\min_{\displaystyle \pi_k} \mathop{\mathrm{pcut}}(\pi_k) &\geq
\underset{\displaystyle X\in\mathbb{R}^{n\times k}}{\min_{\displaystyle X^{T}X=I}}
\mathop{\mathrm{trace}}(X^{T}LX),\\
\min_{\displaystyle \pi_k} \mathop{\mathrm{ncut}}(\pi_k) &\geq
\underset{\displaystyle Y\in\mathbb{R}^{n\times
k}}{\min_{\displaystyle Y^{T}Y=I}}\mathop{\mathrm{trace}}(Y^{T}L_{n}Y).
\end{aligned}$$

3. __Ky-Fan  Theorem.__ Let $A\in \mathbb{R}^{n\times n}$ be a symmetric matrix with eigenvalues $\lambda _1\leq \cdots \leq \lambda_n$. Then

$$
\underset{\displaystyle Z^TZ=I}{\min_{\displaystyle Z\in \mathbb{R}^{n\times
k}}}\mathop{\mathrm{trace}}\left(Z^{T}AZ\right)
=\sum_{i=1}^{k}\lambda_{i}.$$

4. Let $\lambda_1\leq \cdots \leq \lambda_n$ be the eigenvalues of $L$ with eigenvectors $v^{[1]},\cdots ,v^{[k]}$. The solution of the relaxed proportional cut problem is the matrix $X=\begin{bmatrix}v^{[1]} & \cdots & v^{[k]}\end{bmatrix}$, and it holds

$\min\limits_{\displaystyle \pi_k} \mathop{\mathrm{pcut}}(\pi_k)\geq \sum\limits_{i=1}^k \lambda_i.$

5. Let $\mu_1\leq \cdots \leq \mu_n$ be the eigenvalues of $L_n$ with eigenvectors $w^{[1]},\cdots ,w^{[k]}$. The solution of the relaxed normalized cut problem is the matrix $Y=\begin{bmatrix}w^{[1]} & \cdots & w^{[k]}\end{bmatrix}$, and it holds

$\min\limits_{\displaystyle \pi_k} \mathop{\mathrm{ncut}}(\pi_k)\geq \sum\limits_{i=1}^k \mu_i$.

6. It remains to recover the $k$-partition. The k-means algorithm applied to rows of the matrices $X$ or $D^{-1/2}Y$, will compute the $k$ centers and the assignment vector whose $i$-th component denotes the subset $V_j$ to which the vertex $i$ belongs.
"""

# ╔═╡ 639758a5-413d-408a-ace9-3741322c3dc2
md"""

## Examples

### Graph with three clusters
"""

# ╔═╡ d367f447-16e1-4cbc-b05f-b807d4ef955b
begin
	# Sources, targets, and weights
	n=9
	sn=[1,1,1,2,2,2,3,3,5,6,7,7,8]
	tn=[2,3,4,3,4,7,4,5,6,9,8,9,9]
	wn=[2,3,4,4,5,1,6,1,7,1,4,3,2]
	[sn tn wn]
end

# ╔═╡ 9438ace8-4290-4caa-bb34-24c1c8f4f66a
begin
	# What is the optimal tripartition?
	G=Graph(n)
	for i=1:length(sn)
	    add_edge!(G,sn[i],tn[i])
	end
	gplot(G, nodelabel=1:n, edgelabel=wn)
end

# ╔═╡ f05bae5c-8d1c-4f68-a299-fab3af2432ae
begin
	# We define some functions
	function WeightMatrix(src::Array,dst::Array,weights::Array)
	    n=nv(G)
	    sparse([src;dst],[dst;src],[weights;weights],n,n)
	end

	Laplacian(W::AbstractMatrix)=spdiagm(0=>vec(sum(W,dims=2)))-W

	function NormalizedLaplacian(L::AbstractMatrix)
	    D=1.0./sqrt.(diag(L))
	    n=length(D)
	    [L[i,j]*(D[i]*D[j]) for i=1:n, j=1:n]
	end
end

# ╔═╡ 357b8af9-700a-4a2d-9703-575af5b9ab13
begin
	W=WeightMatrix(sn,tn,wn)
	L=Laplacian(W)
	Lₙ=NormalizedLaplacian(L)
end

# ╔═╡ 085b49b4-08a2-4e6e-8cbe-6a011a93b551
Matrix(L)

# ╔═╡ 280049be-261e-49d7-a6fc-2211afe94e06
# Proportional cut. The clustering is visible in
# the components of v₂ and v₃
# E=eigs(L,nev=3,which=:SM)
E=eigen(Matrix(L))

# ╔═╡ 64730ce5-61b4-405d-9a41-81d0253a4713
# Check the assignments
out=kmeans(Matrix(transpose(E.vectors[:,1:3])),3)

# ╔═╡ 692cbb97-2a5f-42c9-b829-e08d5e7479a8
begin
	# Normalized cut
	# Lanczos cannot be used for the "smallest in magnitude"
	# eienvalues of a singular matrix
	# λ,Y=eigs(Ln,nev=3,which=:SM) does not work
	Eₙ=eigen(Lₙ)
	D=sqrt.(diag(L))
	Y=inv(Diagonal(D))*Eₙ.vectors[:,1:3]
	outₙ=kmeans(Matrix(transpose(Y)),3)
end

# ╔═╡ 662b7b01-883f-4d51-b696-26f7559e0211
md"""
### Concentric rings
"""

# ╔═╡ 55d217d5-c31a-431d-bd46-f9f09280cf5b
function plotKpartresult(C::Vector,X::Array)
    k=maximum(C)
    scatter(aspect_ratio=1)
    for j=1:k
        scatter!(X[1,findall(C.==j)],X[2,findall(C.==j)],label="Cluster $j")
    end
	scatter!()
end

# ╔═╡ 3f815c3e-850d-4ca7-af1c-a2ff14e429db
begin
	# Generate concentric rings
	k=4
	Random.seed!(541)
	# Center
	center=[0,0]
	# Radii
	radii=Random.randperm(10)[1:k]
	# Number of points in rings
	sizes=rand(300:500,k)
	center,radii,sizes
end

# ╔═╡ 967b8e40-1f95-11eb-0aa1-8f4592054357
begin
	# Generate points
	X=Array{Float64}(undef,2,sum(sizes))
	csizes=cumsum(sizes)
	# Random angles
	ϕ=2*π*rand(sum(sizes))
	for i=1:csizes[1]
		X[:,i]=center+radii[1]*[cos(ϕ[i]);sin(ϕ[i])] + (rand(2).-0.5)/50
	end
	for j=2:k
		for i=csizes[j-1]+1:csizes[j]
			X[:,i]=center+radii[j]*[cos(ϕ[i]);sin(ϕ[i])] + (rand(2).-0.5)/50
		end
	end
	scatter(X[1,:],X[2,:],title="Concentric rings", aspect_ratio=1,label="Points")
end

# ╔═╡ 708a16c5-f755-490b-abcb-89788813b7d0
begin
	S=pairwise(SqEuclidean(),X,dims=2)
	# S=pairwise(Cityblock(),X)
	β=60
	W₁=exp.(-β*S);
end

# ╔═╡ c298ccfe-f15e-4652-baa7-c122cb22f22f
begin
	L₁=Laplacian(W₁)
	Ln₁=NormalizedLaplacian(L₁);
end

# ╔═╡ 9d8aef37-9640-4d67-9c33-bb6513938ae1
begin
	# Laplacian
	E₁=eigen(L₁)
	sp=sortperm(abs.(E₁.values))[1:k]
	λ₁=E₁.values[sp]
	Y₁=E₁.vectors[:,sp]
	out₁=kmeans(Matrix(transpose(Y₁)),k)
	plotKpartresult(out₁.assignments,X)
end

# ╔═╡ 7c72773d-0915-4ef0-aa7e-7f0474426bfa
begin
	# Normalized Laplacian
	En₁=eigen(Ln₁)
	spn=sortperm(abs.(En₁.values))[1:k]
	λ=En₁.values[spn]
	Yn₁=En₁.vectors[:,spn]
	Yn₁=Diagonal(1.0./sqrt.(diag(L₁)))*Y₁
	outn₁=kmeans(Matrix(transpose(Y₁)),k)
	plotKpartresult(outn₁.assignments,X)
end

# ╔═╡ 68d822d0-1f97-11eb-1a79-8f58785ae13a
L₁[3,1]

# ╔═╡ Cell order:
# ╟─c9cf879e-0b01-4e0a-b31c-2ca5c6f03e63
# ╟─81efb137-7e14-4091-a7f0-c6b72e2d77f7
# ╟─75cd6187-4c51-426c-8a3c-d7b6c922ae0b
# ╟─639758a5-413d-408a-ace9-3741322c3dc2
# ╠═94dde840-54f6-4fe4-846b-bc15f8fca9fa
# ╠═f05bae5c-8d1c-4f68-a299-fab3af2432ae
# ╠═d367f447-16e1-4cbc-b05f-b807d4ef955b
# ╠═9438ace8-4290-4caa-bb34-24c1c8f4f66a
# ╠═357b8af9-700a-4a2d-9703-575af5b9ab13
# ╠═085b49b4-08a2-4e6e-8cbe-6a011a93b551
# ╠═280049be-261e-49d7-a6fc-2211afe94e06
# ╠═64730ce5-61b4-405d-9a41-81d0253a4713
# ╠═692cbb97-2a5f-42c9-b829-e08d5e7479a8
# ╟─662b7b01-883f-4d51-b696-26f7559e0211
# ╠═55d217d5-c31a-431d-bd46-f9f09280cf5b
# ╠═3f815c3e-850d-4ca7-af1c-a2ff14e429db
# ╠═967b8e40-1f95-11eb-0aa1-8f4592054357
# ╠═708a16c5-f755-490b-abcb-89788813b7d0
# ╠═c298ccfe-f15e-4652-baa7-c122cb22f22f
# ╠═7c72773d-0915-4ef0-aa7e-7f0474426bfa
# ╠═9d8aef37-9640-4d67-9c33-bb6513938ae1
# ╠═68d822d0-1f97-11eb-1a79-8f58785ae13a

max_stars_repo_path: src/Module C - Applications/L11 Spectral Graph K-partitioning.jl

max_stars_repo_name: ivanslapnicar/GIAN-Applied-NLA-Course

max_issues_repo_name: ivanslapnicar/GIAN-Applied-NLA-Course

max_forks_repo_name: ivanslapnicar/GIAN-Applied-NLA-Course


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    CopyMapping

A mapping method in which data points are copied directly to the
domain at specified locations.
"""
struct CopyMapping{V1,V2} <: MappingMethod
  orig::V1
  dest::V2
end

CopyMapping(dest) = CopyMapping(nothing, dest)
CopyMapping() = CopyMapping(nothing, nothing)

function map(sdata, sdomain, targetvars, method::CopyMapping)
  @assert targetvars ⊆ name.(variables(sdata)) "target variables must be present in spatial data"

  # dictionary with mappings
  mappings = Dict(var => Dict{Int,Int}() for var in targetvars)

  # retrieve origin and destination indices
  orig = isnothing(method.orig) ? (1:nelms(sdata)) : method.orig
  dest = isnothing(method.dest) ? (1:nelms(sdata)) : method.dest

  @assert length(orig) == length(dest) "invalid mapping specification"

  for i in eachindex(orig, dest)
    # save pair if there is data for variable
    for var in targetvars
      if !ismissing(sdata[var][orig[i]])
        push!(mappings[var], dest[i] => orig[i])
      end
    end
  end

  mappings
end


text:
using Printf

println("Add, subtract, multiply, divide")
@printf("1 + 2 = %d\n", 1 + 2)

@printf("1 - 2 = %d\n", 1 - 2)
@printf("3 * 4 = %d\n", 3 * 4)
@printf("10 / 3 = %.3f\n", 10 / 3)

println("Other operators")
@printf("6^3 = %d\n", 6^3)
@printf("2^0.5 = %f\n", 2^0.5)
@printf("div(11, 4) = %d\n", div(11, 4))
@printf("rem(11, 4) = %d\n", rem(11, 4))
println("divrem(11, 4) = ", divrem(11, 4))
println("7 % 3 = ", 7 % 3)

println("Types and rational numbers")
@printf("typeof(5): %s\n", typeof(5))
@printf("typeof(1.23): %s\n", typeof(1.23))
a = rationalize(5.6)
println("a = ", a)
@printf("typeof(a): %s\n", typeof(a))
println(rationalize(BigInt, 10.3))
println("2 // 3 + 4 // 5 = ", 2//3 + 4//5)
@printf("numerator(a) = %d\n", numerator(a))
@printf("denominator(a) = %d\n", denominator(a))

println("Equality testing")
println(1 == 2)
println(1 + 1 == 2)
println(1 != 2)

println("Inequality testing")
println(1 < 2)
println(1 > 2)
println(1 <= 2)
println(1 <= 1)
println(1 > 2)
println(2 >= 2)

println("Identity testing")
println(1 === 2)
println(1 + 1 === 2)
println(1 !== 2)

println("Type testing")
println(isa(1, Int))
println(isa(1, Float64))
println(isa(1, Char))
println(isa(1, Matrix))

println("Some extra stuff")
println(1_000_000) # An integer
println((1, 000, 000)) # Tuple of integers
println(Inf, "\n", -Inf)
println(Inf + -Inf)
println(Inf * -Inf)
println(Inf + Inf)
println(pi)
println(round(1234, sigdigits=2))
println(round(123.456, digits=2))


text:
const CGDisplayNoErr = 0
const CGEventNoErr = 0
const CGFLOAT_DEFINED = 1
const CGFLOAT_IS_DOUBLE = 0
const CGFLOAT_MAX = 3.4028235E+38
const CGFLOAT_MIN = 1.1754944E-38
const CGGlyphMax = 1
const CGGlyphMin = 0
const CGPDFDataFormatJPEG2000 = 2
const CGPDFDataFormatJPEGEncoded = 1
const CGPDFDataFormatRaw = 0
const CGRectMaxXEdge = 2
const CGRectMaxYEdge = 3
const CGRectMinXEdge = 0
const CGRectMinYEdge = 1
const CGVECTOR_DEFINED = 1
const kCGAnnotatedSessionEventTap = 2
const kCGAnyInputEventType = 4294967295
const kCGAssistiveTechHighWindowLevel = 1500
const kCGAssistiveTechHighWindowLevelKey = 20
const kCGBackingStoreBuffered = 2
const kCGBackingStoreNonretained = 1
const kCGBackingStoreRetained = 0
const kCGBackstopMenuLevel = -20
const kCGBackstopMenuLevelKey = 3
const kCGBaseWindowLevel = -2147483648
const kCGBaseWindowLevelKey = 0
const kCGBitmapAlphaInfoMask = 31
const kCGBitmapByteOrder16Big = 12288
const kCGBitmapByteOrder16Host = 4096
const kCGBitmapByteOrder16Little = 4096
const kCGBitmapByteOrder32Big = 16384
const kCGBitmapByteOrder32Host = 8192
const kCGBitmapByteOrder32Little = 8192
const kCGBitmapByteOrderDefault = 0
const kCGBitmapByteOrderMask = 28672
const kCGBitmapFloatComponents = 256
const kCGBitmapFloatInfoMask = 3840
const kCGBlendModeClear = 16
const kCGBlendModeColor = 14
const kCGBlendModeColorBurn = 7
const kCGBlendModeColorDodge = 6
const kCGBlendModeCopy = 17
const kCGBlendModeDarken = 4
const kCGBlendModeDestinationAtop = 24
const kCGBlendModeDestinationIn = 22
const kCGBlendModeDestinationOut = 23
const kCGBlendModeDestinationOver = 21
const kCGBlendModeDifference = 10
const kCGBlendModeExclusion = 11
const kCGBlendModeHardLight = 9
const kCGBlendModeHue = 12
const kCGBlendModeLighten = 5
const kCGBlendModeLuminosity = 15
const kCGBlendModeMultiply = 1
const kCGBlendModeNormal = 0
const kCGBlendModeOverlay = 3
const kCGBlendModePlusDarker = 26
const kCGBlendModePlusLighter = 27
const kCGBlendModeSaturation = 13
const kCGBlendModeScreen = 2
const kCGBlendModeSoftLight = 8
const kCGBlendModeSourceAtop = 20
const kCGBlendModeSourceIn = 18
const kCGBlendModeSourceOut = 19
const kCGBlendModeXOR = 25
const kCGCaptureNoFill = 1
const kCGCaptureNoOptions = 0
const kCGColorConversionTransformApplySpace = 2
const kCGColorConversionTransformFromSpace = 0
const kCGColorConversionTransformToSpace = 1
const kCGColorSpaceModelCMYK = 2
const kCGColorSpaceModelDeviceN = 4
const kCGColorSpaceModelIndexed = 5
const kCGColorSpaceModelLab = 3
const kCGColorSpaceModelMonochrome = 0
const kCGColorSpaceModelPattern = 6
const kCGColorSpaceModelRGB = 1
const kCGColorSpaceModelUnknown = -1
const kCGConfigureForAppOnly = 0
const kCGConfigureForSession = 1
const kCGConfigurePermanently = 2
const kCGCursorWindowLevel = 2147483630
const kCGCursorWindowLevelKey = 19
const kCGDesktopIconWindowLevel = -2147483603
const kCGDesktopIconWindowLevelKey = 18
const kCGDesktopWindowLevel = -2147483623
const kCGDesktopWindowLevelKey = 2
const kCGDirectMainDisplay = 0
const kCGDisplayAddFlag = 16
const kCGDisplayBeginConfigurationFlag = 1
const kCGDisplayBlendNormal = 0.0
const kCGDisplayBlendSolidColor = 1.0
const kCGDisplayDesktopShapeChangedFlag = 4096
const kCGDisplayDisabledFlag = 512
const kCGDisplayEnabledFlag = 256
const kCGDisplayFadeReservationInvalidToken = 0
const kCGDisplayMirrorFlag = 1024
const kCGDisplayMovedFlag = 2
const kCGDisplayRemoveFlag = 32
const kCGDisplaySetMainFlag = 4
const kCGDisplaySetModeFlag = 8
const kCGDisplayStreamFrameStatusFrameBlank = 2
const kCGDisplayStreamFrameStatusFrameComplete = 0
const kCGDisplayStreamFrameStatusFrameIdle = 1
const kCGDisplayStreamFrameStatusStopped = 3
const kCGDisplayStreamUpdateDirtyRects = 2
const kCGDisplayStreamUpdateMovedRects = 1
const kCGDisplayStreamUpdateReducedDirtyRects = 3
const kCGDisplayStreamUpdateRefreshedRects = 0
const kCGDisplayUnMirrorFlag = 2048
const kCGDockWindowLevel = 20
const kCGDockWindowLevelKey = 7
const kCGDraggingWindowLevel = 500
const kCGDraggingWindowLevelKey = 12
const kCGEncodingFontSpecific = 0
const kCGEncodingMacRoman = 1
const kCGErrorCannotComplete = 1004
const kCGErrorFailure = 1000
const kCGErrorIllegalArgument = 1001
const kCGErrorInvalidConnection = 1002
const kCGErrorInvalidContext = 1003
const kCGErrorInvalidOperation = 1010
const kCGErrorNoneAvailable = 1011
const kCGErrorNotImplemented = 1006
const kCGErrorRangeCheck = 1007
const kCGErrorSuccess = 0
const kCGErrorTypeCheck = 1008
const kCGEventFilterMaskPermitAllEvents = 7
const kCGEventFilterMaskPermitLocalKeyboardEvents = 2
const kCGEventFilterMaskPermitLocalMouseEvents = 1
const kCGEventFilterMaskPermitSystemDefinedEvents = 4
const kCGEventFlagMaskAlphaShift = 65536
const kCGEventFlagMaskAlternate = 524288
const kCGEventFlagMaskCommand = 1048576
const kCGEventFlagMaskControl = 262144
const kCGEventFlagMaskHelp = 4194304
const kCGEventFlagMaskNonCoalesced = 256
const kCGEventFlagMaskNumericPad = 2097152
const kCGEventFlagMaskSecondaryFn = 8388608
const kCGEventFlagMaskShift = 131072
const kCGEventFlagsChanged = 12
const kCGEventKeyDown = 10
const kCGEventKeyUp = 11
const kCGEventLeftMouseDown = 1
const kCGEventLeftMouseDragged = 6
const kCGEventLeftMouseUp = 2
const kCGEventMaskForAllEvents = 18446744073709551615
const kCGEventMouseMoved = 5
const kCGEventMouseSubtypeDefault = 0
const kCGEventMouseSubtypeTabletPoint = 1
const kCGEventMouseSubtypeTabletProximity = 2
const kCGEventNull = 0
const kCGEventOtherMouseDown = 25
const kCGEventOtherMouseDragged = 27
const kCGEventOtherMouseUp = 26
const kCGEventRightMouseDown = 3
const kCGEventRightMouseDragged = 7
const kCGEventRightMouseUp = 4
const kCGEventScrollWheel = 22
const kCGEventSourceGroupID = 44
const kCGEventSourceStateCombinedSessionState = 0
const kCGEventSourceStateHIDSystemState = 1
const kCGEventSourceStateID = 45
const kCGEventSourceStatePrivate = -1
const kCGEventSourceUnixProcessID = 41
const kCGEventSourceUserData = 42
const kCGEventSourceUserID = 43
const kCGEventSuppressionStateRemoteMouseDrag = 1
const kCGEventSuppressionStateSuppressionInterval = 0
const kCGEventSupressionStateRemoteMouseDrag = 1
const kCGEventSupressionStateSupressionInterval = 0
const kCGEventTabletPointer = 23
const kCGEventTabletProximity = 24
const kCGEventTapDisabledByTimeout = 4294967294
const kCGEventTapDisabledByUserInput = 4294967295
const kCGEventTapOptionDefault = 0
const kCGEventTapOptionListenOnly = 1
const kCGEventTargetProcessSerialNumber = 39
const kCGEventTargetUnixProcessID = 40
const kCGFloatingWindowLevel = 3
const kCGFloatingWindowLevelKey = 5
const kCGFontIndexInvalid = 65535
const kCGFontIndexMax = 65534
const kCGFontPostScriptFormatType1 = 1
const kCGFontPostScriptFormatType3 = 3
const kCGFontPostScriptFormatType42 = 42
const kCGGesturePhaseBegan = 1
const kCGGesturePhaseCancelled = 8
const kCGGesturePhaseChanged = 2
const kCGGesturePhaseEnded = 4
const kCGGesturePhaseMayBegin = 128
const kCGGesturePhaseNone = 0
const kCGGlyphMax = 65534
const kCGGradientDrawsAfterEndLocation = 2
const kCGGradientDrawsBeforeStartLocation = 1
const kCGHIDEventTap = 0
const kCGHeadInsertEventTap = 0
const kCGHelpWindowLevel = 200
const kCGHelpWindowLevelKey = 16
const kCGImageAlphaFirst = 4
const kCGImageAlphaLast = 3
const kCGImageAlphaNone = 0
const kCGImageAlphaNoneSkipFirst = 6
const kCGImageAlphaNoneSkipLast = 5
const kCGImageAlphaOnly = 7
const kCGImageAlphaPremultipliedFirst = 2
const kCGImageAlphaPremultipliedLast = 1
const kCGImageByteOrder16Big = 12288
const kCGImageByteOrder16Little = 4096
const kCGImageByteOrder32Big = 16384
const kCGImageByteOrder32Little = 8192
const kCGImageByteOrderMask = 28672
const kCGInterpolationDefault = 0
const kCGInterpolationHigh = 3
const kCGInterpolationLow = 2
const kCGInterpolationMedium = 4
const kCGInterpolationNone = 1
const kCGKeyboardEventAutorepeat = 8
const kCGKeyboardEventKeyboardType = 10
const kCGKeyboardEventKeycode = 9
const kCGLineCapButt = 0
const kCGLineCapRound = 1
const kCGLineCapSquare = 2
const kCGLineJoinBevel = 2
const kCGLineJoinMiter = 0
const kCGLineJoinRound = 1
const kCGMainMenuWindowLevel = 24
const kCGMainMenuWindowLevelKey = 8
const kCGMaxDisplayReservationInterval = 15
const kCGMaximumWindowLevel = 2147483631
const kCGMaximumWindowLevelKey = 14
const kCGMinimumWindowLevel = -2147483643
const kCGMinimumWindowLevelKey = 1
const kCGModalPanelWindowLevel = 8
const kCGModalPanelWindowLevelKey = 10
const kCGMomentumScrollPhaseBegin = 1
const kCGMomentumScrollPhaseContinue = 2
const kCGMomentumScrollPhaseEnd = 3
const kCGMomentumScrollPhaseNone = 0
const kCGMouseButtonCenter = 2
const kCGMouseButtonLeft = 0
const kCGMouseButtonRight = 1
const kCGMouseDownEventMaskingDeadSwitchTimeout = 60.0
const kCGMouseEventButtonNumber = 3
const kCGMouseEventClickState = 1
const kCGMouseEventDeltaX = 4
const kCGMouseEventDeltaY = 5
const kCGMouseEventInstantMouser = 6
const kCGMouseEventNumber = 0
const kCGMouseEventPressure = 2
const kCGMouseEventSubtype = 7
const kCGMouseEventWindowUnderMousePointer = 91
const kCGMouseEventWindowUnderMousePointerThatCanHandleThisEvent = 92
const kCGNormalWindowLevel = 0
const kCGNormalWindowLevelKey = 4
const kCGNullDirectDisplay = 0
const kCGNullWindowID = 0
const kCGNumReservedWindowLevels = 16
const kCGNumberOfEventSuppressionStates = 2
const kCGNumberOfEventSupressionStates = 2
const kCGNumberOfWindowLevelKeys = 21
const kCGOverlayWindowLevel = 102
const kCGOverlayWindowLevelKey = 15
const kCGPDFArtBox = 4
const kCGPDFBleedBox = 2
const kCGPDFCropBox = 1
const kCGPDFMediaBox = 0
const kCGPDFObjectTypeArray = 7
const kCGPDFObjectTypeBoolean = 2
const kCGPDFObjectTypeDictionary = 8
const kCGPDFObjectTypeInteger = 3
const kCGPDFObjectTypeName = 5
const kCGPDFObjectTypeNull = 1
const kCGPDFObjectTypeReal = 4
const kCGPDFObjectTypeStream = 9
const kCGPDFObjectTypeString = 6
const kCGPDFTrimBox = 3
const kCGPathEOFill = 1
const kCGPathEOFillStroke = 4
const kCGPathElementAddCurveToPoint = 3
const kCGPathElementAddLineToPoint = 1
const kCGPathElementAddQuadCurveToPoint = 2
const kCGPathElementCloseSubpath = 4
const kCGPathElementMoveToPoint = 0
const kCGPathFill = 0
const kCGPathFillStroke = 3
const kCGPathStroke = 2
const kCGPatternTilingConstantSpacing = 2
const kCGPatternTilingConstantSpacingMinimalDistortion = 1
const kCGPatternTilingNoDistortion = 0
const kCGPopUpMenuWindowLevel = 101
const kCGPopUpMenuWindowLevelKey = 11
const kCGRenderingIntentAbsoluteColorimetric = 1
const kCGRenderingIntentDefault = 0
const kCGRenderingIntentPerceptual = 3
const kCGRenderingIntentRelativeColorimetric = 2
const kCGRenderingIntentSaturation = 4
const kCGScreenSaverWindowLevel = 1000
const kCGScreenSaverWindowLevelKey = 13
const kCGScreenUpdateOperationMove = 1
const kCGScreenUpdateOperationReducedDirtyRectangleCount = 2147483648
const kCGScreenUpdateOperationRefresh = 0
const kCGScrollEventUnitLine = 1
const kCGScrollEventUnitPixel = 0
const kCGScrollPhaseBegan = 1
const kCGScrollPhaseCancelled = 8
const kCGScrollPhaseChanged = 2
const kCGScrollPhaseEnded = 4
const kCGScrollPhaseMayBegin = 128
const kCGScrollWheelEventDeltaAxis1 = 11
const kCGScrollWheelEventDeltaAxis2 = 12
const kCGScrollWheelEventDeltaAxis3 = 13
const kCGScrollWheelEventFixedPtDeltaAxis1 = 93
const kCGScrollWheelEventFixedPtDeltaAxis2 = 94
const kCGScrollWheelEventFixedPtDeltaAxis3 = 95
const kCGScrollWheelEventInstantMouser = 14
const kCGScrollWheelEventIsContinuous = 88
const kCGScrollWheelEventMomentumPhase = 123
const kCGScrollWheelEventPointDeltaAxis1 = 96
const kCGScrollWheelEventPointDeltaAxis2 = 97
const kCGScrollWheelEventPointDeltaAxis3 = 98
const kCGScrollWheelEventScrollCount = 100
const kCGScrollWheelEventScrollPhase = 99
const kCGSessionEventTap = 1
const kCGStatusWindowLevel = 25
const kCGStatusWindowLevelKey = 9
const kCGTabletEventDeviceID = 24
const kCGTabletEventPointButtons = 18
const kCGTabletEventPointPressure = 19
const kCGTabletEventPointX = 15
const kCGTabletEventPointY = 16
const kCGTabletEventPointZ = 17
const kCGTabletEventRotation = 22
const kCGTabletEventTangentialPressure = 23
const kCGTabletEventTiltX = 20
const kCGTabletEventTiltY = 21
const kCGTabletEventVendor1 = 25
const kCGTabletEventVendor2 = 26
const kCGTabletEventVendor3 = 27
const kCGTabletProximityEventCapabilityMask = 36
const kCGTabletProximityEventDeviceID = 31
const kCGTabletProximityEventEnterProximity = 38
const kCGTabletProximityEventPointerID = 30
const kCGTabletProximityEventPointerType = 37
const kCGTabletProximityEventSystemTabletID = 32
const kCGTabletProximityEventTabletID = 29
const kCGTabletProximityEventVendorID = 28
const kCGTabletProximityEventVendorPointerSerialNumber = 34
const kCGTabletProximityEventVendorPointerType = 33
const kCGTabletProximityEventVendorUniqueID = 35
const kCGTailAppendEventTap = 1
const kCGTextClip = 7
const kCGTextFill = 0
const kCGTextFillClip = 4
const kCGTextFillStroke = 2
const kCGTextFillStrokeClip = 6
const kCGTextInvisible = 3
const kCGTextStroke = 1
const kCGTextStrokeClip = 5
const kCGTornOffMenuWindowLevel = 3
const kCGTornOffMenuWindowLevelKey = 6
const kCGUtilityWindowLevel = 19
const kCGUtilityWindowLevelKey = 17
const kCGWindowBackingCFNumberType = 3
const kCGWindowIDCFNumberType = 3
const kCGWindowImageBestResolution = 8
const kCGWindowImageBoundsIgnoreFraming = 1
const kCGWindowImageDefault = 0
const kCGWindowImageNominalResolution = 16
const kCGWindowImageOnlyShadows = 4
const kCGWindowImageShouldBeOpaque = 2
const kCGWindowListExcludeDesktopElements = 16
const kCGWindowListOptionAll = 0
const kCGWindowListOptionIncludingWindow = 8
const kCGWindowListOptionOnScreenAboveWindow = 2
const kCGWindowListOptionOnScreenBelowWindow = 4
const kCGWindowListOptionOnScreenOnly = 1
const kCGWindowSharingCFNumberType = 3
const kCGWindowSharingNone = 0
const kCGWindowSharingReadOnly = 1
const kCGWindowSharingReadWrite = 2


text:
#  ------------------------------------------------------------------------------------------
#  Copyright (c) Microsoft Corporation. All rights reserved.
#  Licensed under the MIT License (MIT). See LICENSE in the repo root for license information.
#  ------------------------------------------------------------------------------------------

###################################################################################################
# Exports

export @batchdef, @bcast, @batchexec
export create_bcast_batch_resource, submit_batch_job, eval_symbol_get_index!
export replace_symbol_in_expression!, create_batch_resource_from_blob_future
export create_batch_resource_for_blob_future_in_ast, assign_tasks_per_pool

###################################################################################################
# Code collection

# Collect @batchdef expressions in global vars
function append_batchdef_expression(expr)

    linefilter!(expr)
    # If expression is "include" statement, add file to batch resources

    linefilter!(expr)
    if expr.args[1] == :include || expr.args[1] == :fileinclude
        for (client, resource) in zip(__clients__, __resources__)
            create_blob_containers(client["blob_client"], [__container__]) # if not exist
            push!(resource, create_batch_resource_from_file(client["blob_client"], __container__, expr.args[2];
                verbose=__verbose__)[1])
        end
        expr.args[end] = convert(String, split(expr.args[end], "/")[end])    # remove path
    end

    # If expressing is "using" statement, collect separately in __packages__ variable, else append to __expressions__
    if expr.head == :using
        ex0 = :(@everywhere); push!(ex0.args, expr); linefilter!(ex0)
        if isnothing(__packages__)
            global __packages__ = Expr(:block, ex0)
        else
            push!(__packages__.args, ex0)
        end
    elseif expr.args[1] != :fileinclude
        ex0 = :(@everywhere); push!(ex0.args, expr); linefilter!(ex0)
        # Append expression to AST
        if isnothing(__expressions__)
            global __expressions__ = Expr(:block, ex0)
        else
            push!(__expressions__.args, ex0)
        end
    end
end

###################################################################################################
# Submit batch job

# Scan AST and create batch resource file for blob futures
function create_batch_resource_for_blob_future_in_ast(expr, pool_no, blob_futures)

    # Reached leaf
    if typeof(expr) != Expr
        return
    elseif expr.head == :call && length(expr.args) == 1
        return
    end

    for (i, arg) in enumerate(expr.args)
        if typeof(arg) == BlobFuture
            typeof(arg.blob) != BlobRef &&
                throw("BlobFuture does not contain a blob reference. Possibly, fetch! was already called on the reference.")
            create_batch_resource_from_blob_future(arg, pool_no, blob_futures)
        elseif typeof(arg) == Expr
            create_batch_resource_for_blob_future_in_ast(expr.args[i], pool_no, blob_futures)
        end
    end
end

function create_batch_resource_from_blob_future(blob_future::BlobFuture, pool_no, blob_futures)
    for blob in blob_future.blob.name
        push!(blob_futures, create_batch_resource_from_blob(__active_pools__[pool_no]["clients"]["blob_client"], blob_future.container, blob)[1])
    end
end

function create_batch_resource_for_batch_future_in_ast(expr, pool_no, batch_futures)

    # Reached leaf
    if typeof(expr) != Expr
        return
    elseif expr.head == :call && length(expr.args) == 1
        return
    end

    for (i, arg) in enumerate(expr.args)
        if typeof(arg) == BatchFuture
            typeof(arg.blob) != BlobRef &&
                throw("BatchFuture does not contain a blob reference. Possibly, fetch! was already called on the reference.")
            create_batch_resource_from_batch_future(arg, pool_no, batch_futures)
        elseif typeof(arg) == Expr
            create_batch_resource_for_batch_future_in_ast(expr.args[i], pool_no, batch_futures)
        end
    end
end

function create_batch_resource_from_batch_future(batch_future::BatchFuture, pool_no, batch_futures)
    url = create_blob_url(__active_pools__[pool_no]["clients"]["blob_client"],
        batch_future.container, [batch_future.blob.name])
    push!(batch_futures, create_batch_resource_from_blob_url(url, [batch_future.blob.name])[1])
end

function create_batch_task!(expr, pool_no, count, tasks, resources, task_ids, output, app_cmd, options)

    # Append expressions previously tagged via @batchdef
    isnothing(options) ? (task_base = "task_") : (task_base = options.task_name)
    push!(task_ids, Dict("taskname" => join([task_base, count]), "pool" => pool_no))
    if isnothing(__expressions__)
        expressions = expr
    else
        expressions = deepcopy(__expressions__)
        push!(expressions.args, expr)
    end

    # Replace return statement with serialization and create batch output resources
    filenames = []; outfiles = Array{PyObject}(undef, 0)
    funcname = expr.args[1] # name of function that is called remotely
    find_function_in_ast_and_replace_return!(expressions, funcname, filenames)
    for outfile in filenames
        push!(outfiles, create_batch_output_file(__active_pools__[pool_no]["clients"]["blob_client"],
            __active_pools__[pool_no]["credentials"]["_STORAGE_ACCOUNT_NAME"],
            __container__, outfile)...)
    end

    # Collect output blob names in Julia Futures
    future = BlobFuture(__container__, BlobRef(tuple(filenames...)), pool_no)
    push!(output, future)

    # Serialize AST
    filename = join([task_base, count, ".dat"])
    iostream = IOBuffer(); serialize(iostream, linefilter!(expressions))

    # Environment variables
    if ~isnothing(options) && options.reset_mpi == true
        env_num_nodes_per_task = "1"
        env_num_procs_per_node = "1"
    else
        env_num_nodes_per_task = __params__["_NUM_NODES_PER_TASK"]
        env_num_procs_per_node = __params__["_NUM_PROCS_PER_NODE"]
    end
    envs = create_batch_envs(
        ["FILENAME", "JULIA_DEPOT_PATH", "PYTHONPATH", "MPI_RUN", "INTER_NODE_CONNECTION", "NUM_NODES_PER_TASK",
        "NUM_PROCS_PER_NODE", "OMP_NUM_THREADS",], [filename, __params__["_JULIA_DEPOT_PATH"],
        __params__["_PYTHONPATH"], __params__["_MPI_RUN"], __params__["_INTER_NODE_CONNECTION"],
        env_num_nodes_per_task, env_num_procs_per_node, __params__["_OMP_NUM_THREADS"]])

    # Create resource file and append to resource list
    ast_resource = create_batch_resource_from_bytes(__active_pools__[pool_no]["clients"]["blob_client"], __container__, filename, iostream.data; verbose=__verbose__)
    length(__active_pools__[pool_no]["resources"]) > 0 && (ast_resource = vcat(ast_resource, __active_pools__[pool_no]["resources"]))

    # Create resource file for blob futures in AST
    blob_futures = Array{PyCall.PyObject, 1}()
    create_batch_resource_for_blob_future_in_ast(expr, pool_no, blob_futures)
    length(blob_futures) > 0 && (ast_resource = vcat(ast_resource, blob_futures))

    # Create resource file for batch futures in AST
    batch_futures = Array{PyCall.PyObject, 1}()
    create_batch_resource_for_batch_future_in_ast(expr, pool_no, batch_futures)
    length(batch_futures) > 0 && (ast_resource = vcat(ast_resource, batch_futures))

    # Create tasks and submit job
    task_constraints = create_task_constraint(max_task_retry_count=parse(Int, __params__["_NUM_RETRYS"]))
    num_nodes_per_task=parse(Int,__params__["_NUM_NODES_PER_TASK"])
    __params__["_CONTAINER"] == "None" ? (docker_container = nothing) : (docker_container = __params__["_CONTAINER"])
    push!(tasks, create_batch_task(resource_files=vcat(resources, ast_resource), application_cmd=app_cmd, task_constraints=task_constraints,
        environment_variables=envs, output_files=outfiles, taskname=task_ids[end]["taskname"], num_nodes_per_task=num_nodes_per_task,
        docker_container=docker_container))
end

# Split list of task expressions according to schedule
function split_expressions(expression_list, task_list_per_pool)

    # Find pools with tasks in it
    idx = findall(i -> i != 0:0, task_list_per_pool)
    num_pools_active = length(idx)
    expressions_per_pool = Array{Any}(undef, num_pools_active)
    pool_no = Array{Any}(undef, num_pools_active)

    for i=1:num_pools_active
        expressions_per_pool[i] = expression_list[task_list_per_pool[idx[i]]]
        pool_no[i] = idx[i]
    end
    return expressions_per_pool, idx
end

# Assign range of tasks to each batch pool
function assign_tasks_per_pool(num_tasks; strategy="chunk", options=nothing)

    # Distribute tasks among available pools
    ~isnothing(options) && ~isnothing(options.pool) ? pool_no = options.pool : pool_no = nothing
    num_pools = length(__active_pools__)
    if num_pools == 0
        @warn "No active pools found. Tasks will be submitted but cannot be executed."
        num_pools = parse(Int, __params__["_POOL_COUNT"])
    end
    if strategy == "chunk"

        # Split into even chunks
        min_tasks = Int(floor(num_tasks/num_pools))
        leftover_tasks = mod(num_tasks, num_pools)
        tasks_per_worker = Array{Any}(undef, num_pools)
        for i=1:num_pools
            tasks_per_worker[i] = min_tasks
            if i <= leftover_tasks
                tasks_per_worker[i] += 1
            end
        end

        task_list_per_pool = Array{Any}(undef, num_pools)
        if ~isnothing(pool_no)
            # Assign all tasks to specified pool
            for j=1:num_pools
                task_list_per_pool[j] = 0:0
            end
            task_list_per_pool[pool_no] = 1:num_tasks
        else
            count = 0
            for i=1:num_pools
                if count < num_tasks
                    task_list_per_pool[i] = count+1:count+tasks_per_worker[i]
                else
                    task_list_per_pool[i] = 0:0
                end
                count += tasks_per_worker[i]
            end
        end
    else
        throw("Specified stragety not supported.")
    end
    return task_list_per_pool
end

# Submit multi-task batch job
function submit_batch_job(expression_list; options=nothing)

    num_tasks = length(expression_list)

    # CMD, job id and priority
    app_cmd = "/bin/bash -c \'set -e; set -o pipefail; \$AZ_BATCH_TASK_WORKING_DIR/application-cmd; wait\'"
    ~isnothing(options) ? (base_name = options.job_name) : (base_name = __params__["_JOB_ID"])
    job_base = join([base_name, "_", randstring(8)])
    ~isnothing(options) ? (priority = options.priority) : (priority = 0)

    # Split expressions among available batch pools
    task_list_per_pool =  assign_tasks_per_pool(num_tasks; strategy="chunk", options=options)
    expressions_per_pool, pool_numbers = split_expressions(expression_list, task_list_per_pool)

    output = []; task_ids = []; count = 1; job_ids = []
    for (i, expressions) in enumerate(expressions_per_pool)

        pool_no = pool_numbers[i]
        push!(job_ids, join([job_base, "_", i]))
        create_batch_job(__active_pools__[pool_no]["clients"]["batch_client"], job_ids[end],
            __active_pools__[pool_no]["pool_id"]; uses_task_dependencies=false, priority=priority, verbose=__verbose__)
        create_blob_containers(__active_pools__[pool_no]["clients"]["blob_client"], [__container__])

        # Add Julia runtime and cmd to batch resource list
        resources = Array{PyObject}(undef, 0)
        push!(resources, create_batch_resource_from_file(__active_pools__[pool_no]["clients"]["blob_client"],
            __container__, joinpath(dirname(pathof(AzureClusterlessHPC)), "runtime/application-cmd");
            verbose=__verbose__)[1])
        push!(resources, create_batch_resource_from_file(__active_pools__[pool_no]["clients"]["blob_client"],
            __container__, joinpath(dirname(pathof(AzureClusterlessHPC)), "runtime/batch_runtime.jl");
            verbose=__verbose__)[1])

        # Serialize expressions with "using ..." and create batch resource
        if ~isnothing(__packages__)
            iostream = IOBuffer(); serialize(iostream, __packages__)
            push!(resources, create_batch_resource_from_bytes(__active_pools__[pool_no]["clients"]["blob_client"],
                __container__, "packages.dat", iostream.data; verbose=__verbose__)[1])
        end

        # Create tasks for each batch pool
        tasks = []
        @sync begin
            for (j, expr) in enumerate(expressions)
                create_batch_task!(expr, pool_no, count, tasks, resources, task_ids, output, app_cmd, options)
                count += 1
            end
        end
        if ~isnothing(__active_pools__[pool_no]["clients"]["batch_client"])
            __active_pools__[pool_no]["clients"]["batch_client"].task.add_collection(job_ids[end], tasks)
        end
    end
    return BatchController(job_ids, task_ids, length(expression_list), output)
end

#######################################################################################################################
# Capture runtime args and create job submission

# Set index methods for expressions
function setindex!(expr::Expr, val, i)
    expr.args[i] = val
end

function setindex!(expr::Expr, val, i::Array)
    if length(i) == 1
        expr.args[i[1]] = val
    elseif length(i) == 2
        expr.args[i[1]].args[i[2]] = val
    elseif length(i) == 3
        expr.args[i[1]].args[i[2]].args[i[3]] = val
    elseif length(i) == 4
        expr.args[i[1]].args[i[2]].args[i[2]].args[i[4]] = val
    elseif length(i) == 5
        expr.args[i[1]].args[i[2]].args[i[2]].args[i[4]].args[i[5]] = val
    end
end

# Find symbols
function locate_symbols!(expr, level, locations, symbols; symbol_exception=nothing)

    # Reached leaf
    if typeof(expr) != Expr
        return
    elseif expr.head == :call && length(expr.args) == 1
        return
    end

    # Recursively step through AST and find symbols that are neither function nor kwarg names
    (expr.head == :call || expr.head == :kw) ? (start=2) : (start=1)
    for j=start:length(expr.args)
        if typeof(expr.args[j]) == Symbol
            # Dont add iterate from pmap to collection
            if isnothing(symbol_exception) || expr.args[j] != symbol_exception
                push!(symbols, expr.args[j])
                push!(locations, vcat(level, j))
            end
            # To do: don't replace symbols of functions passed as args
        elseif typeof(expr.args[j]) == Expr
            new_level = deepcopy(level)
            push!(new_level, j)
            locate_symbols!(expr.args[j], new_level, locations, symbols; symbol_exception=symbol_exception)
        end
    end
end

# Replace given symbol with value in expression
function replace_symbol_in_expression!(expr, symbol, val)
    if typeof(expr) == Expr
        for (i, arg) in enumerate(expr.args)
            if typeof(arg) == Expr
                replace_symbol_in_expression!(expr.args[i], symbol, val)
            elseif typeof(arg) == Symbol
                if arg == symbol
                    expr.args[i] = val
                end
            end
        end
    end
end

# If symbol is indexed, extract indexed value
function eval_symbol_get_index!(expr)
    if typeof(expr) == Expr
        for (i, arg) in enumerate(expr.args)
            if typeof(arg) == Expr && arg.head == :ref
                expr.args[i] = arg.args[1][arg.args[2:end]...]
            elseif typeof(arg) == Expr
                eval_symbol_get_index!(expr.args[i])
            end
        end
    end
end

# Create list of expressions for multi-task batch job
function create_multi_task_expression_list_pmap(expr; options=nothing)

        # Index symbol and pmap argument collection
        index_symbol = expr.args[2].args[1]
        arg_collection = expr.args[3]

        # Find all symbols and their locations
        level = []; indices = []; symbols = []
        expr_i = deepcopy(expr.args[2].args[2].args[2])
        locate_symbols!(expr_i, level, indices, symbols; symbol_exception=index_symbol)

        # Base quote: deserialize expression + loop over tasks
        iostream = IOBuffer(); AzureClusterlessHPC.serialize(iostream, expr)
        expr_base = quote
            _expr = AzureClusterlessHPC.deserialize(IOBuffer($iostream.data))
            _expr_collection = []
            _index_symbol = _expr.args[2].args[1]
            _arg_collection = $arg_collection
            for (_i, _arg) in enumerate(_arg_collection)
                _expr_i = deepcopy(_expr.args[2].args[2].args[2])
                replace_symbol_in_expression!(_expr_i, _index_symbol, _arg)
            end
        end
        linefilter!(expr_base)

        # Loop over symbols
        for j=1:length(symbols)

            # Replace symbol
            symbol = symbols[j]
            idx = indices[j]
            expr_new = quote
                _expr_i[$idx] = $symbol
            end

            # Insert quote into task loop
            linefilter!(expr_new)
            push!(expr_base.args[5].args[2].args, expr_new.args...)
        end

        # Add expr to collection (inside task loop)
        expr_new = quote
            eval_symbol_get_index!(_expr_i)
            push!(_expr_collection, _expr_i)
        end
        linefilter!(expr_new)
        push!(expr_base.args[5].args[2].args, expr_new.args...)

        # Submit batch job for expression
        expr_new = quote
            @time submit_batch_job(_expr_collection; options=$options)
        end
        linefilter!(expr_new)
        push!(expr_base.args, expr_new.args...)
        return expr_base
end

# Create expression for single-task batch job
function create_single_task_expression_list(expr; options=nothing)

    # Find all symbols and their locations
    level = []; indices = []; symbols = []
    locate_symbols!(expr, level, indices, symbols)

    # Base quote: deserialize expression
    iostream = IOBuffer(); AzureClusterlessHPC.serialize(iostream, expr)
    expr_base = quote
        _expr = AzureClusterlessHPC.deserialize(IOBuffer($iostream.data))
    end
    linefilter!(expr_base)

    # Loop over symbols
    for j=1:length(symbols)

        # Replace symbol
        symbol = symbols[j]
        idx = indices[j]
        expr_new = quote
            _expr[$idx] = $symbol
        end

        # Append to quote
        linefilter!(expr_new)
        push!(expr_base.args, expr_new.args...)
    end

    # Submit batch job for expression
    expr_new = quote
        eval_symbol_get_index!(_expr)   # eval get index functions
        submit_batch_job([_expr]; options=$options)
    end
    linefilter!(expr_new)
    push!(expr_base.args, expr_new.args...)
    return expr_base
end

# Create list of batch tasks from expression
function create_expression_to_submit_batch_job(expr; options=nothing)
    if expr.args[1] == :pmap
        expr_base = create_multi_task_expression_list_pmap(expr; options=options)
    else
        expr_base = create_single_task_expression_list(expr; options=options)
    end
    return expr_base
end

#######################################################################################################################
# Output

# Find "write" statements and collect them in "filenames"
function find_output_files!(expr, filenames)

    # Reached leaf
    if typeof(expr) != Expr
        return
    end

    # Recursively step through AST
    for i=1:length(expr.args)
        if typeof(expr.args[i]) == Expr
            if expr.args[i].args[1] == :write
                filename = randstring(12)   # create new random output filename as original filename may contain symbols
                expr.args[i].args[2] = filename
                push!(filenames, filename)
            else
                find_output_files!(expr.args[i], filenames)
            end
        end
    end
end

# Find return statement in expression and replace w/ serialization
function replace_return_with_serialization!(expr, filelist)
    typeof(expr) != Expr && return

    # Check if return is in current symbol block
    found_return_in_block = false
    idx = nothing
    for (i, arg) in enumerate(expr.args)
        if typeof(arg) == Expr && arg.head == :return
            found_return_in_block = true
            idx = i
        end
    end

    if found_return_in_block
        # Multiple return arguments or expression
        if typeof(expr.args[idx].args[1]) == Expr
            # Loop over return arguments
            for (i, argout) in enumerate(expr.args[idx].args[1].args)

                # Create random filename and add to collection
                filename = randstring(12)
                push!(filelist, filename)

                # Insert serialization at location of return statement
                insert!(expr.args, idx + (3*i-2), :(iostream = open($filename, "w")))
                insert!(expr.args, idx + (3*i-1), :(serialize(iostream, $argout)))
                insert!(expr.args, idx + 3*i, :(close(iostream)))
            end
        else
            argout = expr.args[idx].args[1]

            # Create random filename and add to collection
            filename = randstring(12)
            push!(filelist, filename)

            # Insert serialization at location of return statement
            insert!(expr.args, idx+1, :(iostream = open($filename, "w")))
            insert!(expr.args, idx+2, :(serialize(iostream, $argout)))
            insert!(expr.args, idx+3, :(close(iostream)))
        end
        # Remove return statement from block
        popat!(expr.args, idx)
    else
        # Step through AST
        for i=1:length(expr.args)
            if typeof(expr.args[i]) == Expr
                replace_return_with_serialization!(expr.args[i], filelist)
            end
        end
    end
end

# Replace return statement for given function with serialization
function find_function_in_ast_and_replace_return!(expr, fname, filelist)

    # Reached AST leaf
    typeof(expr) != Expr && return

    if expr.head == :function && expr.args[1].args[1] == fname
        # Found function with given name. Now replace return w/ serialization
        replace_return_with_serialization!(expr, filelist)
    else
        # Step through AST
        for i=1:length(expr.args)
            if typeof(expr.args[i]) == Expr
                find_function_in_ast_and_replace_return!(expr.args[i], fname, filelist)
            end
        end
    end
end

###################################################################################################
# Bcast

# Create batch resource for broadcasted variable
function create_bcast_batch_resource(blob_name, binary; container=nothing)

    # If bcast is overwrite of exisiting resources -> remove original one
    isnothing(container) && (container = __container__)
    for (client, resource) in zip(__clients__, __resources__)
        resource_idx = findall(x -> x.file_path == blob_name, resource)
        if ~isempty(resource_idx)
            popat!(resource, resource_idx[1]) # remove original resource
        end
        upload_bytes_to_container(client["blob_client"], container, blob_name, binary; verbose=__verbose__)
    end

    # Return future w/ blob name
    future = BatchFuture(container, BlobRef(blob_name))

    return future
end

# Broadcasting
function bcast_expression(expr::Symbol; container=nothing)

    filename = string(expr)
    isnothing(container) && (container = __container__)
    for client in __clients__
        create_blob_containers(client["blob_client"], [container]) # if not exist
    end

    # Return expression
    expr_out = quote
        _iobuff = IOBuffer()
        AzureClusterlessHPC.serialize(_iobuff, $expr)   # eval in local scope
        _binary = _iobuff.data

        # Create batch resource
        _blob_name = join([$filename, ".dat"])
        create_bcast_batch_resource(_blob_name, _binary; container=$container)
    end
    return esc(expr_out)
end

###################################################################################################
# Macros

# Macro to upload stuff to batch
"""
    @batchdef expr

 Execute an expression under Main and on the batch workers of (future) batch jobs.

 *Usage*:

 - `@batchdef bar = 1` will define Main.bar on the current process and on all batch workers of
    batch jobs that are executed from the same session.

 See also: [`@batchexec`](@ref), [`batch_clear`](@ref)
"""
macro batchdef(expr)

    # If expression is begin ... end block, split into single expressions
    if expr.head == :block
        for i=1:length(expr.args)
            if typeof(expr.args[i]) == Expr
                append_batchdef_expression(expr.args[i])
            end
        end
    else
        append_batchdef_expression(expr)
    end

    return esc(expr)    # esc makes expression available in caller scope
end

# Broadcasting: upload expr to blob and add to batch resource file
"""
    batch_future = @bcast expr

 Broadcast an expression to the workers of batch jobs that are executed from the current session.
 `@bcast` uploades the tagged expression, creates a resource file for (future) batch jobs and returns
 a `BatchFuture` of the expression.

 *Usage*:

 - `A = randn(2, 2); bfuture = @bcast A` will broadcast `A` to all batch workers/tasks of batch jobs that
    are executed from within the same session.

 - The return argument is a `BatchFuture` that contains the name of the blob to which the expression was uploaded.
    To fetch the data associated with `batch_future`, use the `fetch` or `fetch!` function on it: `A = fetch(batch_future)`.

 See also: [`fetch`](@ref), [`fetch!`](@ref), [`@batchdef`](@ref), [`@batchexec`](@ref), [`batch_clear`](@ref)
"""
macro bcast(expr)
    return bcast_expression(expr; container=nothing)
end

macro bcast(expr, container)
    return bcast_expression(expr; container=container)
end

# Run given expression as batch job
"""
    `batch_controller = @batchexec expr`

 Create a closure around an expression and run it automatically as a batch job, returning a future to the result.
 If the expression is the `pmap` function, it is executed as a multi-task batch job.

 *Usage*:

 - `@batchexec print("Hello world")` executes the `print` function as a (single-task) batch job.

 - `@batchexec pmap(i -> print("Hello from \$i"), 1:10)` executes the `print` function as a multi-task batch job.

 - The return argument is a batch_controller that provides access to return arguments of the batch job,
    as well as basic job and pool functionalities. To fetch the output, use the `fetch` or `fetch!` function on
    it. These operations are blocking and wait for the completion of the batch job.

 See also:  [`BatchController`](@ref), [`fetch`](@ref), [`fetch!`](@ref), [`@batchdef`](@ref), [`batch_clear`](@ref)
 """
macro batchexec(expr)
    expr_out = create_expression_to_submit_batch_job(expr)
    return esc(expr_out)
end

macro batchexec(expr, options)
    expr_out = create_expression_to_submit_batch_job(expr; options=options)
    return esc(expr_out)
end


text:
module Filters
using ..Unwrap
using ..Util
using Polynomials: LaurentPolynomial, Polynomial, coeffs, derivative, fromroots, roots

import Base: *
using LinearAlgebra: I, mul!, rmul!
using Statistics: middle
import ..DSP: filt, filt!, optimalfftfiltlength, os_fft_complexity, SMALL_FILT_CUTOFF
import Compat
using FFTW

include("coefficients.jl")
export FilterCoefficients,
        ZeroPoleGain,
        PolynomialRatio,
        Biquad,
        SecondOrderSections,
        coefa,
        coefb

include("filt.jl")
export  DF2TFilter,
        filtfilt,
        tdfilt,
        tdfilt!,
        fftfilt,
        fftfilt!

include("design.jl")
export  FilterType,
        Butterworth,
        Chebyshev1,
        Chebyshev2,
        Elliptic,
        Lowpass,
        Highpass,
        Bandpass,
        Bandstop,
        analogfilter,
        digitalfilter,
        iirnotch,
        kaiserord,
        FIRWindow,
        resample_filter

include("response.jl")
export  freqresp,
        phaseresp,
        grpdelay,
        impresp,
        stepresp

include("stream_filt.jl")
export  FIRFilter,
        outputlength,
        inputlength,
        reset!,
        resample,
        setphase!,
        timedelay

include("remez_fir.jl")
export  remez,
        RemezFilterType,
        filter_type_bandpass,
        filter_type_differentiator,
        filter_type_hilbert

end


text:
include("../src/geom_module.jl")

module linetest
using ..geom
using PyPlot
PyPlot.pygui(true)

function test_line()

end

end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
module TestPoints

using Test
using TreeParzen

objective(params) = params[:x]^2 + params[:y]^2
space = Dict(:x => HP.Uniform(:x, -10.0, 10.0), :y => HP.Uniform(:y, -10.0, 10.0))
points = [Dict(:x => 0.0, :y => 0.0), Dict(:x => 1.0, :y => 1.0)]

best = fmin(objective, space, 10, points)
@test iszero(best[:x])
@test iszero(best[:y])

# test to catch ArgumentError when number of iterations is lower than the number of points
@test_throws ArgumentError fmin(objective, space, 1, points)

end
true


text:
#=
given a nonnegative matrix with missing entries, find nonnegative integer values for the
missing entries to minimize the spectral or nuclear norm, subject to linear constraints
=#

struct MatrixCompletion <: ExampleInstance
    nuclear_obj::Bool # use nuclear norm objective, else spectral norm
    symmetric::Bool # use symmetric matrix, else rectangular matrix
    nrow::Int # number of rows in matrix
    ncol::Int # not used in symmetric case
    use_nat::Bool # use spectral/nuclear norm cone, else PSD cone formulation
end

function build(inst::MatrixCompletion)
    nrow = inst.nrow
    @assert nrow >= 2
    ncol = (inst.symmetric ? nrow : inst.ncol)
    @assert ncol >= nrow
    # set some options
    sparsity = min(0.7, 1 - inv(sqrt(nrow * ncol))) # fraction of values that are known
    max_rank = div(nrow, 2) # maximum rank of random matrix

    # generate data
    pattern = sprand(Bool, nrow, ncol, sparsity)
    if inst.symmetric
        len = Cones.svec_length(nrow)
        triu!(pattern)
        A0 = rand(0:3, nrow, max_rank)
        A0 = A0 * A0'
    else
        len = nrow * ncol
        A0 = rand(0:3, nrow, max_rank) * rand(0:3, max_rank, ncol)
    end
    (rows, cols, _) = findnz(pattern)
    Avals = [A0[r, c] for (r, c) in zip(rows, cols)]
    num_unknown = len - length(Avals)
    @assert num_unknown >= 1

    # build model
    model = JuMP.Model()
    JuMP.@variable(model, t)
    JuMP.@objective(model, Min, t / nrow)

    JuMP.@variable(model, x[1:num_unknown], Int)
    JuMP.@constraint(model, x .>= 0)
    JuMP.@constraint(model, x .<= 9)
    JuMP.@constraint(model, sum(x) >= length(x))

    X = Matrix{JuMP.AffExpr}(undef, nrow, ncol)
    for (r, c, a) in zip(rows, cols, Avals)
        X[r, c] = a
    end
    k = 1
    for c in 1:ncol, r in 1:(inst.symmetric ? c : nrow)
        if !isassigned(X, r, c)
            X[r, c] = x[k]
            k += 1
        end
    end
    if inst.symmetric
        LinearAlgebra.copytri!(X, 'U')
        add_specnuc = add_symm_specnuc
    else
        add_specnuc = add_nonsymm_specnuc
    end
    @assert k == length(x) + 1
    X ./= nrow # rescale for numerics

    add_specnuc(inst.use_nat, inst.nuclear_obj, t, X, model)

    # save for use in tests
    model.ext[:x] = x
    model.ext[:X] = X

    return model
end

function test_extra(inst::MatrixCompletion, model::JuMP.Model)
    stat = JuMP.termination_status(model)
    @test stat in OPT_OR_LIMIT
    (stat == MOI.OPTIMAL) || return

    tol = eps()^0.2
    x = JuMP.value.(model.ext[:x])
    @test x ≈ round.(Int, x) atol = tol rtol = tol
    # check objective
    X = JuMP.value.(model.ext[:X])
    s = (inst.symmetric ? abs.(eigvals(Symmetric(X, :U))) : svdvals(X))
    snorm = (inst.nuclear_obj ? sum(s) : maximum(s))
    @test JuMP.objective_value(model) * inst.nrow ≈ snorm atol = tol rtol = tol
    return
end


text:
"""
`save_data(df::DataFrame, args...)`

Save a DataFrame as a csv file into the `data` folder. `args...` specify the
subdirecotires in `data`.
"""
function save_data(df::DataFrame, args...)

    # Build a path starting with "data" and continuig based on the passed args
    filename = joinpath("data", args...)

    # If a file with the same name already exists, throw an error
    if isfile(filename)
        error("File already exists!")
    end

    # Create all subdirecotires that do not yet exist
    dir, _ = splitdir(filename)
    if !isdir(dir)
        mkpath(dir)
    end

    # Save the dataframe
    return CSV.write(filename, df)

end

"""
`load_data(args...)`

Load a DataFrame from a csv file located in the `data` folder. `args...` specify
the subdirecotires in `data`.
"""
function load_data(args...)

    # Build a path starting with "data" and continuig based on the passed args
    filename = joinpath(@__DIR__, "..", "data", args...)

    # Load a csv file
    file = CSV.File(filename)

    # Convert the file into a dataframe
    df = DataFrame(file)

    return df

end

"""
`save_figure(fig, args...; kwargs...)`

Save a figure into the `figures` folder. `args...` specify the subdirecotires in
`figures` and `kwargs...` are passed to `fig.savefig`.
"""
function save_figure(fig, args...; kwargs...)

    # Build a path starting with "figures" and continuig based on the passed args
    filename = joinpath("figures", args...)

    # If a file with the same name already exists, throw an error
    if isfile(filename)
        error("File already exists!")
    end

    # Create all subdirecotires that do not yet exist
    dir, _ = splitdir(filename)
    if !isdir(dir)
        mkpath(dir)
    end

    # Save the figure
    return fig.savefig(filename; kwargs...)

end

max_stars_repo_name: vkumpost/stoosc

max_issues_repo_name: vkumpost/stoosc

max_forks_repo_name: vkumpost/stoosc


text:
using Gurobi, MathProgBase

# This problem should error out on Gurobi Versions 6.5.2 and earlier,
# but pass on later versions

m = MathProgBase.LinearQuadraticModel(GurobiSolver(OutputFlag=0))

# Max           z
# s.t.  x + y     >= 0
#       x + y     >= 0
#      2x - y     >= 0
#       x     - z >= 0
#               z >= -10
# x ∈ [0, 1], y ∈ [0, ∞), z \in (-∞, 1]
MathProgBase.loadproblem!(m, [1 1 0; 1 1 0; 2 -1 0; 1 0 -1; 0 0 1], [0, 0, -Inf],[1,Inf,1], [0,0,1], [0,0,0,0,-10], [Inf,Inf,Inf,Inf,Inf], :Max)

MathProgBase.optimize!(m)

MathProgBase.addconstr!(m, [1, 2, 3], [-1.0000019, 1., -1.], 0., Inf)

MathProgBase.optimize!(m)

obj = MathProgBase.getobj(m)
MathProgBase.setobj!(m, obj)

MathProgBase.optimize!(m)