text:
#=
    NEMO: Next Energy Modeling system for Optimization.
    https://github.com/sei-international/NemoMod.jl

    Copyright © 2018: Stockholm Environment Institute U.S.

    File description: Other NEMO functions.
=#

"""
    logmsg(msg::String, suppress=false, dtm=now()::DateTime)

Prints a log message (`msg`) to `STDOUT`. The message is suppressed if `suppress == true`. `dtm`
determines the date and time included in the printed message.

# Examples
```jldoctest
julia> using Dates

julia> logmsg("Test message", false, DateTime(2020))
2020-01-Jan 00:00:00.000 Test message
```
"""
function logmsg(msg::String, suppress=false, dtm=now()::DateTime)
    suppress && return
    println(stdout, Dates.format(dtm, @dateformat_str "YYYY-dd-u HH:MM:SS.sss ") * msg)
end  # logmsg(msg::String)

"""
    getconfig(quiet::Bool = false)

Reads in NEMO's configuration file, which should be in the Julia working directory,
in `ini` format, and named `nemo.ini` or `nemo.cfg`. See the sample file in the NEMO
package's `utils` directory for more information.

# Arguments
- `quiet::Bool = false`: Suppresses low-priority status messages (which are otherwise printed to STDOUT).
"""
function getconfig(quiet::Bool = false)
    configpathini::String = joinpath(pwd(), "nemo.ini")
    configpathcfg::String = joinpath(pwd(), "nemo.cfg")
    configpath::String = ""

    if isfile(configpathini)
        configpath = configpathini
    elseif isfile(configpathcfg)
        configpath = configpathcfg
    else
        return nothing
    end

    try
        local conffile::ConfParse = ConfParse(configpath)
        parse_conf!(conffile)
        logmsg("Read NEMO configuration file at " * configpath * ".", quiet)
        return conffile
    catch
        logmsg("Could not parse NEMO configuration file at " * configpath * ". Please verify format.", quiet)
    end
end  # getconfig(quiet::Bool = false)

"""
    getconfigargs!(configfile::ConfParse,
        calcyears::Array{Int,1},
        varstosavearr::Array{String,1},
        targetprocs::Array{Int,1},
        bools::Array{Bool,1},
        ints::Array{Int,1},
        quiet::Bool)

Loads run-time arguments for `calculatescenario()` from a configuration file.

# Arguments
- `configfile::ConfParse`: Configuration file. This argument is not changed by the function.
- `calcyears::Array{Int,1}`: `calculatescenario()` `calcyears` argument. Values specified
    in the configuration file for `[calculatescenarioargs]/calcyears` (including an empty or null value)
    replace what is in this array.
- `varstosavearr::Array{String,1}`: Array representation of `calculatescenario()` `varstosave` argument.
    Values specified in the configuration file for `[calculatescenarioargs]/varstosave` are added to this array.
- `targetprocs::Array{Int,1}`: `calculatescenario()` `targetprocs` argument. Values specified in the
    configuration file for `[calculatescenarioargs]/targetprocs` are added to this array.
- `bools::Array{Bool,1}`: Array of Boolean arguments for `calculatescenario()`: `restrictvars`, `reportzeros`,
    `continuoustransmission`, `forcemip`, and `quiet`, in that order. Values specified in the configuration file
    for `[calculatescenarioargs]/restrictvars`, `[calculatescenarioargs]/reportzeros`,
    `[calculatescenarioargs]/continuoustransmission`, `[calculatescenarioargs]/forcemip`, and `[calculatescenarioargs]/quiet`
    overwrite what is in this array (provided the values can be parsed as `Bool`).
- `ints::Array{Int,1}`: Array of `Int` arguments for `calculatescenario()`: `numprocs`. Values specified in the
    configuration file for `[calculatescenarioargs]/numprocs` overwrite what is in this array (provided
    the values can be parsed as `Int`).
- `quiet::Bool`: Suppresses low-priority status messages (which are otherwise printed to `STDOUT`).
    This argument is not changed by the function.
"""
function getconfigargs!(configfile::ConfParse, calcyears::Array{Int,1}, varstosavearr::Array{String,1},
    targetprocs::Array{Int,1}, bools::Array{Bool,1}, ints::Array{Int,1}, quiet::Bool)

    if haskey(configfile, "calculatescenarioargs", "calcyears")
        try
            calcyearsconfig = retrieve(configfile, "calculatescenarioargs", "calcyears")
            calcyearsconfigarr::Array{Int,1} = Array{Int,1}()  # calcyearsconfig converted to an Int array

            if typeof(calcyearsconfig) == String
                calcyearsconfigarr = [Meta.parse(calcyearsconfig)]
            else
                # calcyearsconfig should be an array of strings
                calcyearsconfigarr = [Meta.parse(v) for v in calcyearsconfig]
            end

            # This separate operation is necessary in order to modify calcyears
            empty!(calcyears)
            append!(calcyears, calcyearsconfigarr)
            logmsg("Read calcyears argument from configuration file.", quiet)
        catch e
            logmsg("Could not read calcyears argument from configuration file. Error message: " * sprint(showerror, e) * ". Continuing with NEMO.", quiet)
        end
    end

    if haskey(configfile, "calculatescenarioargs", "varstosave")
        try
            varstosaveconfig = retrieve(configfile, "calculatescenarioargs", "varstosave")

            if typeof(varstosaveconfig) == String
                union!(varstosavearr, [lowercase(varstosaveconfig)])
            else
                # varstosaveconfig should be an array of strings
                union!(varstosavearr, [lowercase(v) for v in varstosaveconfig])
            end

            logmsg("Read varstosave argument from configuration file.", quiet)
        catch e
            logmsg("Could not read varstosave argument from configuration file. Error message: " * sprint(showerror, e) * ". Continuing with NEMO.", quiet)
        end
    end

    if haskey(configfile, "calculatescenarioargs", "numprocs")
        try
            ints[1] = Meta.parse(lowercase(retrieve(configfile, "calculatescenarioargs", "numprocs")))
            logmsg("Read numprocs argument from configuration file.", quiet)
        catch e
            logmsg("Could not read numprocs argument from configuration file. Error message: " * sprint(showerror, e) * ". Continuing with NEMO.", quiet)
        end
    end

    if haskey(configfile, "calculatescenarioargs", "targetprocs")
        try
            targetprocsconfig = retrieve(configfile, "calculatescenarioargs", "targetprocs")

            if typeof(targetprocsconfig) == String
                union!(targetprocs, [Meta.parse(targetprocsconfig)])
            else
                # targetprocsconfig should be an array of strings
                union!(targetprocs, [Meta.parse(v) for v in targetprocsconfig])
            end

            logmsg("Read targetprocs argument from configuration file.", quiet)
        catch e
            logmsg("Could not read targetprocs argument from configuration file. Error message: " * sprint(showerror, e) * ". Continuing with NEMO.", quiet)
        end
    end

    if haskey(configfile, "calculatescenarioargs", "restrictvars")
        try
            bools[1] = Meta.parse(lowercase(retrieve(configfile, "calculatescenarioargs", "restrictvars")))
            logmsg("Read restrictvars argument from configuration file.", quiet)
        catch e
            logmsg("Could not read restrictvars argument from configuration file. Error message: " * sprint(showerror, e) * ". Continuing with NEMO.", quiet)
        end
    end

    if haskey(configfile, "calculatescenarioargs", "reportzeros")
        try
            bools[2] = Meta.parse(lowercase(retrieve(configfile, "calculatescenarioargs", "reportzeros")))
            logmsg("Read reportzeros argument from configuration file.", quiet)
        catch e
            logmsg("Could not read reportzeros argument from configuration file. Error message: " * sprint(showerror, e) * ". Continuing with NEMO.", quiet)
        end
    end

    if haskey(configfile, "calculatescenarioargs", "continuoustransmission")
        try
            bools[3] = Meta.parse(lowercase(retrieve(configfile, "calculatescenarioargs", "continuoustransmission")))
            logmsg("Read continuoustransmission argument from configuration file.", quiet)
        catch e
            logmsg("Could not read continuoustransmission argument from configuration file. Error message: " * sprint(showerror, e) * ". Continuing with NEMO.", quiet)
        end
    end

    if haskey(configfile, "calculatescenarioargs", "forcemip")
        try
            bools[4] = Meta.parse(lowercase(retrieve(configfile, "calculatescenarioargs", "forcemip")))
            logmsg("Read forcemip argument from configuration file.", quiet)
        catch e
            logmsg("Could not read forcemip argument from configuration file. Error message: " * sprint(showerror, e) * ". Continuing with NEMO.", quiet)
        end
    end

    if haskey(configfile, "calculatescenarioargs", "quiet")
        try
            bools[5] = Meta.parse(lowercase(retrieve(configfile, "calculatescenarioargs", "quiet")))
            logmsg("Read quiet argument from configuration file. Value in configuration file will be used from this point forward.", quiet)
        catch e
            logmsg("Could not read quiet argument from configuration file. Error message: " * sprint(showerror, e) * ". Continuing with NEMO.", quiet)
        end
    end
end  # getconfigargs!(configfile::ConfParse, varstosavearr::Array{String,1}, targetprocs::Array{Int,1}, bools::Array{Bool,1}, quiet::Bool)

"""
    setsolverparamsfromcfg(configfile::ConfParse, jumpmodel::JuMP.Model, quiet::Bool)

Sets solver parameters specified in a configuration file. Reports which parameters were and were not successfully set.

# Arguments
- `configfile::ConfParse`: Configuration file. Solver parameters should be in `parameters` key within `solver` block,
    specified as follows: `parameter1=value1`, `parameter2=value2`, ...
- `jumpmodel::JuMP.Model`: JuMP model in which solver parameters will be set.
- `quiet::Bool`: Suppresses low-priority status messages (which are otherwise printed to `STDOUT`).
"""
function setsolverparamsfromcfg(configfile::ConfParse, jumpmodel::JuMP.Model, quiet::Bool)
    if haskey(configfile, "solver", "parameters")
        local retrieved_params = retrieve(configfile, "solver", "parameters")  # May be a string (if there's one parameter-value pair) or an array of strings (if there are more than one, comma-delimited parameter-value pairs)
        local params_set::String = ""  # List of parameters that were successfully set; used in status message

        # Add singleton parameters to an array to facilitate processing
        if typeof(retrieved_params) == String
            retrieved_params = [retrieved_params]
        end

        # BEGIN: Process each parameter-value pair in retrieved_params.
        for e in retrieved_params
            local split_p::Array{String,1} = [strip(pp) for pp = split(e, "="; keepempty=false)]  # Current parameter specification split on =; split_p[1] = parameter name, split_p[2] = parameter value

            if length(split_p) != 2
                logmsg("Could not set solver parameter specified in configuration file as " * e * ". Continuing with NEMO.", quiet)
                continue
            end

            local val  # Parsed value for current parameter

            # Int, float, and Boolean values - order is important since an integer value will parse to a float
            for t in [Int, Float64, Bool]
                val = tryparse(t, split_p[2])

                if val != nothing
                    break
                end
            end

            if val == nothing
                # String value
                val = split_p[2]
            end

            # Call set_optimizer_attribute for parameter and value
            try
                set_optimizer_attribute(jumpmodel, split_p[1], val)
                params_set *= split_p[1] * ", "
            catch
                logmsg("Could not set solver parameter specified in configuration file as " * e * ". Continuing with NEMO.", quiet)
            end
        end
        # END: Process each parameter-value pair in retrieved_params.

        # Report parameters that were set
        if length(params_set) != 0
            logmsg("Set following solver parameters using values in configuration file: " * params_set[1:prevind(params_set, lastindex(params_set), 2)] * ".", quiet)
        end
    end  # haskey(configfile, "solver", "parameters")
end  # setsolverparamsfromcfg(configfile::ConfParse, jumpmodel::JuMP.Model, quiet::Bool)

"""
    translatesetabb(a::String)

Translates a set abbreviation into the set's name.

# Examples
```jldoctest
julia> translatesetabb(y)
""YEAR""
```
"""
function translatesetabb(a::String)
    if a == "y"
        return "YEAR"
    elseif a == "t"
        return "TECHNOLOGY"
    elseif a == "f"
        return "FUEL"
    elseif a == "e"
        return "EMISSION"
    elseif a == "m"
        return "MODE_OF_OPERATION"
    elseif a == "r" || a == "rr"
        return "REGION"
    elseif a == "s"
        return "STORAGE"
    elseif a == "l"
        return "TIMESLICE"
    elseif a == "tg1"
        return "TSGROUP1"
    elseif a == "tg2"
        return "TSGROUP2"
    elseif a == "ls"
        return "SEASON"
    elseif a == "ld"
        return "DAYTYPE"
    elseif a == "lh"
        return "DAILYTIMEBRACKET"
    elseif a == "n" || a == "n1" || a == "n2"
        return "NODE"
    elseif a == "tr"
        return "TransmissionLine"
    else
        return a
    end
end  # translatesetabb(a::String)

"""
    keydicts(df::DataFrames.DataFrame, numdicts::Int)

Generates `Dicts` that can be used to restrict JuMP constraints or variables to selected indices
(rather than all values in their dimensions) at creation. Returns an array of the generated `Dicts`.

# Arguments
- `df::DataFrames.DataFrame`: The results of a query that selects the index values.
- `numdicts::Int`: The number of `Dicts` that should be created.  The first `Dict` is for the first field in
        the query, the second is for the first two fields in the query, and so on.  `Dict` keys are arrays
        of the values of the key fields to which the `Dict` corresponds, and `Dict` values are `Sets` of corresponding
        values in the next field of the query (field #2 for the first `Dict`, field #3 for the second `Dict`, and so on).
"""
function keydicts(df::DataFrames.DataFrame, numdicts::Int)
    local returnval = Array{Dict{Array{String,1},Set{String}},1}()  # Function return value

    # Set up empty dictionaries in returnval
    for i in 1:numdicts
        push!(returnval, Dict{Array{String,1},Set{String}}())
    end

    # Populate dictionaries using df
    for row in DataFrames.eachrow(df)
        for j in 1:numdicts
            if !haskey(returnval[j], [row[k] for k = 1:j])
                returnval[j][[row[k] for k = 1:j]] = Set{String}()
            end

            push!(returnval[j][[row[k] for k = 1:j]], row[j+1])
        end
    end

    return returnval
end  # keydicts(df::DataFrames.DataFrame, numdicts::Int)

"""
    keydicts_parallel(df::DataFrames.DataFrame, numdicts::Int,
    targetprocs::Array{Int,1})

Runs `keydicts` using parallel processes identified in `targetprocs`. If there are not at
least 10,000 rows in `df` for each target process, resorts to running `keydicts` on process
that called `keydicts_parallel`.
"""
function keydicts_parallel(df::DataFrames.DataFrame, numdicts::Int, targetprocs::Array{Int,1})
    local returnval = Array{Dict{Array{String,1},Set{String}},1}()  # Function return value
    local availprocs::Array{Int, 1} = intersect(procs(), targetprocs)  # Targeted processes that actually exist
    local np::Int = length(availprocs)  # Number of targeted processes that actually exist
    local dfrows::Int = size(df)[1]  # Number of rows in df

    if np <= 1 || div(dfrows, np) < 10000
        # Run keydicts for entire df
        returnval = keydicts(df, numdicts)
    else
        # Divide operation among available processes
        local blockdivrem = divrem(dfrows, np)  # Quotient and remainder from dividing dfrows by np; element 1 = quotient, element 2 = remainder
        local results = Array{typeof(returnval), 1}(undef, np)  # Collection of results from async processing

        # Dispatch async tasks in main process, each of which performs a remotecall_fetch on an available process. Wrap in sync block
        #   to wait until all async processes finish before proceeding.
        @sync begin
            for p=1:np
                @async begin
                    # In each process, execute keydicts on a block of rows from df
                    results[p] = remotecall_fetch(NemoMod.keydicts, availprocs[p], df[((p-1) * blockdivrem[1] + 1):((p) * blockdivrem[1] + (p == np ? blockdivrem[2] : 0)),:], numdicts)
                end
            end
        end

        # Merge results from async tasks
        for i = 1:numdicts
            push!(returnval, Dict{Array{String,1},Set{String}}())

            for j = 1:np
                returnval[i] = merge(union, returnval[i], results[j][i])
            end
        end
    end

    return returnval
end  # keydicts_parallel(df::DataFrames.DataFrame, numdicts::Int, targetprocs::Array{Int,1})

#= Testing code:
using JuMP, SQLite, IterTools, NullableArrays, DataFrames
dbpath = "C:\\temp\\utopia_2015_08_27.sl3"
db = SQLite.DB(dbpath)
stechnology = dropnull(SQLite.query(db, "select val from technology limit 10")[:val])
@constraintref(abc[1:10])
m = Model()
@variable(m, vtest[stechnology])
createconstraint(m, "TestConstraint", abc, DataFrames.eachrow(SQLite.query(db, "select t.val as t from technology t limit 10")), ["t"], "vtest[t]", ">=", "0")
=#

# Demonstration function - performance is too poor for production use
#=
function createconstraint(model::JuMP.Model, logname::String, constraintref::Array{JuMP.ConstraintRef,1}, rows::DataFrames.DFRowIterator,
    keyfields::Array{String,1}, lh::String, operator::String, rh::String)

    local constraintnum = 1  # Number of next constraint to be added to constraint array
    global row = Void  # Row in subsequent iteration; needs to be defined as global to be accessible in eval statements

    for rw in rows
        row = rw

        for kf in keyfields
            eval(parse(kf * " = row[:" * kf * "]"))
        end

        if operator == "=="
            constraintref[constraintnum] = @constraint(model, eval(parse(lh)) == eval(parse(rh)))
        elseif operator == "<="
            constraintref[constraintnum] = @constraint(model, eval(parse(lh)) <= eval(parse(rh)))
        elseif operator == ">="
            constraintref[constraintnum] = @constraint(model, eval(parse(lh)) >= eval(parse(rh)))
        end

        constraintnum += 1
    end

    logmsg("Created constraint " * logname * ".")
end  # createconstraint
=#

"""
    savevarresults(vars::Array{String,1},
    modelvarindices::Dict{String, Tuple{AbstractArray,Array{String,1}}},
    db::SQLite.DB, solvedtmstr::String, reportzeros::Bool = false, quiet::Bool = false)

Saves model results to a SQLite database using SQL inserts with transaction batching.

# Arguments
- `vars::Array{String,1}`: Names of model variables for which results will be retrieved and saved to database.
- `modelvarindices::Dict{String, Tuple{AbstractArray,Array{String,1}}}`: Dictionary mapping model variable names
    to tuples of (variable, [index column names]).
- `db::SQLite.DB`: SQLite database.
- `solvedtmstr::String`: String to write into solvedtm field in result tables.
- `reportzeros::Bool`: Indicates whether values equal to 0 should be saved.
- `quiet::Bool = false`: Suppresses low-priority status messages (which are otherwise printed to STDOUT).
"""
function savevarresults(vars::Array{String,1}, modelvarindices::Dict{String, Tuple{AbstractArray,Array{String,1}}}, db::SQLite.DB, solvedtmstr::String,
    reportzeros::Bool = false, quiet::Bool = false)
    for vname in intersect(vars, keys(modelvarindices))
        local v = modelvarindices[vname][1]  # Model variable corresponding to vname (technically, a variable container in JuMP)
        local gvv = value.(v)  # Value of v in model solution
        local indices::Array{Tuple}  # Array of tuples of index values for v
        local vals::Array{Float64}  # Array of model results for v

        # Populate indices and vals
        if v isa JuMP.Containers.DenseAxisArray
            # In this case, indices and vals are parallel multidimensional arrays
            indices = collect(Base.product(gvv.axes...))  # Ellipsis splats axes into its values for passing to product()
            vals = gvv.data
        elseif v isa JuMP.Containers.SparseAxisArray
            # In this case, indices and vals are vectors
            indices = collect(keys(gvv.data))
            vals = collect(values(gvv.data))
        end

        # BEGIN: Wrap database operations in try-catch block to allow rollback on error.
        try
            # Begin SQLite transaction
            SQLite.DBInterface.execute(db, "BEGIN")

            # Create target table for v (drop any existing table for v)
            SQLite.DBInterface.execute(db,"drop table if exists " * vname)

            SQLite.DBInterface.execute(db, "create table '" * vname * "' ('" * join(modelvarindices[vname][2], "' text, '") * "' text, 'val' real, 'solvedtm' text)")

            # Insert data from gvv
            for i in eachindex(indices)
                if reportzeros || vals[i] != 0.0
                    SQLite.DBInterface.execute(db, "insert into " * vname * " values('" * join(indices[i], "', '") * "', '" * string(vals[i]) * "', '" * solvedtmstr * "')")
                end
            end

            # Complete SQLite transaction
            SQLite.DBInterface.execute(db, "COMMIT")
            logmsg("Saved results for " * vname * " to database.", quiet)
        catch
            # Rollback db transaction
            SQLite.DBInterface.execute(db, "ROLLBACK")

            # Proceed with normal Julia error sequence
            rethrow()
        end
        # END: Wrap database operations in try-catch block to allow rollback on error.
    end
end  # savevarresults(vars::Array{String,1}, modelvarindices::Dict{String, Tuple{AbstractArray,Array{String,1}}}, db::SQLite.DB, solvedtmstr::String)

"""
    checkactivityupperlimits(db::SQLite.DB, tolerance::Float64, restrictyears::Bool, inyears::String)

For the scenario specified in `db`, checks whether there are:
    1) Any TotalTechnologyAnnualActivityUpperLimit values that are <= tolerance times the maximum
        annual demand
    2) Any TotalTechnologyModelPeriodActivityUpperLimit values that are <= tolerance times the maximum
        annual demand times the number of years in the scenario
If `restrictyears` is true, the maximum annual demand and the TotalTechnologyAnnualActivityUpperLimit values
are selected from the years in `inyears`. `inyears` should be a comma-delimited list of years in
parentheses, with a space at the beginning and end of the string.

Returns a `Tuple` of Booleans whose first value is the result of the first check, and whose second
    value is the result of the second.
"""
function checkactivityupperlimits(db::SQLite.DB, tolerance::Float64, restrictyears::Bool, inyears::String)
    local annual::Bool = true  # Return value indicating whether there are any TotalTechnologyAnnualActivityUpperLimit
        # values that are <= tolerance x maximum annual demand
    local modelperiod::Bool = true  # Return value indicating whether there are any TotalTechnologyModelPeriodActivityUpperLimit
        # values that are <= tolerance x maximum annual demand x # of years in scenario
    local qry::DataFrame  # Working query
    local maxdemand::Float64  # Maximum annual demand

    qry = SQLite.DBInterface.execute(db, "select max(mv) as mv from
    (select max(val) as mv from AccumulatedAnnualDemand_def $(restrictyears ? "where y in" * inyears : "")
    union
    select max(val) as mv from SpecifiedAnnualDemand_def $(restrictyears ? "where y in" * inyears : ""))") |> DataFrame

    maxdemand = qry[!,1][1]

    qry = SQLite.DBInterface.execute(db, "select tau.val from TotalTechnologyAnnualActivityUpperLimit_def tau
        where tau.val / :v1 <= :v2 $(restrictyears ? "and y in" * inyears : "")", [maxdemand, tolerance]) |> DataFrame

    if size(qry)[1] == 0
        annual = false
    end

    qry = SQLite.DBInterface.execute(db, "select tmu.val from TotalTechnologyModelPeriodActivityUpperLimit_def tmu
      where tmu.val / (:v1 * (select count(val) from year)) <= :v2", [maxdemand, tolerance]) |> DataFrame

    if size(qry)[1] == 0
      modelperiod = false
    end

    return (annual, modelperiod)
end  # checkactivityupperlimits(db::SQLite.DB, tolerance::Float64, restrictyears::Bool, inyears::String)

"""
    scenario_calc_queries(dbpath::String, transmissionmodeling::Bool, vproductionbytechnologysaved::Bool,
        vusebytechnologysaved::Bool, restrictyears::Bool, inyears::String)

Returns a `Dict` of query commands used in NEMO's `modelscenario` function. Each key in the return value is
    a query name, and each value is a `Tuple` where:
        - Element 1 = path to NEMO scenario database in which to execute query (taken from `dbpath` argument)
        - Element 2 = query's SQL statement
    The function's arguments other than `dbpath` delimit the set of returned query commands as noted below.

# Arguments
- `dbpath::String`: Path to NEMO scenario database in which query commands should be executed.
- `transmissionmodeling::Bool`: Indicates whether transmission modeling is enabled in `modelscenario`.
    Additional query commands are included in results when transmission modeling is enabled.
- `vproductionbytechnologysaved::Bool`: Indicates whether output variable `vproductionbytechnology`
    will be saved in `modelscenario`. Additional query commands are included in results when this argument
    and `transmissionmodeling` are `true`.
- `vusebytechnologysaved::Bool`: Indicates whether output variable `vusebytechnology`
    will be saved in `modelscenario`. Additional query commands are included in results when this argument
    and `transmissionmodeling` are `true`.
- `restrictyears::Bool`: Indicates whether `modelscenario` is running for selected years only.
- `inyears::String`: SQL IN clause predicate for years selected for `modelscenario`. When `restrictvars`
    is `true`, this argument is used to include filtering by year in query commands."""
function scenario_calc_queries(dbpath::String, transmissionmodeling::Bool, vproductionbytechnologysaved::Bool,
    vusebytechnologysaved::Bool, restrictyears::Bool, inyears::String)

    return_val::Dict{String, Tuple{String, String}} = Dict{String, Tuple{String, String}}()  # Return value for this function; map of query names
    #   to tuples of (DB path, SQL command)

    return_val["queryvrateofactivityvar"] = (dbpath, "with ar as (select r, t, m, y from OutputActivityRatio_def
    where val <> 0 $(restrictyears ? "and y in" * inyears : "")
    union
    select r, t, m, y from InputActivityRatio_def
    where val <> 0 $(restrictyears ? "and y in" * inyears : ""))
    select r.val as r, l.val as l, t.val as t, m.val as m, y.val as y
    from REGION r, TIMESLICE l, TECHNOLOGY t, MODE_OF_OPERATION m, YEAR y, ar
    where ar.r = r.val and ar.t = t.val and ar.m = m.val and ar.y = y.val
    order by r.val, t.val, l.val, y.val")

    return_val["queryvtrade"] = (dbpath, "select r.val as r, rr.val as rr, l.val as l, f.val as f, y.val as y
    from region r, region rr, TIMESLICE l, FUEL f, year y, TradeRoute_def tr
    WHERE
    r.val = tr.r and rr.val = tr.rr and f.val = tr.f and y.val = tr.y
    and tr.r <> tr.rr and tr.val = 1 $(restrictyears ? "and tr.y in" * inyears : "")
    order by r.val, rr.val, f.val, y.val")

    return_val["queryvtradeannual"] = (dbpath, "select r.val as r, rr.val as rr, f.val as f, y.val as y
    from region r, region rr, FUEL f, year y, TradeRoute_def tr
    WHERE
    r.val = tr.r and rr.val = tr.rr and f.val = tr.f and y.val = tr.y
    and tr.r <> tr.rr and tr.val = 1 $(restrictyears ? "and tr.y in" * inyears : "")")

    return_val["queryvrateofproductionbytechnologybymodenn"] = (dbpath, "select r.val as r, ys.l as l, t.val as t, m.val as m, f.val as f, y.val as y,
    cast(oar.val as real) as oar
    from region r, YearSplit_def ys, technology t, MODE_OF_OPERATION m, fuel f, year y, OutputActivityRatio_def oar
    left join TransmissionModelingEnabled tme on tme.r = r.val and tme.f = f.val and tme.y = y.val
    where oar.r = r.val and oar.t = t.val and oar.f = f.val and oar.m = m.val and oar.y = y.val
    and oar.val <> 0
    and ys.y = y.val
    and tme.id is null
    $(restrictyears ? "and y.val in" * inyears : "")
    order by r.val, ys.l, t.val, f.val, y.val")

    return_val["queryvrateofproductionbytechnologynn"] = (dbpath, "select r.val as r, ys.l as l, t.val as t, f.val as f, y.val as y, cast(ys.val as real) as ys
    from region r, YearSplit_def ys, technology t, fuel f, year y,
    (select distinct r, t, f, y
    from OutputActivityRatio_def
    where val <> 0 $(restrictyears ? "and y in" * inyears : "")) oar
    left join TransmissionModelingEnabled tme on tme.r = r.val and tme.f = f.val and tme.y = y.val
    where oar.r = r.val and oar.t = t.val and oar.f = f.val and oar.y = y.val
    and ys.y = y.val
    and tme.id is null
    order by r.val, ys.l, f.val, y.val")

    return_val["queryvproductionbytechnologyannual"] = (dbpath, "select * from (
    select r.val as r, t.val as t, f.val as f, y.val as y, null as n, ys.l as l,
    cast(ys.val as real) as ys
    from region r, technology t, fuel f, year y, YearSplit_def ys,
    (select distinct r, t, f, y
    from OutputActivityRatio_def
    where val <> 0 $(restrictyears ? "and y in" * inyears : "")) oar
    left join TransmissionModelingEnabled tme on tme.r = r.val and tme.f = f.val and tme.y = y.val
    where oar.r = r.val and oar.t = t.val and oar.f = f.val and oar.y = y.val
    and ys.y = y.val
    and tme.id is null
    union all
    select n.r as r, ntc.t as t, oar.f as f, ntc.y as y, ntc.n as n, ys.l as l,
    cast(ys.val as real) as ys
    from NodalDistributionTechnologyCapacity_def ntc, NODE n,
    TransmissionModelingEnabled tme, YearSplit_def ys,
    (select distinct r, t, f, y
    from OutputActivityRatio_def
    where val <> 0 $(restrictyears ? "and y in" * inyears : "")) oar
    where ntc.val > 0
    and ntc.n = n.val
    and tme.r = n.r and tme.f = oar.f and tme.y = ntc.y
    and oar.r = n.r and oar.t = ntc.t and oar.y = ntc.y
    and ntc.y = ys.y
    )
    order by r, t, f, y")

    return_val["queryvrateofusebytechnologybymodenn"] = (dbpath, "select r.val as r, ys.l as l, t.val as t, m.val as m, f.val as f, y.val as y, cast(iar.val as real) as iar
    from region r, YearSplit_def ys, technology t, MODE_OF_OPERATION m, fuel f, year y, InputActivityRatio_def iar
    left join TransmissionModelingEnabled tme on tme.r = r.val and tme.f = f.val and tme.y = y.val
    where iar.r = r.val and iar.t = t.val and iar.f = f.val and iar.m = m.val and iar.y = y.val
    and iar.val <> 0
    and ys.y = y.val
    and tme.id is null
    $(restrictyears ? "and y.val in" * inyears : "")
    order by r.val, ys.l, t.val, f.val, y.val")

    return_val["queryvrateofusebytechnologynn"] = (dbpath, "select
    r.val as r, ys.l as l, t.val as t, f.val as f, y.val as y, cast(ys.val as real) as ys
    from region r, YearSplit_def ys, technology t, fuel f, year y,
    (select distinct r, t, f, y from InputActivityRatio_def where val <> 0 $(restrictyears ? "and y in" * inyears : "")) iar
    left join TransmissionModelingEnabled tme on tme.r = r.val and tme.f = f.val and tme.y = y.val
    where iar.r = r.val and iar.t = t.val and iar.f = f.val and iar.y = y.val
    and ys.y = y.val
    and tme.id is null
    order by r.val, ys.l, f.val, y.val")

    return_val["queryvusebytechnologyannual"] = (dbpath, "select * from (
    select r.val as r, t.val as t, f.val as f, y.val as y, null as n, ys.l as l,
    cast(ys.val as real) as ys
    from region r, technology t, fuel f, year y, YearSplit_def ys,
    (select distinct r, t, f, y from InputActivityRatio_def where val <> 0 $(restrictyears ? "and y in" * inyears : "")) iar
    left join TransmissionModelingEnabled tme on tme.r = r.val and tme.f = f.val and tme.y = y.val
    where iar.r = r.val and iar.t = t.val and iar.f = f.val and iar.y = y.val
    and ys.y = y.val
    and tme.id is null
    union all
    select n.r as r, ntc.t as t, iar.f as f, ntc.y as y, ntc.n as n, ys.l as l,
    cast(ys.val as real) as ys
    from NodalDistributionTechnologyCapacity_def ntc, NODE n,
    TransmissionModelingEnabled tme, YearSplit_def ys,
    (select distinct r, t, f, y from InputActivityRatio_def where val <> 0 $(restrictyears ? "and y in" * inyears : "")) iar
    where ntc.val > 0
    and ntc.n = n.val
    and tme.r = n.r and tme.f = iar.f and tme.y = ntc.y
    and iar.r = n.r and iar.t = ntc.t and iar.y = ntc.y
    and ntc.y = ys.y
    )
    order by r, t, f, y")

    return_val["querycaa5_totalnewcapacity"] = (dbpath, "select cot.r as r, cot.t as t, cot.y as y, cast(cot.val as real) as cot
    from CapacityOfOneTechnologyUnit_def cot where cot.val <> 0 $(restrictyears ? "and cot.y in" * inyears : "")")

    if transmissionmodeling
        return_val["queryvrateofactivitynodal"] = (dbpath, "select ntc.n as n, l.val as l, ntc.t as t, ar.m as m, ntc.y as y
        from NodalDistributionTechnologyCapacity_def ntc, node n,
        	TransmissionModelingEnabled tme, TIMESLICE l,
        (select r, t, f, m, y from OutputActivityRatio_def
        where val <> 0 $(restrictyears ? "and y in" * inyears : "")
        union
        select r, t, f, m, y from InputActivityRatio_def
        where val <> 0 $(restrictyears ? "and y in" * inyears : "")) ar
        where ntc.val > 0
        and ntc.n = n.val
        and tme.r = n.r and tme.f = ar.f and tme.y = ntc.y
        and ar.r = n.r and ar.t = ntc.t and ar.y = ntc.y
        order by ntc.n, ntc.t, l.val, ntc.y")

        return_val["queryvrateofproductionbytechnologynodal"] = (dbpath, "select ntc.n as n, ys.l as l, ntc.t as t, oar.f as f, ntc.y as y,
        	cast(ys.val as real) as ys
        from NodalDistributionTechnologyCapacity_def ntc, YearSplit_def ys, NODE n,
    	TransmissionModelingEnabled tme,
        (select distinct r, t, f, y
        from OutputActivityRatio_def
        where val <> 0 $(restrictyears ? "and y in" * inyears : "")) oar
        where ntc.val > 0
        and ntc.y = ys.y
        and ntc.n = n.val
        and tme.r = n.r and tme.f = oar.f and tme.y = ntc.y
    	and oar.r = n.r and oar.t = ntc.t and oar.y = ntc.y
        order by ntc.n, ys.l, oar.f, ntc.y")

        return_val["queryvrateofusebytechnologynodal"] = (dbpath, "select ntc.n as n, ys.l as l, ntc.t as t, iar.f as f, ntc.y as y,
        	cast(ys.val as real) as ys
        from NodalDistributionTechnologyCapacity_def ntc, YearSplit_def ys, NODE n,
    	TransmissionModelingEnabled tme,
        (select distinct r, t, f, y
        from InputActivityRatio_def
        where val <> 0 $(restrictyears ? "and y in" * inyears : "")) iar
        where ntc.val > 0
        and ntc.y = ys.y
        and ntc.n = n.val
        and tme.r = n.r and tme.f = iar.f and tme.y = ntc.y
    	and iar.r = n.r and iar.t = ntc.t and iar.y = ntc.y
        order by ntc.n, ys.l, iar.f, ntc.y")

        if vproductionbytechnologysaved
            return_val["queryvproductionbytechnologyindices_nodalpart"] = (dbpath, "select distinct n.r as r, ys.l as l, ntc.t as t, oar.f as f, ntc.y as y, null as ys
            from NodalDistributionTechnologyCapacity_def ntc, YearSplit_def ys, NODE n,
            TransmissionModelingEnabled tme,
            (select distinct r, t, f, y
            from OutputActivityRatio_def
            where val <> 0 $(restrictyears ? "and y in" * inyears : "")) oar
            where ntc.val > 0
            and ntc.y = ys.y
            and ntc.n = n.val
            and tme.r = n.r and tme.f = oar.f and tme.y = ntc.y
            and oar.r = n.r and oar.t = ntc.t and oar.y = ntc.y")
        end

        if vusebytechnologysaved
            return_val["queryvusebytechnologyindices_nodalpart"] = (dbpath, "select distinct n.r as r, ys.l as l, ntc.t as t, iar.f as f, ntc.y as y, null as ys
            from NodalDistributionTechnologyCapacity_def ntc, YearSplit_def ys, NODE n,
            TransmissionModelingEnabled tme,
            (select distinct r, t, f, y
            from InputActivityRatio_def
            where val <> 0 $(restrictyears ? "and y in" * inyears : "")) iar
            where ntc.val > 0
            and ntc.y = ys.y
            and ntc.n = n.val
            and tme.r = n.r and tme.f = iar.f and tme.y = ntc.y
            and iar.r = n.r and iar.t = ntc.t and iar.y = ntc.y")
        end

        return_val["queryvtransmissionbyline"] = (dbpath, "select tl.id as tr, ys.l as l, tl.f as f, tme1.y as y, tl.n1 as n1, tl.n2 as n2,
    	tl.reactance as reactance, tme1.type as type, tl.maxflow as maxflow,
        cast(tl.VariableCost as real) as vc, cast(ys.val as real) as ys,
        cast(tl.fixedcost as real) as fc, cast(tcta.val as real) as tcta
        from TransmissionLine tl, NODE n1, NODE n2, TransmissionModelingEnabled tme1,
        TransmissionModelingEnabled tme2, YearSplit_def ys, TransmissionCapacityToActivityUnit_def tcta
        where
        tl.n1 = n1.val and tl.n2 = n2.val
        and tme1.r = n1.r and tme1.f = tl.f
        and tme2.r = n2.r and tme2.f = tl.f
        and tme1.y = tme2.y and tme1.type = tme2.type
    	and ys.y = tme1.y $(restrictyears ? "and ys.y in" * inyears : "")
    	and tcta.r = n1.r and tl.f = tcta.f
        order by tl.id, tme1.y")

        return_val["queryvstorageleveltsgroup1"] = (dbpath, "select ns.n as n, ns.s as s, tg1.name as tg1, ns.y as y
        from nodalstorage ns, TSGROUP1 tg1 $(restrictyears ? "where ns.y in" * inyears : "")")

        return_val["queryvstorageleveltsgroup2"] = (dbpath, "select ns.n as n, ns.s as s, tg1.name as tg1, tg2.name as tg2, ns.y as y
        from nodalstorage ns, TSGROUP1 tg1, TSGROUP2 tg2 $(restrictyears ? "where ns.y in" * inyears : "")")

        return_val["queryvstoragelevelts"] = (dbpath, "select ns.n as n, ns.s as s, l.val as l, ns.y as y
        from nodalstorage ns, TIMESLICE l $(restrictyears ? "where ns.y in" * inyears : "")")
    end  # transmissionmodeling

    return return_val
end  # scenario_calc_queries()

"""
    run_qry(qtpl::Tuple{String, String})

Runs the SQLite database query specified in `qtpl` and returns the result as a DataFrame.
    Element 1 in `qtpl` should be the path to the SQLite database, and element 2 should be the query's
    SQL command. Designed to work with the output of `scenario_calc_queries` in a `map` or `pmap` call."""
function run_qry(qtpl::Tuple{String, String})
    return SQLite.DBInterface.execute(SQLite.DB(qtpl[1]), qtpl[2]) |> DataFrame
end  # run_qry(qtpl::Tuple{String, String, String})

"""
    writescenariomodel(dbpath::String;
        calcyears::Array{Int, 1} = Array{Int, 1}(),
        varstosave::String = "vdemandnn, vnewcapacity, vtotalcapacityannual,
            vproductionbytechnologyannual, vproductionnn, vusebytechnologyannual,
            vusenn, vtotaldiscountedcost",
        numprocs::Int = 0, targetprocs::Array{Int, 1} = Array{Int, 1}(),
        restrictvars::Bool = true, continuoustransmission::Bool = false,
        forcemip::Bool = false, quiet::Bool = false,
        writefilename::String = "nemomodel.bz2",
        writefileformat::MathOptInterface.FileFormats.FileFormat = MathOptInterface.FileFormats.FORMAT_MPS
    )

Writes a file representing the optimization problem for a NEMO scenario. Returns the name of the file.
All arguments except `writefilename` and `writefileformat` function as in
[`calculatescenario`](@ref). `writefilename` and `writefileformat` are described below.

# Arguments

- `writefilename::String`: Name of the output file. If a path is not included in the name,
    the file is written to the Julia working directory. If the name ends in `.gz`, the file
    is compressed with Gzip. If the name ends in `.bz2`, the file is compressed with BZip2.
- `writefileformat::MathOptInterface.FileFormats.FileFormat`: Data format used in the output
    file. Common formats include `MathOptInterface.FileFormats.FORMAT_MPS` (MPS format) and
    `MathOptInterface.FileFormats.FORMAT_LP` (LP format). See the documentation for
    [`MathOptInterface`](https://github.com/jump-dev/MathOptInterface.jl) for additional options.
"""
function writescenariomodel(
    dbpath::String;
    calcyears::Array{Int, 1} = Array{Int, 1}(),
    varstosave::String = "vdemandnn, vnewcapacity, vtotalcapacityannual, vproductionbytechnologyannual, vproductionnn, vusebytechnologyannual, vusenn, vtotaldiscountedcost",
    numprocs::Int = 0,
    targetprocs::Array{Int, 1} = Array{Int, 1}(),
    restrictvars::Bool = true,
    continuoustransmission::Bool = false,
    forcemip = false,
    quiet::Bool = false,
    writefilename::String = "nemomodel.bz2",
    writefileformat::MathOptInterface.FileFormats.FileFormat = MathOptInterface.FileFormats.FORMAT_MPS
    )

    try
        modelscenario(dbpath; calcyears=calcyears, varstosave=varstosave, numprocs=numprocs, targetprocs=targetprocs, restrictvars=restrictvars,
            continuoustransmission=continuoustransmission, forcemip=forcemip, quiet=quiet, writemodel=true,
            writefilename=writefilename, writefileformat=writefileformat)
    catch e
        println("NEMO encountered an error with the following message: " * sprint(showerror, e) * ".")
        println("To report this issue to the NEMO team, please submit an error report at https://leap.sei.org/support/. Please include in the report a list of steps to reproduce the error and the error message.")
    end
end  # writescenariomodel


text:
# --------------------------------  Dynamic System State --------------------------------
"""
    DSState{T}
    DSState(x::AbstractVector)

Dynamic system state of type `T` (type of the elements of the state vector, ie. `Float64`).
"""
abstract type DSState{T} end

eltype(::Type{<:DSState{T}}) where {T} = T

"""
    get_x(dss::DSState)
    get_x(dssd::DSStateDot)

Return the state vector `x` of the dynamic system.
"""
get_x(dss::DSState) = dss.x

"""
    get_n_states(dss::DSState)
    get_n_states(dssd::DSStateDot)

Return the length of the state vector `x` of the dynamic system.
"""
get_n_states(dss::DSState) = length(dss.x)

# ----------------------- Methods to be implemented by each subtype -----------------------
"""
    get_x_names(dss::DSState)
    get_x_names(dssd::DSStateDot)

Return an array of symbols with the name for each variable in the state vector.
"""
function get_x_names end

"""
    state_eqs(dss::DSState, time, mass, inertia, forces, moments, h)

Calculate the derivative of the state vector using the state equations.

# Arguments

- `dss::DSState`: dynamic system state.
- `time::Number`: time (s).
- `mass::Number`: mass (kg).
- `inertia::AbstractMatrix`: 3x3 inertia tensor (kg·m^2).
- `forces::AbstractVector`: body axis forces [Fx, Fy, Fz] (N).
- `moments::AbstractVector`: body axis moments (L, M, N) (N·m).
- `h::AbstractVector`:  Additional angular momentum contributions such as those coming from
 spinning rotors (kg·m²/s).

 # Returns
 - `dssd::DSStateDot`: dynamic system state dot.
"""
function state_eqs end

"""
    get_earth_positon(dss::DSState)
    get_earth_positon(dssd::DSStateDot)

Return Earth position [x, y, z] (m).
"""
function get_earth_position end

"""
    get_height(dss::DSState)
    get_height(dssd::DSStateDot)

Return height (m).
"""
function get_height end

"""
    get_euler_angles(dss::DSState)
    get_euler_angles(dssd::DSStateDot)

Return Euler angles (ψ, θ, ϕ) (rad).
"""
function get_euler_angles end

"""
    get_tas(Dss::DState)
    get_tas(dssd::DSStateDot)

Return true air speed (m/s).
"""
function get_tas end

"""
    get_α(dss::DSState)
    get_α(dssd::DSStateDot)

Get angle of attack (rad).
"""
function get_α end

"""
    get_β(dss::DSState)
    get_β(dssd::DSStateDot)

Get angle of side slip (rad).
"""
function get_β end

"""
    get_tasαβ(dss:DSState)
    get_tasαβ(dssd::DSStateDot)

Get [tas (m/s), α (rad), β (rad)].
"""
function get_tasαβ end

"""
    get_body_velocity(dss::DState)
    get_body_velocity(dssd::DSStateDot)

Get velocity in body axis [u, v, w] (m/s).
"""
function get_body_velocity end

"""
    get_horizon_velocity(dss::DSState)
    get_horizon_velocity(dssd::DSStateDot)

Get velocity in horizon axis [Vn, Ve, Vd] (m/s).
"""
function get_horizon_velocity end

"""
    get_ang_vel_body(dss::DSState)
    get_ang_vel_body(dssd::DSStateDot)

Get angular velocity in body axis [p, q, r] (rad/s).
"""
function get_ang_vel_body end

"""
    get_euler_angles_rates(dss::DSState)
    get_euler_angles_rates(dssd::DSStateDot)

Get euler angles rates [ψ_dot, θ_dot, ϕ_dot] (rad/s).
"""
function get_euler_angles_rates end

"""
    get_engine_power(dss:DSState)
    get_engine_power(dssd::DSStateDot)

Get engine power (%).
"""
function get_engine_power end

# ------------------------------- Dynamic System State  Dot ------------------------------
"""
    DSStateDot{S, N, T}
    DSStateDot(dss::DSState, xdot::AbstractArray)

Time derivative, `x_dot`, of the state vector `x`.

- `S`: type of the associated `DSState`.
- `N`: number of elements of the `x_dot` vector.
- `T`: type of the elements of `x_dot` vector.
"""
struct DSStateDot{S, N, T}
    dss::S
    xdot::SVector{N, T}
    DSStateDot(dss::S, xdot::V) where {S<:DSState, V<:SVector} =
        new{typeof(dss), get_n_states(dss), eltype(dss)}(
            dss,
            SVector{get_n_states(dss), eltype(dss)}(xdot)
        )
end

function DSStateDot(dss::DSState, xdot::AbstractArray)
    return DSStateDot(dss, SVector{get_n_states(dss), eltype(dss)}(xdot))
end

"""
    get_xdot(dssd::DSStateDot)

Get the time derivative of state vector, `x_dot`.
"""
get_xdot(dssd::DSStateDot) = dssd.xdot

"""
    get_ds_state(dssd::DSStateDot)

Get the `DSState`.
"""
get_ds_state(dssd::DSStateDot) = dssd.dss

# --------------------- Methoos reexported from coposition with DSState --------------------
get_x(dssd::DSStateDot) = get_x(dssd.dss)
get_x_names(dssd::DSStateDot) = get_x_names(dssd.dss)
get_n_states(dssd::DSStateDot) = get_n_states(dssd.dss)
get_earth_position(dssd::DSStateDot) = get_earth_position(dssd.dss)
get_height(dssd::DSStateDot) = get_height(dssd.dss)
get_euler_angles(dssd::DSStateDot) = get_euler_angles(dssd.dss)
get_tas(dssd::DSStateDot) = get_tas(dssd.dss)
get_α(dssd::DSStateDot) = get_α(dssd.dss)
get_β(dssd::DSStateDot) = get_β(dssd.dss)
get_tasαβ(dssd::DSStateDot) = get_tasαβ(dssd.dss)
get_body_velocity(dssd::DSStateDot) = get_body_velocity(dssd.dss)
get_horizon_velocity(dssd::DSStateDot) = get_horizon_velocity(dssd.dss)
get_ang_vel_body(dssd::DSStateDot) = get_ang_vel_body(dssd.dss)
get_euler_angles_rates(dssd::DSStateDot) = get_euler_angles_rates(dssd.dss)
get_engine_power(dssd::DSStateDot) = get_engine_power(dssd.dss)

# ----------------------- Methods to be implemented by each subtype -----------------------
"""
    get_tas_dot(dssd::DSStateDot)

Get derivative of true air speed (m/s²).
"""
function get_tas_dot end

"""
    get_α_dot(dssd::DSStateDot)

Get angle of attack derivative (rad/s).
"""
function get_α_dot end

"""
    get_β_dot(dssd::DSStateDot)

Get derivative of angles of sideslip (rad/s).
"""
function get_β_dot end

"""
    get_tasαβ_dot(dssd::DSStateDot)

Get [tas_dot (m/s²), α_dot (rad/s), β_dot (rad/s)].
"""
function get_tasαβ_dot end

"""
    get_accel_body(dssd::DSStateDot)

Get accelerations in body axis [u_dot, v_dot, w_dot] (m/s²).
"""
function get_accel_body end

"""
    get_ang_accel_body(dssd::DSStateDot)

Get angular acceleration in body axis [p_dot, q_dot, r_dot] (rad/s²).
"""
function get_ang_accel_body end

"""
    get_engine_power_dot(dssd::DSStateDot)

Get derivative of engine power (s⁻¹).
"""
function get_engine_power_dot end

"""
    get_trimmer_cost(dssd::DSStateDot)

Get the cost function evaluation for the given `DSStateDot`.
"""
function get_trimmer_cost end

ext: jl


text:
module ThompsonSampling

export ThompsonSampler

include("thompson.jl")

end # module


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "GPI2"
version = v"1.5.1"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/cc-hpc-itwm/GPI-2/archive/refs/tags/v$(version).tar.gz", "4dac7e9152694d2ec4aefd982a52ecc064a8cb8f2c9eab0425428127c3719e2e")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd GPI-2-*
./autogen.sh
./configure --prefix=$prefix --build=${MACHTYPE} --host=${target} --without-fortran --with-ethernet
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("x86_64", "linux"; libc = "glibc"),
    Platform("aarch64", "linux"; libc = "glibc"),
    Platform("powerpc64le", "linux"; libc = "glibc")
]

# The products that we will ensure are always built
products = [
    LibraryProduct("libGPI2-stats", :libGPI2_stats),
    LibraryProduct("libGPI2-dbg", :libGPI2_dbg),
    LibraryProduct("libGPI2", :libGPI2),
    ExecutableProduct("gaspi_logger", :gaspi_logger),
    # The following are shell scripts that are required for GPI-2 to work,
    # thus they are added as `FileProduct`s
    FileProduct("bin/gaspi_run", :gaspi_run),
    FileProduct("bin/gaspi_cleanup", :gaspi_cleanup),
    FileProduct("bin/ssh.spawner", :ssh_spawner),
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")


text:
@testset "Caccessor + Cconst" begin
	@cstruct CconstStruct {
		i::Cint
		ci::Cconst{Cint}

		s::@cstruct {
			i::Cint
			ci::Cconst{Cint}
		}

		cs::Cconst{@cstruct {
			i::Cint
			ci::Cconst{Cint}
		}}

		u::@cunion {
			i::Cint
			ci::Cconst{Cint}
		}

		cu::Cconst{@cunion {
			i::Cint
			ci::Cconst{Cint}
		}}

		(as::{}[1])::@cstruct {
			i::Cint
			ci::Cconst{Cint}
		}

		(acs::{}[1])::Cconst{@cstruct {
			i::Cint
			ci::Cconst{Cint}
		}}

		@cstruct {
			si::Cint
			sci::Cconst{Cint}
		}

		Cconst{@cstruct {
			csi::Cint
			csci::Cconst{Cint}
		}}

		@cunion {
			ui::Cint
			uci::Cconst{Cint}
		}

		Cconst{@cunion {
			cui::Cint
			cuci::Cconst{Cint}
		}}
	}
	@test sizeof(CconstStruct) == sizeof(Cint)*18

	x = CconstStruct()
	@test typeof(x.i) === Cint
	@test typeof(x.ci) === Cint
	@test typeof(x.s) <: Caccessor{<:Cstruct}
	@test typeof(x.s.i) === Cint
	@test typeof(x.s.ci) === Cint
	@test typeof(x.cs) <: Caccessor{<:Cconst{<:Cstruct}}
	@test typeof(x.cs.i) === Cint
	@test typeof(x.cs.ci) === Cint
	@test typeof(x.u) <: Caccessor{<:Cunion}
	@test typeof(x.u.i) === Cint
	@test typeof(x.u.ci) === Cint
	@test typeof(x.cu) <: Caccessor{<:Cconst{<:Cunion}}
	@test typeof(x.cu.i) === Cint
	@test typeof(x.cu.ci) === Cint
	@test typeof(x.as) <: Caccessor{<:Carray{<:Cstruct}}
	@test typeof(x.as[1]) <: Caccessor{<:Cstruct}
	@test typeof(x.as[1].i) === Cint
	@test typeof(x.as[1].ci) === Cint
	@test typeof(x.acs) <: Caccessor{<:Carray{<:Cconst{<:Cstruct}}}
	@test typeof(x.acs[1]) <: Caccessor{<:Cconst{<:Cstruct}}
	@test typeof(x.acs[1].i) === Cint
	@test typeof(x.acs[1].ci) === Cint
	@test typeof(x.si) === Cint
	@test typeof(x.sci) === Cint
	@test typeof(x.csi) === Cint
	@test typeof(x.csci) === Cint
	@test typeof(x.ui) === Cint
	@test typeof(x.uci) === Cint
	@test typeof(x.cui) === Cint
	@test typeof(x.cuci) === Cint

	f(x) =
		x.i + x.ci +
		x.s.i + x.s.ci +
		x.cs.i + x.cs.ci +
		x.u.i + x.u.ci +
		x.cu.i + x.cu.ci +
		# x.as[1].i + x.as[1].ci +   # TODO: fix these to not allocate
		# x.acs[1].i + x.acs[1].ci +
		x.si + x.sci +
		x.csi + x.csci +
		x.ui + x.uci +
		x.cui + x.cuci
	@test f(x) == 0
	@test @allocated(f(x)) == 0

	@test x.i == 0
	x.i = 1
	@test x.i == 1

	@test x.ci == 0
	@test_throws ErrorException x.ci = 1
	@test x.ci == 0

	@test x.s[] == typeof(x.s[])()
	@test x.s.i == 0
	@test x.s.ci == 0
	x.s = typeof(x.s[])(i = 2, ci = 3)

	@test x.s.i == 2
	x.s.i = 22
	@test x.s.i == 22
	x.s[].i = 222
	@test x.s.i == 22

	@test x.s.ci == 3
	@test_throws ErrorException x.s.ci = 33
	@test x.s.ci == 3
	@test_throws ErrorException x.s[].ci = 333
	@test x.s.ci == 3

	@test x.cs[] == typeof(x.cs[])()
	@test x.cs.i == 0
	@test x.cs.ci == 0
	@test_throws ErrorException x.cs = typeof(x.cs[])(i = 4, ci = 5)

	@test x.cs.i == 0
	@test_throws ErrorException x.cs.i = 44
	@test x.cs.i == 0
	@test_throws ErrorException x.cs[].i = 444
	@test x.cs.i == 0

	@test x.cs.ci == 0
	@test_throws ErrorException x.cs.ci = 55
	@test x.cs.ci == 0
	@test_throws ErrorException x.cs[].ci = 555
	@test x.cs.ci == 0

	@test x.u[] == typeof(x.u[])()

	@test x.u.i == 0
	x.u = typeof(x.u[])(i = 6)
	@test x.u.i == 6
	x.u.i = 66
	@test x.u.i == 66
	x.u[].i = 666
	@test x.u.i == 66

	@test x.u.ci == 66
	x.u = typeof(x.u[])(ci = 7)
	@test x.u.ci == 7
	@test_throws ErrorException x.u.ci = 77
	@test x.u.ci == 7
	@test_throws ErrorException x.u[].ci = 777
	@test x.u.ci == 7

	@test x.cu[] == typeof(x.cu[])()

	@test x.cu.i == 0
	@test_throws ErrorException x.cu = typeof(x.cu[])(i = 8)
	@test x.cu.i == 0
	@test_throws ErrorException x.cu.i = 88
	@test x.cu.i == 0
	@test_throws ErrorException x.cu[].i = 888
	@test x.cu.i == 0

	@test x.cu.ci == 0
	@test_throws ErrorException x.cu = typeof(x.cu[])(ci = 9)
	@test x.cu.ci == 0
	@test_throws ErrorException x.cu.ci = 99
	@test x.cu.ci == 0
	@test_throws ErrorException x.cu[].ci = 999
	@test x.cu.ci == 0

	@test x.si == 0
	x.si = 1
	@test x.si == 1

	@test x.sci == 0
	@test_throws ErrorException x.sci = 2
	@test x.sci == 0

	@test x.csi == 0
	@test_throws ErrorException x.csi = 3
	@test x.csi == 0

	@test x.csci == 0
	@test_throws ErrorException x.csci = 4
	@test x.csci == 0

	@test x.ui == 0
	x.ui = 5
	@test x.ui == 5

	@test x.uci == 5
	@test_throws ErrorException x.uci = 6
	@test x.uci == 5

	@test x.cui == 0
	@test_throws ErrorException x.cui = 7
	@test x.cui == 0

	@test x.csci == 0
	@test_throws ErrorException x.csci = 8
	@test x.csci == 0

	@test f(x) == 51

	x = Cconst(CconstStruct())
	@test typeof(x.i) === Cint
	@test typeof(x.ci) === Cint
	@test typeof(x.s) <: Caccessor{<:Cconst{<:Cstruct}}
	@test typeof(x.s.i) === Cint
	@test typeof(x.s.ci) === Cint
	@test typeof(x.cs) <: Caccessor{<:Cconst{<:Cstruct}}
	@test typeof(x.cs.i) === Cint
	@test typeof(x.cs.ci) === Cint
	@test typeof(x.u) <: Caccessor{<:Cconst{<:Cunion}}
	@test typeof(x.u.i) === Cint
	@test typeof(x.u.ci) === Cint
	@test typeof(x.cu) <: Caccessor{<:Cconst{<:Cunion}}
	@test typeof(x.cu.i) === Cint
	@test typeof(x.cu.ci) === Cint
	@test typeof(x.as) <: Caccessor{<:Carray{<:Cconst{<:Cstruct}}}
	@test typeof(x.as[1]) <: Caccessor{<:Cconst{<:Cstruct}}
	@test typeof(x.as[1].i) === Cint
	@test typeof(x.as[1].ci) === Cint
	@test typeof(x.acs) <: Caccessor{<:Carray{<:Cconst{<:Cstruct}}}
	@test typeof(x.acs[1]) <: Caccessor{<:Cconst{<:Cstruct}}
	@test typeof(x.acs[1].i) === Cint
	@test typeof(x.acs[1].ci) === Cint
	@test typeof(x.si) === Cint
	@test typeof(x.sci) === Cint
	@test typeof(x.csi) === Cint
	@test typeof(x.csci) === Cint
	@test typeof(x.ui) === Cint
	@test typeof(x.uci) === Cint
	@test typeof(x.cui) === Cint
	@test typeof(x.cuci) === Cint

	@test f(x) == 0
	@test @allocated(f(x)) == 0

	@test x.i == 0
	@test_throws ErrorException x.i = 1
	@test x.i == 0

	@test x.s.i == 0
	@test_throws ErrorException x.s.i = 2
	@test x.s.i == 0

	@test x.u.i == 0
	@test_throws ErrorException x.u.i = 3
	@test x.u.i == 0

	@test x.si == 0
	@test_throws ErrorException x.si = 4
	@test x.si == 0

	@test x.ui == 0
	@test_throws ErrorException x.ui = 5
	@test x.ui == 0

	@test f(x) == 0
	@test @allocated(f(x)) == 0
end


text:
# Return x, or default value if it's nothing
defaultTo(x, default) = x ≠ nothing ? x : default

# Clear collection (seems to be pretty efficient)
clear(collection) = filter(_ -> false, collection)
clear!(collection) = filter!(_ -> false, collection)

# To string + add brackets if it's an expression
wrap(x) = string(x)
wrap(ex::Expr) = "($(string(ex)))"

# Map expression (skips first argument for CALL expressions)
walk(f, leaf) = f(leaf)
function walk(f, ex::Expr)
    new_ex = copy(ex)
    if ex.head == :call
        op, args = new_ex.args[1], new_ex.args[2:end]
        new_ex.args = [op, [walk(f, x) for x ∈ args]...]
    else
        println(ex)
        map!(x -> walk(f, x), new_ex.args)
    end
    new_ex
end

# Find an element in expression
contains(x, ex) = x == ex
function contains(x, ex::Expr)
    found = false
    walk(y -> found |= y==x, ex)
    found
end

# Split a collection by predicate
function partition(xs, pred)
    a, b = [], []
    for x ∈ xs
        push!(pred(x) ? a : b, x)
    end
    a, b
end

# Flatten nested array
flatten(x) = x
flatten(xs::Array) = vcat([flatten(x) for x ∈ xs]...)


text:
using Pkg

# source: https://discourse.julialang.org/t/how-to-find-out-the-version-of-a-package-from-its-module/37755/10?u=findmyway
pkgproject(m::Module) = Pkg.Operations.read_project(Pkg.Types.projectfile_path(pkgdir(m)))
pkgversion(m::Module) = pkgproject(m).version

struct InstrumentationInfo
    name::String
    version::VersionNumber
end

InstrumentationInfo(name::Module) =  InstrumentationInfo(name, pkgversion(name))
InstrumentationInfo(name::String) =  InstrumentationInfo(name, v"0.0.0-dev")


text:
using ParameterizedFunctions, DiffEqBase
using Test, LinearAlgebra, InteractiveUtils

using SpecialFunctions

### ODE Macros

println("Build some examples")
f_t = @ode_def SymCheck begin # Checks for error due to symbol on 1
  dx = x
  dy = -c*y + d*x*y*t^2
end a b c d

@test DiffEqBase.__has_syms(f_t)

f_t2 = @ode_def SymCheck2 begin # Checks for error due to symbol on 1
  dx = 1
  dy = -c*y + d*x*y*t^2
end a b c d

f_t3 = @ode_def ExprCheck begin # Checks for error due to symbol on 1
  dx = a*x - b*x*y
  dy = -c*y + d*x*y
end a b c d # Change to π after unicode fix

f = @ode_def_all LotkaVolterra begin
  dx = a*x - b*x*y
  dy = -c*y + d*x*y
end a b c d

f_2 = @ode_def LotkaVolterra3 begin
  dx = a*x - b^2*x*y
  dy = -c*y + d*x*y
end a b c d

println("Test Values")
t = 1.0
u = [2.0,3.0]
p = [1.5,1,3,1]
du = zeros(2)
grad = similar(du)
J = zeros(2,2)
iJ= zeros(2,2)
iW= zeros(2,2)
f(du,u,p,t)
@test du == [-3.0,-3.0]
@test du == f(u,p,t)
f_t(du,u,p,t)
@test du == [2.0,-3.0]
f_t2(du,u,p,t)
@test du == [1.0,-3.0]

println("Test t-gradient")
f.tgrad(grad,u,p,t)
@test grad == zeros(2)
f_t.tgrad(grad,u,p,t)
@test grad == [0.0;12.0]

println("Test Jacobians")
f.jac(J,u,p,t)
@test J  == [-1.5 -2.0
             3.0 -1.0]
#@test f.jac(u, p, t) == [-1.5 -2.0; 3.0 -1.0]

println("Test Inv Rosenbrock-W")
f.invW(iW,u,p,2.0,t)
@test minimum(iW - inv(I - 2*J) .< 1e-10)

f.invW_t(iW,u,p,2.0,t)
@test minimum(iW - inv(I/2 - J) .< 1e-10)

println("Parameter Jacobians")
pJ = zeros(Float64,2,4)
f.paramjac(pJ,u,[2.0;2.5;3.0;1.0],t)
@test pJ == [2.0 -6.0  0.0 0.0
             0.0  0.0 -3.0 6.0]

@code_llvm DiffEqBase.has_jac(f)

println("Test booleans")
@test DiffEqBase.has_jac(f) == true
@test DiffEqBase.has_paramjac(f) == true

@code_llvm DiffEqBase.has_paramjac(f)

println("Test difficult differentiable")
NJ = @ode_def DiffDiff begin
  dx = a*x - b*x*y
  dy = -c*y + erf(x*y/d)
end a b c d
NJ(du,u,[1.5,1,3,4],t)
@test du == [-3.0;-3*3.0 + erf(2.0*3.0/4)]
@test du == NJ(u, [1.5,1,3,4], t)
# NJ(Val{:jac},t,u,J) # Currently gives E not defined, will be fixed by the next SymEgine

test_fail(x,y,d) = erf(x*y/d)
println("Test non-differentiable")
NJ = @ode_def NoJacTest begin
  dx = a*x - μ*x*y
  dy = -c*y + test_fail(x,y,d)
end a μ c d
NJ(du,u,[1.5,1,3,4],t)
@test du == [-3.0;-3*3.0 + erf(2.0*3.0/4)]
@test du == NJ(u,[1.5,1,3,4],t)
@test_throws MethodError NJ.jac(iJ,u,p,t)
# NJ(Val{:jac},t,u,J) # Currently gives E not defined, will be fixed by the next SymEgine

println("Test anonymous definition")

f_t_noname = @ode_def begin # Checks for error due to symbol on 1
  dx = x
  dy = -c*y + d*x*y*t^2
end a b c d

@test DiffEqBase.__has_syms(f_t_noname)

f = @ode_def begin
  dx = a*x - b*x*y
  dy = -c*y + d*x*y
end a b c d
@test_nowarn f([0.1,0.2], [1,2], [1,2,3,4], 1)

# Test failures of derivatives should not have #undef
sir_ode = @ode_def SIRModel begin
  dS = -b*S*I
  dI =  b*S*I - g*I
  dR =  g*I
end b g


text:
# Collision Detection routines based on Minkowski techniques.

export
        ConvexPolygon,
        CPAMemory,
        CollisionCheckResult,

        LineSegment,

        to_oriented_bounding_box!,
        get_oriented_bounding_box,
        is_colliding,
        is_potentially_colliding,
        get_collision_time,
        get_first_collision,
        get_time_and_dist_of_closest_approach,
        is_collision_free,
        get_distance,
        get_edge,
        get_side

######################################

function cyclic_shift_left!(arr::AbstractVector, d::Int, n::Int=length(a))
    #=
    Perform a cyclic rotation of the elements in the array, in place
        d = amount to shift
        n = length of array (if we only want to work with first n elements)
    =#

    for i in 1 : gcd(d, n)
    # move i-th values of blocks

        temp = arr[i]
        j = i
        while true
            k = j + d
            if k > n
                k = k - n
            end
            if k == i
                break
            end
            arr[j] = arr[k]
            j = k
        end
        arr[j] = temp
    end

    arr
end

######################################

immutable LineSegment
    a::VecE2
    b::VecE2
end
@compat Base.:+(seg::LineSegment, v::VecE2) = LineSegment(seg.a + v, seg.b + v)
@compat Base.:-(seg::LineSegment, v::VecE2) = LineSegment(seg.a - v, seg.b - v)

function get_polar_angle(seg::LineSegment)
    θ = atan2(seg.b.y - seg.a.y, seg.b.x - seg.a.x)
    if θ < 0.0
        θ += 2π
    end
    θ
end

function get_distance(seg::LineSegment, p::VecE2)

    ab = seg.b - seg.a
    pb = p - seg.a

    denom = abs2(ab)
    if denom == 0.0
        return 0.0
    end

    r = dot(ab, pb)/denom

    if r ≤ 0.0
        abs(p - seg.a)
    elseif r ≥ 1.0
        abs(p - seg.b)
    else
        abs(p - (seg.a + r*ab))
    end
end

get_side(seg::LineSegment, p::VecE2) = sign((seg.b.x-seg.a.x) * (p.y-seg.a.y) - (seg.b.y-seg.a.y) * (p.x-seg.a.x))

######################################

function get_signed_area(pts::Vector{VecE2}, npts::Int = length(pts))

    # https://en.wikipedia.org/wiki/Shoelace_formula
    # sign of -1 means clockwise, sign of 1 means counterclockwise

    retval = pts[npts].x*pts[1].y - pts[1].x*pts[npts].y
    for i in 1 : npts-1
        retval += pts[i].x * pts[i+1].y
        retval -= pts[i+1].x * pts[i].y
    end

    retval / 2
end
function get_edge(pts::Vector{VecE2}, i::Int, npts::Int=length(pts))
    a = pts[i]
    b = i+1 ≤ npts ? pts[i+1] : pts[1]
    LineSegment(a,b)
end

######################################

type ConvexPolygon
    pts::Vector{VecE2} # ordered counterclockwise along polygon boundary s.t. first edge has minimum polar angle in [0,2π]
    npts::Int # number of pts that are currently used (since we preallocate a longer array)

    ConvexPolygon(npts::Int) = new(Array(VecE2, npts), 0)
end

Base.start(poly::ConvexPolygon) = 1
Base.done(poly::ConvexPolygon, i::Int) = i > length(poly)
Base.next(poly::ConvexPolygon, i::Int) = (poly.pts[i], i+1)

Base.length(poly::ConvexPolygon) = poly.npts
Base.isempty(poly::ConvexPolygon) = poly.npts == 0
get_edge(P::ConvexPolygon, i::Int) = get_edge(P.pts, i, P.npts)
get_signed_area(poly::ConvexPolygon) = get_signed_area(poly.pts, poly.npts)
function Base.empty!(poly::ConvexPolygon)
    poly.npts = 0
    poly
end
function Base.copy!(dest::ConvexPolygon, src::ConvexPolygon)
    dest.npts = src.npts
    copy!(dest.pts, 1, src.pts, 1, src.npts)
    dest
end
function Base.push!(poly::ConvexPolygon, v::VecE2)
    poly.pts[poly.npts+1] = v
    poly.npts += 1
    poly
end
function Base.contains(poly::ConvexPolygon, v::VecE2)

    # does NOT include pts on the physical boundary

    previous_side = 0

    for i in 1 : length(poly)
        seg = get_edge(poly, i)
        affine_segment = seg.b - seg.a
        affine_point = v - seg.a
        current_side = round(Int, sign(cross(affine_segment, affine_point))) # sign indicates side
        if current_side == 0
            # outside or over an edge
            return false
        elseif previous_side == 0
            # first segment
            previous_side = current_side
        elseif previous_side != current_side
            # all most be on the same side
            return false
        end
    end
    true
end
Base.contains(poly::ConvexPolygon, v::VecSE2) = contains(poly, convert(VecE2, v))
function Base.show(io::IO, poly::ConvexPolygon)
    @printf(io, "ConvexPolygon: len %d (max %d pts)\n", length(poly), length(poly.pts))
    for i in 1 : length(poly)
        print(io, "\t"); show(io, poly.pts[i])
        print(io, "\n")
    end
end

AutomotiveDrivingModels.get_center(poly::ConvexPolygon) = sum(poly.pts) / poly.npts
function get_distance(poly::ConvexPolygon, v::VecE2; solid::Bool=true)
    if solid && contains(poly, v)
        0.0
    else
        min_dist = Inf
        for i in 1 : length(poly)
            seg = get_edge(poly, i)
            min_dist = min(min_dist, get_distance(seg, v))
        end
        min_dist
    end
end

function ensure_pts_sorted_by_min_polar_angle!(poly::ConvexPolygon, npts::Int=poly.npts)

    @assert(npts ≥ 3)
    @assert(sign(get_signed_area(poly)) == 1) # must be counter-clockwise

    # ensure that edges are sorted by minimum polar angle in [0,2π]

    angle_start = Inf
    index_start = -1
    for i in 1 : npts
        seg = get_edge(poly.pts, i)

        θ = atan2(seg.b.y - seg.a.y, seg.b.x - seg.a.x)
        if θ < 0.0
            θ += 2π
        end

        if θ < angle_start
            angle_start = θ
            index_start = i
        end
    end

    if index_start != 1
        cyclic_shift_left!(poly.pts, index_start-1, npts)
    end
    poly
end
function shift!(poly::ConvexPolygon, v::VecE2)
    for i in 1 : length(poly)
        poly.pts[i] += v
    end
    ensure_pts_sorted_by_min_polar_angle!(poly)
    poly
end
function rotate!(poly::ConvexPolygon, θ::Float64)
    for i in 1 : length(poly)
        poly.pts[i] = Vec.rot(poly.pts[i], θ)
    end
    ensure_pts_sorted_by_min_polar_angle!(poly)
    poly
end
function mirror!(poly::ConvexPolygon)
    for i in 1 : length(poly)
        poly.pts[i] = -poly.pts[i]
    end
    ensure_pts_sorted_by_min_polar_angle!(poly)
    poly
end

function minkowksi_sum!(retval::ConvexPolygon, P::ConvexPolygon, Q::ConvexPolygon)

    #=
    For two convex polygons P and Q in the plane with m and n vertices, their Minkowski sum is a
    convex polygon with at most m + n vertices and may be computed in time O (m + n) by a very simple procedure,
    which may be informally described as follows.

    Assume that the edges of a polygon are given and the direction, say, counterclockwise, along the polygon boundary.
    Then it is easily seen that these edges of the convex polygon are ordered by polar angle.
    Let us merge the ordered sequences of the directed edges from P and Q into a single ordered sequence S.
    Imagine that these edges are solid arrows which can be moved freely while keeping them parallel to their original direction.
    Assemble these arrows in the order of the sequence S by attaching the tail of the next arrow to the head of the previous arrow.
    It turns out that the resulting polygonal chain will in fact be a convex polygon which is the Minkowski sum of P and Q.
    =#

    empty!(retval)

    index_P = 1
    index_Q = 1
    θp = get_polar_angle(get_edge(P, index_P))
    θq = get_polar_angle(get_edge(Q, index_Q))

    while index_P ≤ length(P) || index_Q ≤ length(Q)
        # select next edge with minimum polar angle

        if θp == θq
            seg_p = get_edge(P, index_P)
            seg_q = get_edge(Q, index_Q)
            O = isempty(retval) ? P.pts[1] + Q.pts[1] : retval.pts[retval.npts]
            push!(retval, O + seg_p.b - seg_p.a + seg_q.b - seg_q.a)
            index_P += 1
            θp = index_P ≤ length(P) ? get_polar_angle(get_edge(P, index_P)) : Inf
            index_Q += 1
            θq = index_Q ≤ length(Q) ? get_polar_angle(get_edge(Q, index_Q)) : Inf
        elseif θp ≤ θq
            seg = get_edge(P, index_P)
            O = isempty(retval) ? P.pts[1] + Q.pts[1] : retval.pts[retval.npts]
            push!(retval, O + seg.b - seg.a)
            index_P += 1
            θp = index_P ≤ length(P) ? get_polar_angle(get_edge(P, index_P)) : Inf
        else
            seg = get_edge(Q, index_Q)
            O = isempty(retval) ? P.pts[1] + Q.pts[1] : retval.pts[retval.npts]
            push!(retval, O + seg.b - seg.a)
            index_Q += 1
            θq = index_Q ≤ length(Q) ? get_polar_angle(get_edge(Q, index_Q)) : Inf
        end
    end

    ensure_pts_sorted_by_min_polar_angle!(retval)

    retval
end
function minkowski_difference!(retval::ConvexPolygon, P::ConvexPolygon, Q::ConvexPolygon)

    #=
    The minkowski difference is what you get by taking the minkowski sum of a shape and the mirror of another shape.
    So, your second shape gets flipped about the origin (all of its points are negated).

    The idea is that you do a binary operation on two shapes to get a new shape,
    and if the origin (the zero vector) is inside that shape, then they are colliding.
    =#

    minkowksi_sum!(retval, P, mirror!(Q))
end

function is_colliding(P::ConvexPolygon, Q::ConvexPolygon, temp::ConvexPolygon=ConvexPolygon(length(P) + length(Q)))
    minkowski_difference!(temp, P, Q)
    contains(temp, VecE2(0,0))
end
function get_distance(P::ConvexPolygon, Q::ConvexPolygon, temp::ConvexPolygon=ConvexPolygon(length(P) + length(Q)))
    minkowski_difference!(temp, P, Q)
    get_distance(temp, VecE2(0,0))
end

function to_oriented_bounding_box!(retval::ConvexPolygon, center::VecSE2, len::Float64, wid::Float64)

    @assert(len > 0)
    @assert(wid > 0)
    @assert(!isnan(center.θ))
    @assert(!isnan(center.x))
    @assert(!isnan(center.y))

    x = polar(len/2, center.θ)
    y = polar(wid/2, center.θ+π/2)

    C = convert(VecE2,center)
    retval.pts[1] =  x - y + C
    retval.pts[2] =  x + y + C
    retval.pts[3] = -x + y + C
    retval.pts[4] = -x - y + C
    retval.npts = 4

    AutomotiveDrivingModels.ensure_pts_sorted_by_min_polar_angle!(retval)

    retval
end
get_oriented_bounding_box(center::VecSE2, len::Float64, wid::Float64) = to_oriented_bounding_box!(ConvexPolygon(4), center, len, wid)
function to_oriented_bounding_box!(retval::ConvexPolygon, veh::Vehicle, center::VecSE2 = get_center(veh))

    # get an oriented bounding box at the vehicle's position

    to_oriented_bounding_box!(retval, center, veh.def.length, veh.def.width)

    retval
end
get_oriented_bounding_box(veh::Vehicle, center::VecE2 = convert(VecE2, get_center(veh))) = get_oriented_bounding_box(ConvexPolygon(4), veh, center)

######################################

get_future_position(ray::VecSE2, speed::Float64, t::Float64) = convert(VecE2, ray) + polar(speed*t, ray.θ)
function is_colliding(ray::VecSE2, seg::LineSegment)
    v₁ = convert(VecE2, ray) - seg.a
    v₂ = seg.b - seg.a
    v₃ = polar(1.0, ray.θ + π/2)

    denom = dot(v₂, v₃)

    if abs(denom) > 0.0
        t₁ = cross(v₂, v₁) / denom # time for ray (0 ≤ t₁)
        t₂ = dot(v₁, v₃) / denom # time for segment (0 ≤ t₂ ≤ 1)
        return 0.0 ≤ t₁ && 0.0 ≤ t₂ ≤ 1.0
    else
        # denom is zero if the segment and the ray are parallel
        # only collide if they are perfectly aligned
        return are_collinear(ray, seg.a, seg.b)
    end
end
function get_collision_time(ray::VecSE2, seg::LineSegment, ray_speed::Float64)

    o = convert(VecE2, ray)
    v₁ = o - seg.a
    v₂ = seg.b - seg.a
    v₃ = polar(1.0, ray.θ + π/2)

    denom = dot(v₂, v₃)

    if abs(denom) > 0.0
        d₁ = cross(v₂, v₁) / denom # time for ray (0 ≤ t₁)
        t₂ = dot(v₁, v₃) / denom # time for segment (0 ≤ t₂ ≤ 1)
        if 0.0 ≤ d₁ && 0.0 ≤ t₂ ≤ 1.0
            return d₁/ray_speed
        end
    else
        # denom is zero if the segment and the ray are parallel
        # only collide if they are perfectly aligned
        if are_collinear(ray, seg.a, seg.b)
            dist_a = abs2(seg.a - o)
            dist_b = abs2(seg.b - o)
            return sqrt(min(dist_a, dist_b)) / ray_speed
        end
    end

    NaN
end
function get_time_and_dist_of_closest_approach(ray::VecSE2, seg::LineSegment, ray_speed::Float64, if_no_col_skip_eval::Bool=false)

    o = convert(VecE2, ray)
    v₁ = o - seg.a
    v₂ = seg.b - seg.a
    v₃ = polar(1.0, ray.θ + π/2)

    denom = dot(v₂, v₃)

    if abs(denom) > 0.0

        d₁ = cross(v₂, v₁) / denom # time for ray (0 ≤ d₁)
        t₂ = dot(v₁, v₃) / denom # time for segment (0 ≤ t₂ ≤ 1)

        if 0.0 ≤ d₁ && 0.0 ≤ t₂ ≤ 1.0
            t = d₁/ray_speed
            return (t, 0.0)
        else
            # no collision, get time of closest approach to each endpoint
            if !if_no_col_skip_eval
                ray_pt = polar(1.0, ray.θ)
                projA = proj(seg.a - o, ray_pt, VecE2)
                projB = proj(seg.b - o, ray_pt, VecE2)
                tA = max(abs(projA)/ray_speed * sign(dot(ray_pt, projA)), 0.0)
                tB = max(abs(projB)/ray_speed * sign(dot(ray_pt, projB)), 0.0)
                pA = get_future_position(ray, ray_speed, tA)
                pB = get_future_position(ray, ray_speed, tB)
                distA = abs2(seg.a - pA)
                distB = abs2(seg.b - pB)
                if distA < distB
                    return (tA, sqrt(distA))
                else
                    return (tB, sqrt(distB))
                end
            else
                return (Inf, Inf)
            end
        end
    else
        # denom is zero if the segment and the ray are parallel
        # only collide if they are perfectly aligned
        if are_collinear(ray, seg.a, seg.b)
            dist_a = abs2(seg.a - o)
            dist_b = abs2(seg.b - o)
            t = sqrt(min(dist_a, dist_b)) / ray_speed
            return (t, 0.0)
        else
            if !if_no_col_skip_eval
                ray_pt = polar(1.0, ray.θ)
                projA = proj(seg.a - o, ray_pt, VecE2)
                projB = proj(seg.b - o, ray_pt, VecE2)
                tA = max(abs(projA)/ray_speed * sign(dot(ray_pt, projA)), 0.0)
                tB = max(abs(projB)/ray_speed * sign(dot(ray_pt, projB)), 0.0)
                pA = get_future_position(ray, ray_speed, tA)
                pB = get_future_position(ray, ray_speed, tB)
                distA = abs2(seg.a - pA)
                distB = abs2(seg.b - pB)
                if distA < distB
                    return (tA, distA)
                else
                    return (tB, distB)
                end
            else
                return (Inf,Inf)
            end
        end
    end
end
function is_colliding(ray::VecSE2, poly::ConvexPolygon)
    # collides if at least one of the segments collides with the ray

    for i in 1 : length(poly)
        seg = get_edge(poly, i)
        if is_colliding(ray, seg)
            return true
        end
    end
    false
end
function get_collision_time(ray::VecSE2, poly::ConvexPolygon, ray_speed::Float64)
    min_col_time = Inf
    for i in 1 : length(poly)
        seg = get_edge(poly, i)
        col_time = get_collision_time(ray, seg, ray_speed)
        if !isnan(col_time) && col_time < min_col_time
            min_col_time = col_time
        end
    end
    min_col_time
end

######################################

immutable CPAMemory
    vehA::ConvexPolygon # bounding box for vehicle A
    vehB::ConvexPolygon # bounding box for vehicle B
    mink::ConvexPolygon # minkowski bounding box

    CPAMemory() = new(ConvexPolygon(4), ConvexPolygon(4), ConvexPolygon(8))
end
is_colliding(mem::CPAMemory) = is_colliding(mem.vehA, mem.vehB, mem.mink)
get_distance(mem::CPAMemory) = get_distance(mem.vehA, mem.vehB, mem.mink)
function get_time_and_dist_of_closest_approach(a::Vehicle, b::Vehicle, mem::CPAMemory=CPAMemory())

    to_oriented_bounding_box!(mem.vehA, a)
    to_oriented_bounding_box!(mem.vehB, b)
    minkowksi_sum!(mem.mink, mem.vehA, mem.vehB)

    rel_pos = convert(VecE2, b.state.posG) - a.state.posG
    rel_velocity = polar(b.state.v, b.state.posG.θ) - polar(a.state.v, a.state.posG.θ)
    ray_speed = hypot(rel_velocity)
    ray = VecSE2(rel_pos, atan2(rel_velocity))

    if contains(mem.mink, convert(VecE2, ray))
        return (0.0, 0.0)
    end

    best_t_CPA = NaN
    best_d_CPA = Inf
    if_no_col_skip_eval = false

    for i in 1 : length(mem.mink)
        seg = get_edge(mem.mink, i)
        t_CPA, d_CPA = get_time_and_dist_of_closest_approach(ray, seg, ray_speed, if_no_col_skip_eval)
        if d_CPA < best_d_CPA
            best_t_CPA = t_CPA
            best_d_CPA = d_CPA
            if d_CPA == 0.0
                if_no_col_skip_eval = true
            end
        end
    end

    (best_t_CPA, best_d_CPA)
end

_bounding_radius(veh::Vehicle) = sqrt(veh.def.length*veh.def.length/4 + veh.def.width*veh.def.width/4)

"""
A fast collision check to remove things clearly not colliding
"""
function is_potentially_colliding(A::Vehicle, B::Vehicle)
    Δ² = abs2(A.state.posG - B.state.posG)
    r_a = _bounding_radius(A)
    r_b = _bounding_radius(B)
    Δ² ≤ r_a*r_a + 2*r_a*r_b + r_b*r_b
end

function is_colliding(A::Vehicle, B::Vehicle, mem::CPAMemory=CPAMemory())
    if is_potentially_colliding(A, B)
        to_oriented_bounding_box!(mem.vehA, A)
        to_oriented_bounding_box!(mem.vehB, B)
        return is_colliding(mem)
    end
    false
end
function get_distance(A::Vehicle, B::Vehicle, mem::CPAMemory=CPAMemory())
    to_oriented_bounding_box!(mem.vehA, A)
    to_oriented_bounding_box!(mem.vehB, B)
    get_distance(mem)
end

immutable CollisionCheckResult
    is_colliding::Bool
    A::Int # index of 1st vehicle
    B::Int # index of 2nd vehicle
    # TODO: penetration vector?
end

"""
Loops through the scene and finds the first collision between a vehicle and scene[target_index]
"""
function get_first_collision(scene::Scene, target_index::Int, mem::CPAMemory=CPAMemory())
    A = target_index
    vehA = scene[A]
    to_oriented_bounding_box!(mem.vehA, vehA)
    for (B,vehB) in enumerate(scene)
        if B != A
            to_oriented_bounding_box!(mem.vehB, vehB)
            if is_potentially_colliding(vehA, vehB) && is_colliding(mem)
                return CollisionCheckResult(true, A, B)
            end
        end
    end

    CollisionCheckResult(false, A, 0)
end

"""
Loops through the scene and finds the first collision between any two vehicles
"""
function get_first_collision(scene::Scene, vehicle_indeces::AbstractVector{Int}, mem::CPAMemory=CPAMemory())

    N = length(vehicle_indeces)
    for (a,A) in enumerate(vehicle_indeces)
        vehA = scene[A]
        to_oriented_bounding_box!(mem.vehA, vehA)
        for b in a +1 : length(vehicle_indeces)
            B = vehicle_indeces[b]
            vehB = scene[B]
            if is_potentially_colliding(vehA, vehB)
                to_oriented_bounding_box!(mem.vehB, vehB)
                if is_colliding(mem)
                    return CollisionCheckResult(true, A, B)
                end
            end
        end
    end

    CollisionCheckResult(false, 0, 0)
end
get_first_collision(scene::Scene, mem::CPAMemory=CPAMemory()) = get_first_collision(scene, 1:length(scene), mem)
is_collision_free(scene::Scene, mem::CPAMemory=CPAMemory()) = !(get_first_collision(scene, mem).is_colliding)
is_collision_free(scene::Scene, vehicle_indeces::AbstractVector{Int}, mem::CPAMemory=CPAMemory()) = get_first_collision(scene, vehicle_indeces, mem).is_colliding


text:
include("sawtooth.jl")
include("clvs.jl")
include("lss.jl")
using Test
function test_dsawtooth()
	s = rand()
	m = 10
	u_trj = rand(m)
	eps = 1.e-4
	du_fd = (sawtooth(u_trj .+ eps,
						s, 1)[2,:] -
				sawtooth(u_trj .- eps,
						s, 1)[2,:])./(2*eps)
	du = dsawtooth(u_trj, s)
	@assert all(isapprox.(du_fd, du)) == 1
end
function test_les()
	s = 0.01
	m = 1
	n = 25000
	u0 = rand(m)
	u_trj = sawtooth(u0, s, n)[:,1]
	n = n+1
	du_trj = reshape(dsawtooth(u_trj, s),1,1,n)
	les, clv_trj = clvs(du_trj,1)
	println(les)
	@assert all(isapprox.(les, log(2),
						  rtol=1.e-1)) == 1
end
function test_lss()
	s = 0.01
    m = 1
    n = 5000
    u0 = rand(m)
	d = 1
    d_u = 1
	dJ = ones(d,n+1)
	n_samples = 10
	dJds = zeros(n_samples)
	vsh = zeros(d,n+1)

	for i=1:n_samples
		println("Starting LSS, sample ", i)
		u_trj = sawtooth(u0, s, n)[:,1]
		du_trj = reshape(dsawtooth(u_trj, s),d,d,n+1)
		X = reshape(perturbation(u_trj,s),d,n+1)
		f = zeros(1,n+1)
		J = u_trj
		u_trj = reshape(u_trj, n+1, d)
		y, dJds[i] = lss(u_trj,
						du_trj, X, f, J, dJ,
						  s, d_u)
		println(dJds[i])
		vsh .= y
		u0 .= u_trj[end,:]
	end
	println("Sensitivities are : ", dJds)
	@test isapprox((sum(dJds)/n_samples),-0.61,rtol=0.1)
	return vsh, dJds
end

max_stars_repo_name: nishaChandramoorthy/energies

max_issues_repo_name: nishaChandramoorthy/energies

max_forks_repo_name: nishaChandramoorthy/energies


text:
module ConsoleInput

include("./utils.jl")

DlmType = Union{
    AbstractString,
    AbstractChar,
    Regex
}

function readInt(io::IO=stdin, delimiter::DlmType=" ")
    readline(io) |>
    x -> split(x, delimiter) |>
    x -> parse.(Int64, x) |>
    process
end

function readString(io::IO=stdin, delimiter::DlmType=" ")
    readline(io) |>
    x -> split(x, delimiter) |>
    process
end

function readGeneral(type, io::IO=stdin, delimiter::DlmType=" ")
    readline(io) |>
    x -> split(x, delimiter) |>
    x -> parse.(type, x) |>
    process
end

end # module


text:
module FrostHelperEntityRainbowifyController

using ..Ahorn, Maple

@mapdef Entity "FrostHelper/EntityRainbowifyController" EntityRainbowifyController(x::Integer, y::Integer, types::String="")

const placements = Ahorn.PlacementDict(
	"Entity Rainbowifier (Frost Helper)" => Ahorn.EntityPlacement(
		EntityRainbowifyController
	)
)

const sprite = "editor/FrostHelper/EntityRainbowifyController"

function Ahorn.selection(entity::EntityRainbowifyController)
	x, y = Ahorn.position(entity)

	return Ahorn.getSpriteRectangle(sprite, x, y)
end

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::EntityRainbowifyController, room::Maple.Room)
	x, y = Ahorn.position(entity)

	Ahorn.drawSprite(ctx, sprite, x, y)
end

end

max_stars_repo_name: JaThePlayer/FrostHelper

max_issues_repo_name: JaThePlayer/FrostHelper

max_forks_repo_name: JaThePlayer/FrostHelper


text:
function QRDecomposition(arg0::RealMatrix)
    return QRDecomposition((RealMatrix,), arg0)
end

function QRDecomposition(arg0::RealMatrix, arg1::jdouble)
    return QRDecomposition((RealMatrix, jdouble), arg0, arg1)
end

function get_h(obj::QRDecomposition)
    return jcall(obj, "getH", RealMatrix, ())
end

function get_q(obj::QRDecomposition)
    return jcall(obj, "getQ", RealMatrix, ())
end

function get_qt(obj::QRDecomposition)
    return jcall(obj, "getQT", RealMatrix, ())
end

function get_r(obj::QRDecomposition)
    return jcall(obj, "getR", RealMatrix, ())
end

function get_solver(obj::QRDecomposition)
    return jcall(obj, "getSolver", DecompositionSolver, ())
end


text:
"""
$(DocStringExtensions.README)

# Imports
$(DocStringExtensions.IMPORTS)

# Exports
$(DocStringExtensions.EXPORTS)
"""
module PNML

# Width for printing.
if !haskey(ENV, "COLUMNS")
    ENV["COLUMNS"] = 180
end

using EzXML
using MLStyle: @match
using DocStringExtensions
using AutoHashEquals
using PrettyPrinting
using AbstractTrees
using LabelledArrays
using Reexport

include("config.jl")

#include("Base/PnmlBase.jl") #TODO sub module test
#@reexport using .PnmlBase

include("xmlutils.jl")
include("id.jl")
include("interfaces.jl")
include("types.jl")
include("exceptions.jl")
include("utils.jl")

include("pnmltypes.jl")
@reexport using .PnmlTypes

include("IR/intermediate.jl")

include("Parse/parseutils.jl")
include("Parse/anyelement.jl")
include("Parse/parse.jl")
include("Parse/graphics.jl")
include("Parse/declarations.jl")
include("Parse/toolspecific.jl")
include("Parse/maps.jl")

include("Net/petrinet.jl")
include("Net/simplenet.jl")
include("Net/hlnet.jl")

export @xml_str, PnmlDict,
    parse_str, parse_file, parse_pnml, parse_node,
    PnmlException, MissingIDException, MalformedException
end