text:
module Shapefile

import GeoInterface, DBFTables, Tables

struct Rect
    left::Float64
    bottom::Float64
    right::Float64
    top::Float64
end

struct Interval
    left::Float64
    right::Float64
end

struct Point <: GeoInterface.AbstractPoint
    x::Float64
    y::Float64
end

struct PointM <: GeoInterface.AbstractPoint
    x::Float64
    y::Float64
    m::Float64  # measure
end

struct PointZ <: GeoInterface.AbstractPoint
    x::Float64
    y::Float64
    z::Float64
    m::Float64  # measure
end

struct Polyline <: GeoInterface.AbstractMultiLineString
    MBR::Rect
    parts::Vector{Int32}
    points::Vector{Point}
end

struct PolylineM <: GeoInterface.AbstractMultiLineString
    MBR::Rect
    parts::Vector{Int32}
    points::Vector{Point}
    measures::Vector{Float64}
end

struct PolylineZ <: GeoInterface.AbstractMultiLineString
    MBR::Rect
    parts::Vector{Int32}
    points::Vector{Point}
    zvalues::Vector{Float64}
    measures::Vector{Float64}
end

struct Polygon <: GeoInterface.AbstractMultiPolygon
    MBR::Rect
    parts::Vector{Int32}
    points::Vector{Point}
end

Base.show(io::IO, p::Polygon) =
    print(io, "Polygon(", length(p.points), " Points)")

struct PolygonM <: GeoInterface.AbstractMultiPolygon
    MBR::Rect
    parts::Vector{Int32}
    points::Vector{Point}
    measures::Vector{Float64}
end

struct PolygonZ <: GeoInterface.AbstractMultiPolygon
    MBR::Rect
    parts::Vector{Int32}
    points::Vector{Point}
    zvalues::Vector{Float64}
    measures::Vector{Float64}
end

struct MultiPoint <: GeoInterface.AbstractMultiPoint
    MBR::Rect
    points::Vector{Point}
end

struct MultiPointM <: GeoInterface.AbstractMultiPoint
    MBR::Rect
    points::Vector{Point}
    measures::Vector{Float64}
end

struct MultiPointZ <: GeoInterface.AbstractMultiPoint
    MBR::Rect
    points::Vector{Point}
    zvalues::Vector{Float64}
    measures::Vector{Float64}
end

struct MultiPatch <: GeoInterface.AbstractGeometry
    MBR::Rect
    parts::Vector{Int32}
    parttypes::Vector{Int32}
    points::Vector{Point}
    zvalues::Vector{Float64}
    # measures::Vector{Float64}  # (optional)
end

const SHAPETYPE = Dict{Int32,DataType}(
    0 => Missing,
    1 => Point,
    3 => Polyline,
    5 => Polygon,
    8 => MultiPoint,
    11 => PointZ,
    13 => PolylineZ,
    15 => PolygonZ,
    18 => MultiPointZ,
    21 => PointM,
    23 => PolylineM,
    25 => PolygonM,
    28 => MultiPointM,
    31 => MultiPatch,
)

mutable struct Handle{T<:Union{<:GeoInterface.AbstractGeometry,Missing}}
    code::Int32
    length::Int32
    version::Int32
    shapeType::Int32
    MBR::Rect
    zrange::Interval
    mrange::Interval
    shapes::Vector{T}
end

Base.length(shp::Handle) = length(shp.shapes)

function Base.read(io::IO, ::Type{Rect})
    minx = read(io, Float64)
    miny = read(io, Float64)
    maxx = read(io, Float64)
    maxy = read(io, Float64)
    Rect(minx, miny, maxx, maxy)
end

function Base.read(io::IO, ::Type{Point})
    x = read(io, Float64)
    y = read(io, Float64)
    Point(x, y)
end

function Base.read(io::IO, ::Type{PointM})
    x = read(io, Float64)
    y = read(io, Float64)
    m = read(io, Float64)
    PointM(x, y, m)
end

function Base.read(io::IO, ::Type{PointZ})
    x = read(io, Float64)
    y = read(io, Float64)
    z = read(io, Float64)
    m = read(io, Float64)
    PointZ(x, y, z, m)
end

function Base.read(io::IO, ::Type{Polyline})
    box = read(io, Rect)
    numparts = read(io, Int32)
    numpoints = read(io, Int32)
    parts = Vector{Int32}(undef, numparts)
    read!(io, parts)
    points = Vector{Point}(undef, numpoints)
    read!(io, points)
    Polyline(box, parts, points)
end

function Base.read(io::IO, ::Type{PolylineM})
    box = read(io, Rect)
    numparts = read(io, Int32)
    numpoints = read(io, Int32)
    parts = Vector{Int32}(undef, numparts)
    read!(io, parts)
    points = Vector{Point}(undef, numpoints)
    read!(io, points)
    mrange = Vector{Float64}(undef, 2)
    read!(io, mrange)
    measures = Vector{Float64}(undef, numpoints)
    read!(io, measures)
    PolylineM(box, parts, points, measures)
end

function Base.read(io::IO, ::Type{PolylineZ})
    box = read(io, Rect)
    numparts = read(io, Int32)
    numpoints = read(io, Int32)
    parts = Vector{Int32}(undef, numparts)
    read!(io, parts)
    points = Vector{Point}(undef, numpoints)
    read!(io, points)
    zrange = Vector{Float64}(undef, 2)
    read!(io, zrange)
    zvalues = Vector{Float64}(undef, numpoints)
    read!(io, zvalues)
    mrange = Vector{Float64}(undef, 2)
    read!(io, mrange)
    measures = Vector{Float64}(undef, numpoints)
    read!(io, measures)
    PolylineZ(box, parts, points, zvalues, measures)
end

function Base.read(io::IO, ::Type{Polygon})
    box = read(io, Rect)
    numparts = read(io, Int32)
    numpoints = read(io, Int32)
    parts = Vector{Int32}(undef, numparts)
    read!(io, parts)
    points = Vector{Point}(undef, numpoints)
    read!(io, points)
    Polygon(box, parts, points)
end

function Base.read(io::IO, ::Type{PolygonM})
    box = read(io, Rect)
    numparts = read(io, Int32)
    numpoints = read(io, Int32)
    parts = Vector{Int32}(undef, numparts)
    read!(io, parts)
    points = Vector{Point}(undef, numpoints)
    read!(io, points)
    mrange = Vector{Float64}(undef, 2)
    read!(io, mrange)
    measures = Vector{Float64}(undef, numpoints)
    read!(io, measures)
    PolygonM(box, parts, points, measures)
end

function Base.read(io::IO, ::Type{PolygonZ})
    box = read(io, Rect)
    numparts = read(io, Int32)
    numpoints = read(io, Int32)
    parts = Vector{Int32}(undef, numparts)
    read!(io, parts)
    points = Vector{Point}(undef, numpoints)
    read!(io, points)
    zrange = Vector{Float64}(undef, 2)
    read!(io, zrange)
    zvalues = Vector{Float64}(undef, numpoints)
    read!(io, zvalues)
    mrange = Vector{Float64}(undef, 2)
    read!(io, mrange)
    measures = Vector{Float64}(undef, numpoints)
    read!(io, measures)
    PolygonZ(box, parts, points, zvalues, measures)
end

function Base.read(io::IO, ::Type{MultiPoint})
    box = read(io, Rect)
    numpoints = read(io, Int32)
    points = Vector{Point}(undef, numpoints)
    read!(io, points)
    MultiPoint(box, points)
end

function Base.read(io::IO, ::Type{MultiPointM})
    box = read(io, Rect)
    numpoints = read(io, Int32)
    points = Vector{Point}(undef, numpoints)
    read!(io, points)
    mrange = Vector{Float64}(undef, 2)
    read!(io, mrange)
    measures = Vector{Float64}(undef, numpoints)
    read!(io, measures)
    MultiPointM(box, points, measures)
end

function Base.read(io::IO, ::Type{MultiPointZ})
    box = read(io, Rect)
    numpoints = read(io, Int32)
    points = Vector{Point}(undef, numpoints)
    read!(io, points)
    zrange = Vector{Float64}(undef, 2)
    read!(io, zrange)
    zvalues = Vector{Float64}(undef, numpoints)
    read!(io, zvalues)
    mrange = Vector{Float64}(undef, 2)
    read!(io, mrange)
    measures = Vector{Float64}(undef, numpoints)
    read!(io, measures)
    MultiPointZ(box, points, zvalues, measures)
end

function Base.read(io::IO, ::Type{MultiPatch})
    box = read(io, Rect)
    numparts = read(io, Int32)
    numpoints = read(io, Int32)
    parts = Vector{Int32}(undef, numparts)
    read!(io, parts)
    parttypes = Vector{Int32}(undef, numparts)
    read!(io, parttypes)
    points = Vector{Point}(undef, numpoints)
    read!(io, points)
    zrange = Vector{Float64}(undef, 2)
    read!(io, zrange)
    zvalues = Vector{Float64}(undef, numpoints)
    read!(io, zvalues)
    # mrange = Vector{Float64}(2)
    # read!(io, mrange)
    # measures = Vector{Float64}(numpoints)
    # read!(io, measures)
    MultiPatch(box, parts, parttypes, points, zvalues) #,measures)
end

function Base.read(io::IO, ::Type{Handle})
    code = bswap(read(io, Int32))
    read!(io, Vector{Int32}(undef, 5))
    fileSize = bswap(read(io, Int32))
    version = read(io, Int32)
    shapeType = read(io, Int32)
    MBR = read(io, Rect)
    zmin = read(io, Float64)
    zmax = read(io, Float64)
    mmin = read(io, Float64)
    mmax = read(io, Float64)
    jltype = SHAPETYPE[shapeType]
    shapes = Vector{Union{jltype,Missing}}(undef, 0)
    file = Handle(
        code,
        fileSize,
        version,
        shapeType,
        MBR,
        Interval(zmin, zmax),
        Interval(mmin, mmax),
        shapes,
    )
    while (!eof(io))
        num = bswap(read(io, Int32))
        rlength = bswap(read(io, Int32))
        shapeType = read(io, Int32)
        if shapeType === Int32(0)
            push!(shapes, missing)
        else
            push!(shapes, read(io, jltype))
        end
    end
    file
end

function Base.:(==)(a::Rect, b::Rect)
    a.left == b.left &&
    a.bottom == b.bottom && a.right == b.right && a.top == b.top
end

include("table.jl")
include("geo_interface.jl")
include("shx.jl")

end # module


text:
# This file is part of EndemicModel project, which is licensed under BDS-3-Clause.
# See file LICENSE.md for more information.

const _OPTIONS = Dict{Symbol,Any}()

"""
    option(key::Symbol)

Retrieve default option value associated with the given `key`. Returns `default` if no
option is found. See the documentation of `option(::Val{key})` for a given `key`.
"""
function option(key::Symbol, default = missing)
    get(_OPTIONS, key, default)
end

function option(::Val, default = missing)
    default
end

"""
    option!(key::Symbol, value)

Set default option `value` associated with the given `key`. Note that `value` must be
of the type defined with `@defoption`.
"""
function option!(key::Symbol, value)
    option!(Val(key), value)
end

macro defoption(key::Symbol, value, docstring::String, valuetype)
    docstring = """
    option(:$(key))::$(valuetype)

$(docstring)
"""
    key = QuoteNode(key)
    quote
        _OPTIONS[$(key)] = $(value)
        @doc($(docstring), function option(::Val{$(key)})::$(esc(valuetype))
            _OPTIONS[$(key)]
        end)
        function option!(::Val{$(key)}, value::$(esc(valuetype)))
            _OPTIONS[$(key)] = value
        end
    end
end

@defoption(
    cache_directory,
    joinpath(home(), p".cache/EndemicModel/"),
    "Directory where cache files are saved.",
    AbstractPath
)

@defoption(
    parameter_directory,
    joinpath(parent(@__PATH__), p"input/"),
    "Directory to look for input parameters.",
    PathDesignator
)

@defoption(
    cache_filename,
    "Sources",
    "Default name of cached data source file.",
    PathDesignator
)

@defoption(cache_data_type, :csv, "Default format of cached data source.", Symbol)

@defoption(
    sources_to_update,
    [:csse, :brasil_io, :world_population, :world_gdp_per_capita, :total_tests],
    "Sources to import from on update.",
    AbstractVector{Symbol}
)

@defoption(cache_nsamples, 7, "Number of cache samples to keep.", Integer)

@defoption(
    cache_table_name,
    "DATABASE",
    "Default name of the table that contains cached imported data.",
    AbstractString
)

@defoption(force_update, false, "If true, force updating data when importing.", Bool)

@defoption(update_period, Day(1), "Frequency of update.", Period)

@defoption(
    database_directory,
    joinpath(parent(@__PATH__), p"database/"),
    "Directory where databases are exported.",
    PathDesignator
)

@defoption(
    database_filename,
    "Database",
    "Default name of exported database file.",
    PathDesignator
)

@defoption(database_data_type, :csv, "Default format of database files.", Symbol)

@defoption(
    model_directory,
    joinpath(parent(@__PATH__), p"model/"),
    "Directory where models are exported.",
    PathDesignator
)

@defoption(
    model_filename,
    "Database",
    "Default name of exported model file.",
    PathDesignator
)

@defoption(model_data_type, :csv, "Default format of model files.", Symbol)

@defoption(
    parameters_subdirectory,
    p"PARAMETERS",
    "Subdirectory where model parameters are exported.",
    PathDesignator
)

@defoption(
    parameters_filename,
    "Parameters",
    "Default name of exported model parameters file.",
    PathDesignator
)

@defoption(parameters_data_type, :csv, "Default format of model parameters files.", Symbol)

@defoption(
    minimum_infected,
    1.0,
    "Minimum number of confirmed people to be considered in estimates.",
    Real
)

@defoption(
    plot_columns,
    [:confirmed, :infected, :exposed, :active, :recovered, :deaths],
    "Columns of data to be ploted by `model_plot`.",
    Vector{Symbol},
)

@defoption(
    plot_date_format,
    dateformat"dd/mm/yyyy",
    "Format of dates in plots.",
    DateFormat,
)

ext: jl


text:
println("LATHE.STATS TESTS")
println("...........")
using Lathe.stats: NormalDist, T_Dist, mean, std
# Distributions
@testset "Distributions" begin
x = [5, 10, 15]
dist = NormalDist(x)
zero_std = dist.apply(5)
@test zero_std == dist.apply(5)
end


text:
using Test, FinancialMonteCarlo, DualNumbers;
@show "HestonModel"
S0 = 100.0;
K = 100.0;
r = 0.02;
T = 1.0;
d = 0.01;
D = 90.0;

Nsim = 10000;
Nstep = 30;
sigma = 0.2;
sigma_zero = dual(0.2, 1.0);
kappa = 0.01;
theta = 0.03;
lambda = 0.01;
rho_ = 0.0;
mc = MonteCarloConfiguration(Nsim, Nstep);
mc1 = MonteCarloConfiguration(Nsim, Nstep, FinancialMonteCarlo.AntitheticMC());
toll = 0.8;

rfCurve = ZeroRate(r);

FwdData = Forward(T)
EUData = EuropeanOption(T, K)
AMData = AmericanOption(T, K)
BarrierData = BarrierOptionDownOut(T, K, D)
AsianData1 = AsianFloatingStrikeOption(T)
AsianData2 = AsianFixedStrikeOption(T, K)
Model = HestonProcess(sigma, sigma_zero, lambda, kappa, rho_, theta, Underlying(S0, d));

@show FwdPrice = pricer(Model, rfCurve, mc, FwdData);
@show EuPrice = pricer(Model, rfCurve, mc, EUData);
@show AmPrice = pricer(Model, rfCurve, mc, AMData);
@show BarrierPrice = pricer(Model, rfCurve, mc, BarrierData);
@show AsianPrice1 = pricer(Model, rfCurve, mc, AsianData1);
@show AsianPrice2 = pricer(Model, rfCurve, mc, AsianData2);

@test abs(FwdPrice - 98.8790717426047) < toll
@test abs(EuPrice - 7.9813557592924) < toll
@test abs(BarrierPrice - 7.006564636309922) < toll
@test abs(AsianPrice1 - 4.242256952013707) < toll

@show FwdPrice = pricer(Model, rfCurve, mc1, FwdData);
@show EuPrice = pricer(Model, rfCurve, mc1, EUData);
@show AmPrice = pricer(Model, rfCurve, mc1, AMData);
@show BarrierPrice = pricer(Model, rfCurve, mc1, BarrierData);
@show AsianPrice1 = pricer(Model, rfCurve, mc1, AsianData1);
@show AsianPrice2 = pricer(Model, rfCurve, mc1, AsianData2);
tollanti = 0.8
@test abs(FwdPrice - 98.8790717426047) < tollanti
@test abs(EuPrice - 7.9813557592924) < tollanti
@test abs(BarrierPrice - 7.006564636309922) < tollanti
@test abs(AsianPrice1 - 4.242256952013707) < tollanti


text:
#=##############################################################################################
Filename: solve_rotor_nondimensional.jl
Author: Ryan Anderson
Contact: rymanderson@gmail.com
README: define an `Action` object to solve a CCBlade rotor
=###############################################################################################

"""
    solve_rotor_nondimensional(aircraft, parameters, freestream, environment, step_range, stepi, step_symbol) <: Action

# Arguments:

* `aircraft::Aircraft`: `Aircraft` system object
* `parameters<:Parameters`: `Parameters` struct
* `freestream::Freestream`: `Freestream` object
* `environment::Environment` `Environment` object
* `step_range::AbstractArray`: array of times for which the simulation is run
* `stepi::Int`: index of the current step
* `step_symbol::String`: defines the step, e.g. `alpha` or `time`

`parameters <: Parameters` requires the following elements:

* `omegas::Array{Float64,2}`: [i,j]th element is the rotational velocities for the ith rotor at the jth step
* `Js::Vector{Vector{Float64}}`: a vector of length `length(step_range)` containing a vector of advance ratios for each rotor
* `Ts::Array{Float64,2}`: each [i,j]th element is the thrust of the ith rotor at the jth step
* `Qs::Array{Float64,2}`: each [i,j]th element is the torque of the ith rotor at the jth step
* `CTs::Array{Float64,2}`: each [i,j]th element is the thrust coefficient of the ith rotor at the jth step
* `CQs::Array{Float64,2}`: each [i,j]th element is the torque coefficient of the ith rotor at the jth step
* `ηs::Array{Float64,2}`: each [i,j]th element is the propulsive efficiency of the ith rotor at the jth step

"""
function solve_rotor_nondimensional(aircraft, parameters, freestream, environment, step_range, stepi, step_symbol)

    Js = view(parameters.Js,:,stepi)
    Ts = view(parameters.Ts,:,stepi)
    Qs = view(parameters.Qs,:,stepi)
    Ps = view(parameters.Ps,:,stepi)
    CTs = view(parameters.CTs,:,stepi)
    CQs = view(parameters.CQs,:,stepi)
    ηs = view(parameters.ηs,:,stepi)
    us = parameters.us[stepi]
    vs = parameters.vs[stepi]
    omegas = view(parameters.omegas,:,stepi)

    solverotorsnondimensional!(Js, Ts, Qs, Ps, CTs, CQs, ηs, us, vs, aircraft.rotor_system, omegas, freestream, environment)

    return false
end

"""
    solve_rotor_nondimensional(aircraft, step_range) <: Action

Method returns initialized elements required for the `parameters <: Parameters` struct during simulation.

# Arguments:

* `aircraft::Aircraft` : system to be simulated
* `step_range::AbstractArray` : defines each step of the simulation

# Returns:

* `omegas::Vector{Float64}` : a vector of rotational speeds in rad/s at the current step
* `Js::Array{Float64,2}` : each [i,j]th element is the advance ratio of the ith rotor at the jth step
* `Ts::Array{Float64,2}` : each [i,j]th element is the thrust of the ith rotor at the jth step
* `Qs::Array{Float64,2}` : each [i,j]th element is the torque of the ith rotor at the jth step
* `Ps::Array{Float64,2}`: [i,j]th element is the power of the ith rotor at the jth step
* `CTs::Array{Float64,2}` : each [i,j]th element is the thrust coefficient of the ith rotor at the jth step
* `CQs::Array{Float64,2}` : each [i,j]th element is the torque coefficient of the ith rotor at the jth step
* `ηs::Array{Float64,2}` : each [i,j]th element is the propulsive efficiency of the ith rotor at the jth step
* `us::Vector{Vector{Vector{Float64}}}` : each [i][j][k]th element is the axial induced velocity at ith step of the jth rotor at the kth radial section
* `vs::Vector{Vector{Vector{Float64}}}` : each [i][j][k]th element is the swirl induced velocity at ith step of the jth rotor at the kth radial section

"""
function solve_rotor_nondimensional(aircraft, step_range)

    omegas, Js, Ts, Qs, Ps, us, vs = solve_rotor(aircraft, step_range)
    CTs = deepcopy(Qs)
    CQs = deepcopy(Qs)
    ηs = deepcopy(Qs)

    return omegas, Js, Ts, Qs, Ps, CTs, CQs, ηs, us, vs
end

max_stars_repo_name: byuflowlab/AircraftSystems

max_issues_repo_name: byuflowlab/AircraftSystems

max_forks_repo_name: byuflowlab/AircraftSystems


text:
abstract type Argument end

struct PositionalArgument <: Argument
    name::AbstractString
    type
end

struct KeywordArgument <: Argument
    name::AbstractString
    default
end

function KeywordArgument(decl::AbstractString)
    name_split = split(decl, "=")
    name = first(name_split)
    default = length(name_split) == 2 ? last(name_split) : nothing
    KeywordArgument(name, default)
end
function PositionalArgument(decl::AbstractString)
    name_split = split(decl, "::")
    name = first(name_split)
    type = length(name_split) == 2 ? last(name_split) : nothing
    PositionalArgument(name, type)
end

Base.convert(T::Type{KeywordArgument}, arg::PositionalArgument) = T(arg.name, nothing)

Base.names(arg::AbstractArray{<: Argument}) = getproperty.(arg, :name)
types(arg::AbstractArray{PositionalArgument}) = getproperty.(arg, :type)

struct Signature
    name::AbstractString
    args::AbstractArray{PositionalArgument}
    kwargs::AbstractArray{KeywordArgument}
end

argnames(sig::Signature) = [getproperty.(sig.args, :name)..., getproperty.(sig.kwargs, :name)...]
argtypes(sig::Signature) = getproperty.(sig.args, :type)

function Signature(m::Method)
    args_and_types = Base.arg_decl_parts(m)[2][2:end]
    arg_names = getindex.(args_and_types, 1)
    arg_types = map(x -> isempty(x) ? nothing : x, getindex.(args_and_types, 2))
    args = PositionalArgument.(arg_names, arg_types)
    kwargs = KeywordArgument.(String.(Base.kwarg_decl(m)))
    sig = Signature(String(m.name), args, kwargs)
end

"""Get argument names from the first method of function f
"""
function Signature(f)
    ms = methods(f)
    length(ms) > 1 && @warn "Function $f has $(length(ms)) methods"
    Signature(first(ms))
end

typed_field(name, type) = isnothing(type) || type == "" ? name : join([name, type], "::")
function decompose_field_decl(typed_field)
    parts = split(typed_field, "::")
    length(parts) == 1 ? parts[1] => nothing : parts[1] => parts[2]
end

Base.string(arg::KeywordArgument) = arg.name * (isnothing(arg.default) ? "" : "=" * string(arg.default))
Base.string(arg::PositionalArgument) = typed_field(arg.name, arg.type)
Base.join(args::AbstractArray{<: Argument}) = join_args(string.(args))
Base.show(io::IO, sig::Signature) = print(io, sig.name, "(", join(sig.args), isempty(sig.kwargs) ? "" : ("; " * join(sig.kwargs)), ")")

generate(sig::Signature) = join([sig.name, "(", join(sig.args), isempty(sig.kwargs) ? "" : ("; " * join(sig.kwargs)), ")"])


text:
using Test
using LinearAlgebra
using Printf
using Random

Base.show(io::IO, f::Float64) = @printf(io, "%20.13e", f)

import Juqbox

include("evalGrad.jl")
include("test-stormer-verlet.jl")

refDir = "reference_solutions"

# test Stormer-Verlet time stepping
refFile = refDir * "/" * "err-mat-ref.jld2"
success = timestep_convergence(refFile)
println("time stepping test: pass=", success)
@test success

# test traceobjgrad for various setups
caseNames =["rabi", "swap02", "flux", "cnot2", "cnot3"]
caseDir = "cases"
nCases = length(caseNames)
pass = BitArray(undef, nCases)
pass .= false

q = 0
for case in caseNames
    global q += 1

     # setup the testcase (assign pcof0, params and wa
    juliaFile = caseDir * "/" * case * "-setup.jl"
    include(juliaFile);

    # evaluate fidelity
    caseRefFile = refDir * "/" * case * "-ref.jld2"
    pass[q] = evalObjGrad(pcof0, params, wa, caseRefFile)

    println("case=", case, " pass=", pass[q])
    @test pass[q]

end


text:
@testset "tuples" begin
    for n = 1:4
        for t in unique(combinations(repeat(1:n, n), n))
            @test all(Bennu.tuplesort(tuple(t...)) .== sort(t))
            @test all(t[collect(Bennu.tuplesortperm(tuple(t...)))] .== sort(t))
        end

        for t in permutations(1:n)
            @test all(Bennu.tuplesortperm(tuple(t...)) .== sortperm(t))
        end

        @test Bennu.permutationtuples(n) == Tuple(Tuple.(collect(permutations(1:n))))
    end

    @test (1, 2, false) == Bennu.minmaxflip(1, 2)
    @test (1, 1, false) == Bennu.minmaxflip(1, 1)
    @test (1, 2, true) == Bennu.minmaxflip(2, 1)
end


text:
module JUDI4Cloud

import Base.vcat, Base.+

using AzureClusterlessHPC, Reexport, PyCall
@reexport using JUDI

export init_culsterless

_judi_defaults = Dict("_POOL_ID"                => "JudiPool",
                    "_POOL_VM_SIZE"           => "Standard_E4s_v3",
                    "_VERBOSE"                => "0",
                    "_NODE_OS_OFFER"          => "ubuntu-server-container",
                    "_NODE_OS_PUBLISHER"      => "microsoft-azure-batch",
                    "_CONTAINER"              => ENV["CONTAINER_IMAGE"],
                    "_NODE_COUNT_PER_POOL"    => ENV["NODE_COUNT_PER_POOL"],
                    "_NUM_RETRYS"             => "1",
                    "_POOL_COUNT"             => ENV["POOL_COUNT"],
                    "_PYTHONPATH"             => "/usr/local/lib/python3.8/dist-packages",
                    "_JULIA_DEPOT_PATH"       => "/root/.julia",
                    "_OMP_NUM_THREADS"        => "4",
                    "_NODE_OS_SKU"            => "20-04-lts")

"""
Define auto scale formula to avoid idle pools
"""
auto_scale_formula(x) = """// Get pending tasks for the past 15 minutes.
\$samples = \$ActiveTasks.GetSamplePercent(TimeInterval_Minute * 15);
// If we have fewer than 70 percent data points, we use the last sample point, otherwise we use the maximum of last sample point and the history average.
\$tasks = \$samples < 70 ? max(0, \$ActiveTasks.GetSample(1)) :
max( \$ActiveTasks.GetSample(1), avg(\$ActiveTasks.GetSample(TimeInterval_Minute * 15)));
// If number of pending tasks is not 0, set targetVM to pending tasks, otherwise 25% of current dedicated.
\$targetVMs = \$tasks > 0 ? \$tasks : max(0, \$TargetDedicatedNodes / 4);
// The pool size is capped at NWORKERS, if target VM value is more than that, set it to NWORKERS.
cappedPoolSize = $x;
\$TargetDedicatedNodes = max(0, min(\$targetVMs, cappedPoolSize));
// Set node deallocation mode - keep nodes active only until tasks finish
\$NodeDeallocationOption = taskcompletion;"""

len_vm(s::String) = 1
len_vm(s::Array{String, 1}) = len(s)
len_vm(s) = throw(ArgumentError("`vm_size` must be a String Array{String, 1}"))

function init_culsterless(nworkers=2; credentials=nothing, vm_size="Standard_E8s_v3",
                                      pool_name="JudiPool", verbose=0, nthreads=4,
                                      auto_scale=true, kw...)
    # Check input
    npool = len_vm(vm_size)
    blob_name = lowercase("$(pool_name)tmp")
    # Update verbosity and parameters
    @eval(AzureClusterlessHPC, global __verbose__ =  Bool($verbose))
    global AzureClusterlessHPC.__params__["_NODE_COUNT_PER_POOL"] = "$(nworkers)"
    global AzureClusterlessHPC.__params__["_POOL_ID"] = "$(pool_name)"
    global AzureClusterlessHPC.__params__["_POOL_COUNT"] = "$(npool)"
    global AzureClusterlessHPC.__params__["_POOL_VM_SIZE"] = vm_size
    global AzureClusterlessHPC.__params__["_OMP_NUM_THREADS"] = "$(nthreads)"
    global AzureClusterlessHPC.__params__["_VERBOSE"] = "$(verbose)"
    global AzureClusterlessHPC.__params__["_BLOB_CONTAINER"] = blob_name

    if !isnothing(credentials)
        # reinit everything
        isfile(credentials) || throw(FileNotFoundError(credentials))
        creds = AzureClusterlessHPC.JSON.parsefile(credentials)
        @eval(AzureClusterlessHPC, global __container__ = $blob_name)
        @eval(AzureClusterlessHPC, global __credentials__ = [$creds])
        @eval(AzureClusterlessHPC, global __resources__ = [[] for i=1:length(__credentials__)])
        @eval(AzureClusterlessHPC, global __clients__ = create_clients(__credentials__, batch=true, blob=true))
    end
    # Create pool with idle autoscale. This will be much more efficient with a defined image rather than docker.
    batch = pyimport("azure.batch")
    container_registry = batch.models.ContainerRegistry(
        registry_server = ENV["ACR_REGISTRY_SERVER"],
        user_name = ENV["ACR_USERNAME"],
        password = ENV["ACR_PASSWORD"]
    )

    create_pool(container_registry=container_registry)

    # Export JUDI on azure
    eval(macroexpand(JUDI4Cloud, quote @batchdef using Distributed, JUDI end))
    include(joinpath(@__DIR__, "batch_defs.jl"))
    include(joinpath(@__DIR__, "modeling.jl"))
end

"""
    finalize_culsterless()

Finalize the clusterless job and deletes all resources (pool, tmp container, jobs)
"""
function finalize_culsterless()
    delete_all_jobs()
    delete_container()
    try delete_pool() catch; nothing end
end

function __init__()
    merge!(AzureClusterlessHPC.__params__, _judi_defaults)
    atexit(finalize_culsterless)
end

end # module


text:
module CommunalHelperDreamFlagSwitchGate

using ..Ahorn, Maple
using Ahorn.CommunalHelper

@mapdef Entity "CommunalHelper/MaxHelpingHand/DreamFlagSwitchGate" DreamFlagSwitchGate(
    x::Integer,
    y::Integer,
    width::Integer=Maple.defaultBlockWidth,
    height::Integer=Maple.defaultBlockHeight,
    featherMode::Bool=false,
    oneUse::Bool=false,
    refillCount::Integer=-1,
    below::Bool=false,
    persistent::Bool=false,
    flag::String="flag_touch_switch",
    icon::String="vanilla",
    inactiveColor::String="5FCDE4",
    activeColor::String="FFFFFF",
    finishColor::String="F141DF",
    shakeTime::Number=0.5,
    moveTime::Number=1.8,
    moveEased::Bool=true,
    allowReturn::Bool=false,
    moveSound::String="event:/game/general/touchswitch_gate_open",
    finishedSound::String="event:/game/general/touchswitch_gate_finish",
)

const placements = Ahorn.PlacementDict()

if isdefined(Ahorn, :MaxHelpingHand)
    placements["Dream Flag Switch Gate (Communal Helper, max480's Helping Hand)"] = Ahorn.EntityPlacement(
        DreamFlagSwitchGate,
        "rectangle",
        Dict{String,Any}(),
        function (entity)
            entity.data["nodes"] = [(
                Int(entity.data["x"]) + Int(entity.data["width"]),
                Int(entity.data["y"]),
            )]
        end,
    )
else
    @warn "Max480's Helping Hand not detected: CommunalHelper+MaxHelpingHand plugins not loaded."
end

function getIconSprite(entity::DreamFlagSwitchGate)
    icon = get(entity.data, "icon", "vanilla")

    iconResource = "objects/switchgate/icon00"
    if icon != "vanilla"
        iconResource = "objects/MaxHelpingHand/flagSwitchGate/$(icon)/icon00"
    end

    return Ahorn.getSprite(iconResource, "Gameplay")
end

Ahorn.nodeLimits(entity::DreamFlagSwitchGate) = 1, 1

Ahorn.minimumSize(entity::DreamFlagSwitchGate) = 16, 16
Ahorn.resizable(entity::DreamFlagSwitchGate) = true, true

Ahorn.editingOrder(entity::DreamFlagSwitchGate) = String[
    "x",
    "y",
    "width",
    "height",
    "flag",
    "inactiveColor",
    "activeColor",
    "finishColor",
    "hitSound",
    "moveSound",
    "finishedSound",
    "shakeTime",
    "moveTime",
]

Ahorn.editingOptions(entity::DreamFlagSwitchGate) =
    Dict{String,Any}("icon" => Ahorn.MaxHelpingHandFlagSwitchGate.bundledIcons)

function Ahorn.selection(entity::DreamFlagSwitchGate)
    x, y = Ahorn.position(entity)
    stopX, stopY = Int.(entity.data["nodes"][1])

    width = Int(get(entity.data, "width", 8))
    height = Int(get(entity.data, "height", 8))

    return [
        Ahorn.Rectangle(x, y, width, height),
        Ahorn.Rectangle(stopX, stopY, width, height),
    ]
end

function Ahorn.renderSelectedAbs(ctx::Ahorn.Cairo.CairoContext, entity::DreamFlagSwitchGate)
    startX, startY = Int(entity.data["x"]), Int(entity.data["y"])
    stopX, stopY = Int.(entity.data["nodes"][1])

    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    renderDreamBlock(ctx, stopX, stopY, width, height, entity.data)

    Ahorn.drawArrow(ctx, startX + width / 2, startY + height / 2, stopX + width / 2, stopY + height / 2, Ahorn.colors.selection_selected_fc, headLength=6)

    iconSprite = getIconSprite(entity)
    Ahorn.drawImage(ctx, iconSprite, stopX + div(width - iconSprite.width, 2), stopY + div(height - iconSprite.height, 2))
end

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::DreamFlagSwitchGate)
    x, y = Ahorn.position(entity)

    width = Int(get(entity.data, "width", 8))
    height = Int(get(entity.data, "height", 8))

    renderDreamBlock(ctx, x, y, width, height, entity.data)

    iconSprite = getIconSprite(entity)
    Ahorn.drawImage(ctx, iconSprite, x + div(width - iconSprite.width, 2), y + div(height - iconSprite.height, 2))
end

end

max_stars_repo_name: vitellaryjr/CommunalHelper

max_issues_repo_name: vitellaryjr/CommunalHelper

max_forks_repo_name: vitellaryjr/CommunalHelper


text:
# Stable Dictionary
# =================

mutable struct StableDict{K, V} <: AbstractDict{K, V}
    slots::Vector{Int}
    keys::Vector{K}
    vals::Vector{V}
    used::Int
    nextidx::Int

    function StableDict{K, V}() where {K, V}
        size = 16
        slots = zeros(Int, size)
        keys = Vector{K}(undef, size)
        vals = Vector{V}(undef, size)
        return new{K,V}(slots, keys, vals, 0, 1)
    end

    function StableDict(dict::StableDict{K, V}) where {K, V}
        copy = StableDict{K, V}()
        for (k, v) in dict
            copy[k] = v
        end
        return copy
    end
end

function StableDict(kvs::Pair{K, V}...) where {K, V}
    dict = StableDict{K, V}()
    for (k, v) in kvs
        dict[k] = v
    end
    return dict
end

function StableDict{K, V}(kvs) where {K, V}
    dict = StableDict{K, V}()
    for (k, v) in kvs
        dict[k] = v
    end
    return dict
end

function StableDict(kvs)
    return StableDict([Pair(k, v) for (k, v) in kvs]...)
end

function StableDict()
    return StableDict{Any, Any}()
end

function Base.copy(dict::StableDict)
    return StableDict(dict)
end

function Base.length(dict::StableDict)
    return dict.used
end

function Base.haskey(dict::StableDict, key)
    _, j = indexes(dict, convert(keytype(dict), key))
    return j > 0
end

function Base.getindex(dict::StableDict, key)
    _, j = indexes(dict, convert(keytype(dict), key))
    if j == 0
        throw(KeyError(key))
    end
    return dict.vals[j]
end

function Base.get!(dict::StableDict, key, default)
    if haskey(dict, key)
        return dict[key]
    end
    val = convert(valtype(dict), default)
    dict[key] = val
    return val
end

function Base.get!(f::Function, dict::StableDict, key)
    if haskey(dict, key)
        return dict[key]
    end
    val = convert(valtype(dict), f())
    dict[key] = val
    return val
end

function Base.setindex!(dict::StableDict, val, key)
    k = convert(keytype(dict), key)
    v = convert(valtype(dict), val)
    @label index
    i, j = indexes(dict, k)
    if j == 0
        if dict.nextidx > lastindex(dict.keys)
            expand!(dict)
            @goto index
        end
        dict.keys[dict.nextidx] = k
        dict.vals[dict.nextidx] = v
        dict.slots[i] = dict.nextidx
        dict.used += 1
        dict.nextidx += 1
    else
        dict.slots[i] = j
        dict.keys[j] = k
        dict.vals[j] = v
    end
    return dict
end

function Base.delete!(dict::StableDict, key)
    k = convert(keytype(dict), key)
    i, j = indexes(dict, k)
    if j > 0
        dict.slots[i] = -j
        dict.used -= 1
    end
    return dict
end

function Base.pop!(dict::StableDict)
    if isempty(dict)
        throw(ArgumentError("empty"))
    end
    i = dict.slots[argmax(dict.slots)]
    key = dict.keys[i]
    val = dict.vals[i]
    delete!(dict, key)
    return key => val
end

function Base.iterate(dict::StableDict)
    if length(dict) == 0
        return nothing
    end
    if dict.used == dict.nextidx - 1
        keys = dict.keys[1:dict.used]
        vals = dict.vals[1:dict.used]
    else
        idx = sort!(dict.slots[dict.slots .> 0])
        @assert length(idx) == length(dict)
        keys = dict.keys[idx]
        vals = dict.vals[idx]
    end
    return (keys[1] => vals[1]), (2, keys, vals)
end

function Base.iterate(dict::StableDict, st)
    i = st[1]
    if i > length(st[2])
        return nothing
    end
    return (st[2][i] => st[3][i]), (i + 1, st[2], st[3])
end

function hashindex(key, sz)
    return (reinterpret(Int, hash(key)) & (sz-1)) + 1
end

function indexes(dict, key)
    sz = length(dict.slots)
    h = hashindex(key, sz)
    i = 0
    while i < sz
        j = mod1(h + i, sz)
        k = dict.slots[j]
        if k == 0
            return j, k
        elseif k > 0 && isequal(dict.keys[k], key)
            return j, k
        end
        i += 1
    end
    return 0, 0
end

function expand!(dict)
    sz = length(dict.slots)
    newsz = sz * 2
    newslots = zeros(Int, newsz)
    resize!(dict.keys, newsz)
    resize!(dict.vals, newsz)
    for i in 1:sz
        j = dict.slots[i]
        if j > 0
            k = hashindex(dict.keys[j], newsz)
            while newslots[mod1(k, newsz)] != 0
                k += 1
            end
            newslots[mod1(k, newsz)] = j
        end
    end
    dict.slots = newslots
    return dict
end


text:
if !isdefined(:TimeTestExecution)
    const TimeTestExecution = false
end

module RegularizedRegressionTests

startclocktime = time()
include("helper.jl")

import Compat.String

my_tests = [
    "test_loss_functions.jl",

    "test_matrix_dataset.jl",

    "test_scad_and_mcp_regression.jl",
    "test_slope_regression.jl",
    "test_glmnet_lasso_and_ridge_regression.jl",
]

if Main.TimeTestExecution

# readstring not available pre 0.5:
if VERSION < v"0.5.0"
    readstring(stream) = readall(stream)
end

function get_git_remote_and_branch()
    lines = split(readstring(`git remote -v show`), "\n")
    remote = match(r"[a-z0-9]+\s+([^\s]+)", lines[1]).captures[1]
    branch = strip(readstring(`git rev-parse --abbrev-ref HEAD`))
    commit = strip(readstring(`git rev-parse HEAD`))
    return remote, branch, commit
end

gitremote, gitbranch, gitcommit = get_git_remote_and_branch()
gitstr = gitremote * "/" * gitbranch * "/" * gitcommit[1:6]
versionstr = string(VERSION)

using DataFrames

TestTimingFileName = "test/timing_testing.csv"

if isfile("test/timing_testing.csv")
    timing_data = readtable("test/timing_testing.csv")
else
    timing_data = DataFrame(TimeStamp = AbstractString[], Julia = AbstractString[], Git = AbstractString[], TestFile = AbstractString[], Elapsed = Float64[])
end

end

using CPUTime

numtestfiles = 0
starttime = CPUtime_us()
@testset "RegularizedRegression test suite" begin

for t in my_tests
    if Main.TimeTestExecution
        CPUtic()

        # Including the test file runs the tests in there...
        include(t)

        elapsed = CPUtoq()
        datestr = Libc.strftime("%Y%m%d %H:%M.%S", time())
        push!(timing_data, [datestr, versionstr, gitstr, t, elapsed])
    else
        include(t)
    end
    numtestfiles += 1
    print("."); flush(STDOUT);
end

end
elapsed = float(CPUtime_us() - starttime)/1e6

if Main.TimeTestExecution
    datestr = Libc.strftime("%Y%m%d %H:%M.%S", time())
    using SHA
    hash = bytes2hex(sha512(join(map(fn -> readstring(open(joinpath("test", fn))), my_tests))))[1:16]
    push!(timing_data, [datestr, versionstr, gitstr, "TOTAL TIME for $(length(my_tests)) test files, $(hash)", elapsed])
    writetable(TestTimingFileName, timing_data)
    println("Wrote $(nrow(timing_data)) rows to file $TestTimingFileName")
end

elapsedclock = time() - startclocktime
println("Tested $(numtestfiles) files in $(round(elapsedclock, 1)) seconds.")

end # module RegRegTests


text:
# Routines to compute surface velocities and their jumps

@inline function surface_velocity_jump!(Δus::VectorData{N},sys::NavierStokes{NX,NY,N},t::Real) where {NX,NY,N}
    assign_velocity!(sys.Vb,sys.bodies,sys.motions,t)
    surface_velocity_jump!(Δus,sys.Vb,sys)
    return Δus
  end

@inline surface_velocity_jump!(Δus::VectorData{N},us::VectorData{N},
                               ::NavierStokes{NX,NY,N,MT,FS,InternalFlow}) where {NX,NY,N,MT,FS} =
                               (Δus .= -us; Δus)

@inline surface_velocity_jump!(Δus::VectorData{N},us::VectorData{N},
                               ::NavierStokes{NX,NY,N,MT,FS,ExternalFlow}) where {NX,NY,N,MT,FS} =
                               (Δus .= us; Δus)

@inline surface_velocity_jump!(Δus::VectorData{N},us::VectorData{N},
                               ::NavierStokes{NX,NY,N,MT,FS,ExternalInternalFlow}) where {NX,NY,N,MT,FS} =
                               (Δus .= 0.0; Δus)

@inline function surface_velocity!(ūs::VectorData{N},sys::NavierStokes{NX,NY,N},t::Real) where {NX,NY,N}
    assign_velocity!(sys.Vb,sys.bodies,sys.motions,t)
    surface_velocity!(ūs,sys.Vb,sys)
    return ūs
  end

@inline surface_velocity!(ūs::VectorData{N},us::VectorData{N},
                               ::NavierStokes{NX,NY,N,MT,FS,ExternalInternalFlow}) where {NX,NY,N,MT,FS} =
                               (ūs .= us; ūs)

@inline surface_velocity!(ūs::VectorData{N},us::VectorData{N},
                               ::NavierStokes{NX,NY,N,MT,FS,SD}) where {NX,NY,N,MT,FS,SD} =
                               (ūs .= 0.5*us; ūs)

@inline function relative_surface_velocity!(ūsr::VectorData{N},sys::NavierStokes{NX,NY,N},t::Real) where {NX,NY,N}
  assign_velocity!(sys.Vb,sys.bodies,sys.motions,t)
  surface_velocity!(ūsr,sys.Vb,sys)
  ūsr .-= sys.Vb
end


text:
function EnumeratedDistribution(arg0::List)
    return EnumeratedDistribution((List,), arg0)
end

function check_and_normalize(::Type{EnumeratedDistribution}, arg0::Vector{jdouble})
    return jcall(EnumeratedDistribution, "checkAndNormalize", Vector{jdouble}, (Vector{jdouble},), arg0)
end

function equals(obj::Object, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_class(obj::Object)
    return jcall(obj, "getClass", Class, ())
end

function get_pmf(obj::EnumeratedDistribution)
    return jcall(obj, "getPmf", List, ())
end

function hash_code(obj::Object)
    return jcall(obj, "hashCode", jint, ())
end

function notify(obj::Object)
    return jcall(obj, "notify", void, ())
end

function notify_all(obj::Object)
    return jcall(obj, "notifyAll", void, ())
end

function probability(obj::EnumeratedDistribution, arg0::Object)
    return jcall(obj, "probability", jdouble, (Object,), arg0)
end

function to_string(obj::Object)
    return jcall(obj, "toString", JString, ())
end

function wait(obj::Object)
    return jcall(obj, "wait", void, ())
end

function wait(obj::Object, arg0::jlong)
    return jcall(obj, "wait", void, (jlong,), arg0)
end

function wait(obj::Object, arg0::jlong, arg1::jint)
    return jcall(obj, "wait", void, (jlong, jint), arg0, arg1)
end


text:
include("load.jl")

using Dates
using TapirSchedulersBenchmarks
suite = TapirSchedulersBenchmarks.BenchFib.setup([30])
bench = suite["N=30"]["ws"]

git_status = read(setenv(`git status`; dir = @__DIR__), String)
git_show = read(setenv(`git show --no-patch`; dir = @__DIR__), String)

stats = @timed for i in 1:1000
    @show (i, now())
    @time run(bench)
end

@info(
    "Finished successfully",
    git_status = Text(git_status),
    git_show = Text(git_show),
    time = Text(string(floor(Int, stats.time ÷ 60), " minutes")),
    VERSION,
    Base.GIT_VERSION_INFO,
)


text:
# replay_buffer.jl
# original: https://github.com/denizyuret/Knet.jl/blob/master/examples/reinforcement-learning/dqn/replay_buffer.jl

import Base.length
import Base.push!

"""
Replay buffer implementation based on https://github.com/openai/baselines/blob/master/baselines/deepq/replay_buffer.py
"""
mutable struct ReplayBuffer
    size::Int
    storage::Vector{Any}
    next_idx::Int
end

"""
Constructor
"""
ReplayBuffer(size::Int) = ReplayBuffer(size, Any[], 1)

"""
Length of the buffer
"""
length(buf::ReplayBuffer) = length(buf.storage)

"""
Tells whether there is enough data in the buffer or not
"""
can_sample(buf::ReplayBuffer, batch_size::Int) = batch_size <= length(buf)

"""
Add new transition to the buffer
"""
function push!(buf::ReplayBuffer, obs_t, action, reward, obs_tp1, done)
    data = (obs_t, action, reward, obs_tp1, done)

    if buf.next_idx > length(buf)
        push!(buf.storage, data)
    else
        buf.storage[buf.next_idx] = data
    end

    buf.next_idx = mod1(buf.next_idx + 1, buf.size)
end

function encode_sample(buf::ReplayBuffer, idxes; stack=1)
    bs = length(idxes)
    inpsize = size(buf.storage[1][1])[1:(end-1)] # get the dimensions of the input
    stacksize = nothing
    if length(inpsize) < 3
        stacksize = (inpsize[1]*stack, inpsize[2:end]...)
    else
        stacksize = (inpsize[1], inpsize[2], inpsize[3]*stack, inpsize[4:end]...)
    end

    obses_t = zeros(Float32, stacksize..., bs)
    actions = zeros(Int, bs)
    rewards = zeros(Float32, bs)
    obses_tp1 = zeros(Float32, stacksize..., bs)
    dones = zeros(Float32, 1, bs)

    indx = 0
    for ind in idxes
        indx += 1
        for i=-(stack-1):0
            curr = ind + i
            curr = curr < 1 ? 1 : curr
            data = buf.storage[curr] #use available data instead of filling the history with 0

            obs_t, action, reward, obs_tp1, done = data
            if length(inpsize) < 3
                t_t = size(obs_t, 1)
                obses_t[t_t*(i-1+stack)+1:t_t*(i-1+stack)+t_t, (Colon() for _=2:ndims(obs_t)-1)..., indx] = obs_t #stack frames
                t_tp1 = size(obs_tp1, 1)
                obses_tp1[t_tp1*(i-1+stack)+1:t_tp1*(i-1+stack)+t_tp1, (Colon() for _=2:ndims(obs_tp1)-1)..., indx] = obs_tp1 #stack frames
            else
                t_t = size(obs_t, 3)
                obses_t[:, :, t_t*(i-1+stack)+1:t_t*(i-1+stack)+t_t, (Colon() for _=4:ndims(obs_t)-1)..., indx] = obs_t #stack frames
                t_tp1 = size(obs_tp1, 3)
                obses_tp1[:, :, t_tp1*(i-1+stack)+1:t_tp1*(i-1+stack)+t_tp1, (Colon() for _=4:ndims(obs_tp1)-1)..., indx] = obs_tp1 #stack frames
            end
            if i == 0
                actions[indx] = action
                rewards[indx] = reward
                dones[1, indx] = done ? 0.0 : 1.0
            end
        end
    end
    obses_t, actions, rewards, obses_tp1, dones
end

function encode_recent(buf::ReplayBuffer, obs_t; stack::Int=1)
    if length(buf) != 0
        recent_idx = max(1, (buf.next_idx - 1) % (buf.size+1))
        obses_t, _, _, _,_ = encode_sample(buf, [recent_idx]; stack=stack)
        inpsize = size(obs_t)[1:(end-1)] # get the dimensions of the input
        i = 0
        if length(inpsize) < 3
            t = size(obs_t, 1)
            obses_t[t*(i-1+stack)+1:t*(i-1+stack)+t, (Colon() for _=2:ndims(obs_t)-1)..., 1] = obs_t #stack frames
        else
            t = size(obs_t, 3)
            obses_t[:, :, t*(i-1+stack)+1:t*(i-1+stack)+t, (Colon() for _=4:ndims(obs_t)-1)..., 1] = obs_t #stack frames
        end
    else
        bs = 1
        inpsize = size(obs_t)[1:end-1] # get the dimensions of the input
        stacksize = nothing
        if length(inpsize) < 3
            stacksize = (inpsize[1]*stack, inpsize[2:end]...)
        else
            stacksize = (inpsize[1], inpsize[2], inpsize[3]*stack, inpsize[4:end]...)
        end

        obses_t = zeros(Float32, stacksize..., bs)
        for i = -(stack-1):0
            if length(inpsize) < 3
                t = size(obs_t, 1)
                obses_t[t*(i-1+stack)+1:t*(i-1+stack)+t, (Colon() for _=2:ndims(obs_t)-1)..., 1] = obs_t #stack frames
            else
                t = size(obs_t, 3)
                obses_t[:, :, t*(i-1+stack)+1:t*(i-1+stack)+t, (Colon() for _=4:ndims(obs_t)-1)..., 1] = obs_t #stack frames
            end
        end
    end
    return obses_t
end

function sample_batch(buf::ReplayBuffer, batchsize::Int; stack::Int=1)
    idxes = randperm(length(buf))[1:batchsize]
    return encode_sample(buf, idxes; stack=stack)
end

max_stars_repo_name: JuliaTagBot/JuliaGymDemo

max_issues_repo_name: JuliaTagBot/JuliaGymDemo

max_forks_repo_name: JuliaTagBot/JuliaGymDemo


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule MPItrampoline_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("MPItrampoline")
JLLWrappers.@generate_main_file("MPItrampoline", UUID("f1f71cc9-e9ae-5b93-9b94-4fe0e1ad3748"))
end  # module MPItrampoline_jll


text:
#const libstemmer = joinpath(dirname(@__FILE__),"..","deps","usr","lib", "libstemmer."*Libdl.dlext)
#@BinDeps.load_dependencies [:libstemmer=>:libstemmer]

##
# character encodings supported by libstemmer
const UTF_8         = "UTF_8"
const ISO_8859_1    = "ISO_8859_1"
const CP850         = "CP850"
const KOI8_R        = "KOI8_R"

##
# lists the stemmer algorithms loaded
function stemmer_types()
    cptr = ccall((:sb_stemmer_list, libstemmer), Ptr{Ptr{UInt8}}, ())
    (C_NULL == cptr) && error("error getting stemmer types")

    stypes = AbstractString[]
    i = 1
    while true
        name_ptr = unsafe_load(cptr, i)
        (C_NULL == name_ptr) && break
        push!(stypes, unsafe_string(name_ptr))
        i += 1
    end
    stypes
end

mutable struct Stemmer
    cptr::Ptr{Cvoid}
    alg::String
    enc::String

    function Stemmer(stemmer_type, charenc=UTF_8)
        cptr = ccall((:sb_stemmer_new, libstemmer),
                    Ptr{Cvoid},
                    (Ptr{UInt8}, Ptr{UInt8}),
                    String(stemmer_type), String(charenc))

        if cptr == C_NULL
            if charenc == UTF_8
                error("stemmer '$(stemmer_type)' is not available")
            else
                error("stemmer '$(stemmer_type)' is not available for encoding '$(charenc)'")
            end
        end

        stm = new(cptr, stemmer_type, charenc)
        finalizer(release, stm)
        stm
    end
end

show(io::IO, stm::Stemmer) = println(io, "Stemmer algorithm:$(stm.alg) encoding:$(stm.enc)")

function release(stm::Stemmer)
    (C_NULL == stm.cptr) && return
    ccall((:sb_stemmer_delete, libstemmer), Cvoid, (Ptr{Cvoid},), stm.cptr)
    stm.cptr = C_NULL
    nothing
end

function stem(stemmer::Stemmer, bstr::AbstractString)
    sres = ccall((:sb_stemmer_stem, libstemmer),
                Ptr{UInt8},
                (Ptr{UInt8}, Ptr{UInt8}, Cint),
                stemmer.cptr, bstr, sizeof(bstr))
    (C_NULL == sres) && error("error in stemming")
    slen = ccall((:sb_stemmer_length, libstemmer), Cint, (Ptr{Cvoid},), stemmer.cptr)
    bytes = unsafe_wrap(Array, sres, Int(slen), own=false)
    String(copy(bytes))
end

function stem_all(stemmer::Stemmer, lang::S, sentence::AbstractString) where S <: Language
    tokens = TextAnalysis.tokenize(lang, sentence)
    stemmed = stem(stemmer, tokens)
    join(stemmed, ' ')
end

function stem(stemmer::Stemmer, words::Array)
    l::Int = length(words)
    ret = Array{String}(undef, l)
    for idx in 1:l
        ret[idx] = stem(stemmer, words[idx])
    end
    ret
end

function stemmer_for_document(d::AbstractDocument)
    Stemmer(lowercase(Languages.english_name(language(d))))
end

function stem!(d::AbstractDocument)
    stemmer = stemmer_for_document(d)
    stem!(stemmer, d)
    release(stemmer)
end

stem!(stemmer::Stemmer, d::FileDocument) = error("FileDocument cannot be modified")

function stem!(stemmer::Stemmer, d::StringDocument)
    stemmer = stemmer_for_document(d)
    d.text = stem_all(stemmer, language(d), d.text)
    nothing
end

function stem!(stemmer::Stemmer, d::TokenDocument)
    d.tokens = stem(stemmer, d.tokens)
    nothing
end

function stem!(stemmer::Stemmer, d::NGramDocument)
    for token in keys(d.ngrams)
        new_token = stem(stemmer, token)
        if new_token != token
            if haskey(d.ngrams, new_token)
                d.ngrams[new_token] = d.ngrams[new_token] + d.ngrams[token]
            else
                d.ngrams[new_token] = d.ngrams[token]
            end
            delete!(d.ngrams, token)
        end
    end
end

function stem!(crps::Corpus)
    stemmer = stemmer_for_document(crps.documents[1])
    for doc in crps
        stem!(stemmer, doc)
    end
    release(stemmer)
end


text:
import Flux

export NormalHyperprior

"""
    Hyperprior

An abstract hyperprior.
"""
abstract type Hyperprior end

"""
    NormalHyperprior

An isotropic multivariate normal hyperprior, parameterized by mean and
element-wise standard deviation.
"""
struct NormalHyperprior{V<:AbstractVector} <: Hyperprior
  mean   :: V
  stddev :: V
end

Flux.trainable(k::NormalHyperprior) = ()
Flux.@functor NormalHyperprior

"""
    (hp::NormalHyperprior)(x::AbstractVector)

Evaluates the normal hyperprior at ``x``.
"""
function (hp::NormalHyperprior)(x::AbstractVector)
  sum(((x .- hp.mean) ./ hp.stddev).^2; dims=1) ./ 2
end


text:
module Day22

export day22

function parseinput(filename)
    deck = Dict{Int,Vector{Int}}()
    player = 0
    for line in readlines(filename)
        if line == ""
            continue
        elseif length(line) >= 6
            player = parse(Int, line[8])
        else
            push!(get!(deck, player, Int[]), parse(Int, line))
        end
    end
    deck
end

function day22(::Val{:a}, filename)
    deck = parseinput(filename)
    looser = maximum(deck[1]) > maximum(deck[2]) ? 2 : 1
    winner = maximum(deck[1]) > maximum(deck[2]) ? 1 : 2
    round = 1
    while !isempty(deck[looser])
        a = deck[1][1]
        deck[1] = deck[1][2:end]
        b = deck[2][1]
        deck[2] = deck[2][2:end]
        if a > b
            push!(deck[1],a)
            push!(deck[1],b)
        else
            push!(deck[2],b)
            push!(deck[2],a)
        end
        round += 1
    end
    foldl((acc, (i,x)) -> acc + i*x, enumerate(reverse(deck[winner])), init = 0)
end

function play(deck)
    seen_decks = Set{Tuple{Int,Vector{Int}}}()
    newdeck = Dict{Int,Vector{Int}}()
    winner = 0
    while length(deck[1]) != 0 && length(deck[2]) != 0

        if (1,deck[1]) in seen_decks || (2,deck[2]) in seen_decks
            return (1,deck)
        else
            a = deck[1][1]
            b = deck[2][1]
            if length(deck[1][2:end]) >= a && length(deck[2][2:end]) >= b
                newdeck[1] = copy(deck[1][2:1+a])
                newdeck[2] = copy(deck[2][2:1+b])
                (winner, ds) = play(newdeck)
            else
                winner = a > b ? 1 : 2
            end
            if winner == 1
                push!(deck[1],a)
                push!(deck[1],b)
            else
                push!(deck[2],b)
                push!(deck[2],a)

            end
            push!(seen_decks, (1, deck[1]))
            push!(seen_decks, (2, deck[2]))

            deck[1] = deck[1][2:end]
            deck[2] = deck[2][2:end]

        end
    end
    return (winner, deck)
end

function day22(::Val{:b}, filename)
    deck = parseinput(filename)

    (winner, deck) = play(deck)

    display(deck)

    foldl((acc, (i,x)) -> acc + i*x, enumerate(reverse(deck[winner])), init = 0)
end

end


text:
using LinearAlgebra: checksquare

"""
    quadratic_expansion(A::IntervalMatrix, α::Real, β::Real)

Compute the quadratic expansion of an interval matrix, ``αA + βA^2``, using
interval arithmetics.

### Input

- `A` -- interval matrix
- `α` -- linear coefficient
- `β` -- quadratic coefficient

### Output

An interval matrix that encloses ``B := αA + βA^2``.

### Algorithm

This a variation of the algorithm in [1, Section 6]. If ``A = (aᵢⱼ)`` and
``B := αA + βA^2 = (bᵢⱼ)``, the idea is to compute each ``bᵢⱼ`` by factoring
out repeated expressions (thus the term *single-use expressions*).

First, let ``i = j``. In this case,

```math
bⱼⱼ = β\\sum_\\{k, k ≠ j} a_{jk} a_{kj} + (α + βa_{jj}) a_{jj}.
```

Now consider ``i ≠ j``. Then,

```math
bᵢⱼ = β\\sum_\\{k, k ≠ i, k ≠ j} a_{ik} a_{kj} + (α + βa_{ii} + βa_{jj}) a_{ij}.
```

[1] Kosheleva, Kreinovich, Mayer, Nguyen. Computing the cube of an interval
matrix is NP-hard. SAC 2005.
"""
function quadratic_expansion(A::IntervalMatrix, α::Real, β::Real)
    B = similar(A)
    n = checksquare(A)

    # case i = j
    @inbounds for j in 1:n
        B[j, j] = quadratic_expansion(A[j, j], α, β)
        for k in 1:n
            k == j && continue
            B[j, j] += β * (A[j, k] * A[k, j])
        end
    end

    # case i ≠ j
    @inbounds for j in 1:n
        for i in 1:n
            i == j && continue
            B[i, j] = A[i, j] * (α + β*(A[j, j] + A[i, i]))
            for k in 1:n
                (k == i || k == j) && continue
                B[i, j] += β * (A[i, k] * A[k, j])
            end
        end
    end
    return B
end

function quadratic_expansion(x::Interval, α::Real, β::Real)
    iszero(β) && return α * x

    return ((2 * β * x + α) ^ 2 - α ^ 2) / (4 * β)
end

function _truncated_exponential_series(A::IntervalMatrix{T}, t, p::Integer;
                                       n=checksquare(A)) where {T}
    if p == 0
        # index i = 0 (identity matrix)
        return IntervalMatrix(Interval(one(T)) * I, n)
    elseif p == 1
        # index i = 1
        S = A * t
    else
        # indices i = 1 and i = 2
        S = quadratic_expansion(A, t, t^2/2)
    end

    # index i = 0, (identity matrix, added implicitly)
    for i in 1:n
        S[i, i] += one(T)
    end

    if p < 3
        return S
    end

    # indices i >= 3
    pow = IntervalMatrixPower(A)
    increment!(pow)
    fact_num = t^2
    fact_denom = 2
    for i in 3:p
        fact_num *= t
        fact_denom *= i
        Aⁱ = increment!(pow)
        S += Aⁱ * (fact_num / fact_denom)
    end

    return S
end

"""
    exp_overapproximation(A::IntervalMatrix{T, Interval{T}}, t, p) where {T}

Overapproximation of the exponential of an interval matrix.

### Input

- `A` -- interval matrix
- `t` -- non-negative time value
- `p` -- order of the approximation

### Algorithm

See Theorem 1 in *Reachability Analysis of Linear Systems with Uncertain
Parameters and Inputs* by M. Althoff, O. Stursberg, M. Buss.
"""
function exp_overapproximation(A::IntervalMatrix{T, Interval{T}}, t, p) where {T}
    n = checksquare(A)
    S = _truncated_exponential_series(A, t, p; n=n)
    E = _exp_remainder(A, t, p; n=n)
    return S + E
end

# Implementation of Prop. 1 in Althoff, Matthias, Bruce H. Krogh, and Olaf Stursberg.
# "Analyzing reachability of linear dynamic systems with parametric uncertainties."
# Modeling, Design, and Simulation of Systems with Uncertainties. Springer, Berlin, Heidelberg, 2011. 69-94.
function _exp_remainder(A::IntervalMatrix{T}, t, p; n=checksquare(A)) where {T}
    C = max.(abs.(inf(A)), abs.(sup(A)))
    # compute Q = I + Ct + (Ct)^2/2! + ... + (Ct)^p/p!
    Q = Matrix(Diagonal(ones(T, n)))

    tⁱ = 1
    i! = 1
    Cⁱ = copy(C)
    for i in 1:p
        i! *= i
        tⁱ *= t
        Q += Cⁱ * tⁱ/i!
        Cⁱ *= C
    end
    M = exp(C*t)
    Y = M - Q
    Γ = IntervalMatrix(fill(zero(T)±one(T), (n, n)))
    E = Γ * Y
    return E
end

# Estimates the sum of the series in the matrix exponential. See Theorem 1
# in [1] Althoff, Matthias, Olaf Stursberg, and Martin Buss.
# Reachability analysis of nonlinear systems with uncertain parameters using conservative linearization.
# 2008 47th IEEE Conference on Decision and Control. IEEE, 2008.
function _exp_remainder_series(A::IntervalMatrix{T}, t, p; n=checksquare(A)) where {T}
    nA = opnorm(A, Inf)
    c = nA * t / (p + 2)
    @assert c < 1 "the remainder of the matrix exponential could not be " *
        "computed because a convergence condition is not satisfied: $c ≥ 1 " *
        "but it should be smaller than 1; try choosing a larger order"
    Γ = IntervalMatrix(fill(zero(T)±one(T), (n , n)))
    return Γ * ((nA*t)^(p+1) * (1/factorial(p + 1) * 1/(1-c)))
end

"""
    horner(A::IntervalMatrix{T}, K::Integer; [validate]::Bool=true)

Compute the matrix exponential using the Horner scheme.

### Input

- `A` -- interval matrix
- `K` -- number of expansions in the Horner scheme
- `validate` -- (optional; default: `true`) option to validate the precondition
                of the algorithm

### Algorithm

We use the algorithm in [1, Section 4.2].

[1] Goldsztejn, Alexandre, Arnold Neumaier. "On the exponentiation of interval
matrices". Reliable Computing. 2014.
"""
function horner(A::IntervalMatrix{T}, K::Integer;
                       validate::Bool=true) where {T}
    if validate
        nA = opnorm(A, Inf)
        c = K + 2
        if c <= nA
            throw(ArgumentError("the precondition for the " *
                "Horner-scheme algorithm is not satisfied: $c <= $nA; " *
                "try choosing a larger order"))
        end
    end
    if K <= 0
        throw(ArgumentError("the Horner evaluation requires a positive " *
            "number of expansions but received $K"))
    end

    n = checksquare(A)
    Iₙ = IntervalMatrix(Interval(one(T)) * I, n)
    H = Iₙ + A/K
    for i in (K-1):-1:1
        H = Iₙ + A / i * H
    end

    # remainder; the paper uses a less precise computation here
    R = _exp_remainder(A, one(T), K)

    return H + R
end

"""
    scale_and_square(A::IntervalMatrix{T}, l::Integer, t, p;
                     [validate]::Bool=true)

Compute the matrix exponential using scaling and squaring.

### Input

- `A` -- interval matrix
- `l` -- scaling-and-squaring order
- `t` -- non-negative time value
- `p` -- order of the approximation
- `validate` -- (optional; default: `true`) option to validate the precondition
                of the algorithm

### Algorithm

We use the algorithm in [1, Section 4.3], which first scales `A` by factor
``2^{-l}``, computes the matrix exponential for the scaled matrix, and then
squares the result ``l`` times.

```math
    \\exp(A * 2^{-l})^{2^l}
```

[1] Goldsztejn, Alexandre, Arnold Neumaier. "On the exponentiation of interval
matrices". Reliable Computing. 2014.
"""
function scale_and_square(A::IntervalMatrix{T}, l::Integer, t, p;
                          validate::Bool=true) where {T}
    if validate
        nA = opnorm(A, Inf) * t
        c = (p + 2) * 2.0^l
        if c <= nA
            throw(ArgumentError("the precondition for the " *
                "scaling-and-squaring algorithm is not satisfied: $c <= $nA; " *
                "try choosing a larger order"))
        end
    end

    A_scaled = A / Interval(T(2))^l
    E = exp_overapproximation(A_scaled, t, p)
    for i in 1:l
        E = square(E)
    end
    return E
end

"""
    exp_underapproximation(M::IntervalMatrix{T, Interval{T}}, t, p) where {T}

Overapproximation of the exponential of an interval matrix.

### Input

- `A` -- interval matrix
- `t` -- non-negative time value
- `p` -- order of the approximation

### Algorithm

See Theorem 2 in *Reachability Analysis of Linear Systems with Uncertain
Parameters and Inputs* by M. Althoff, O. Stursberg, M. Buss.
"""
function exp_underapproximation(A::IntervalMatrix{T, Interval{T}}, t, p) where {T}
    @assert p > 1 "the order $p < 2 is not supported"
    n = checksquare(A)

    Y = zeros(n, n)
    LA = inf(A)
    Aⁱl = LA^2
    Z = zeros(n, n)
    RA = sup(A)
    Aⁱr = RA^2
    fact_num = t^2
    fact_denom = 2
    for i in 3:p
        fact_num *= t
        fact_denom *= i
        fact = fact_num / fact_denom
        Aⁱl *= LA
        Y += Aⁱl * fact
        Aⁱr *= RA
        Z += Aⁱr * fact
    end

    B = IntervalMatrix{T}(undef, n , n)
    @inbounds for j in 1:n
        for i in 1:n
            minYZ = min(Y[i, j], Z[i, j])
            maxYZ = max(Y[i, j], Z[i, j])
            B[i, j] = Interval(minYZ, maxYZ)
        end
    end

    W = quadratic_expansion(A, t, t^2/2)
    res = W + B

    # add identity matrix implicitly
    for i in 1:n
        @inbounds res[i, i] += one(T)
    end
    return res
end


text:
timeline, modules = replay("/home/kfischer/.local/share/rr/llc-43/")
for i = 1:2000
RR.step!(current_session(timeline))
icxx"$timeline->maybe_add_reverse_exec_checkpoint(rr::ReplayTimeline::LOW_OVERHEAD);"
end

using Gallium: RemotePtr

h,base,sym = Gallium.lookup_sym(timeline, modules, :malloc)
hook_addr = Gallium.RemoteCodePtr(base + ObjFileBase.deref(sym).st_value)
jit, callbacks = create_remote_jit(timeline, hook_addr)
icxx"$callbacks->session = $(pointer_from_objref(timeline));"

task = current_task(current_session(timeline))

# Ok, now allocate space (2GB) for the data buffer
data_buffer_start = 0x70001000
data_buffer_size = 2*1024*1024*1024
icxx"""
rr::AutoRemoteSyscalls remote($(current_task(current_session(timeline))));
remote.infallible_mmap_syscall($data_buffer_start, $data_buffer_size,
    PROT_READ | PROT_WRITE,
    MAP_ANONYMOUS | MAP_PRIVATE | MAP_FIXED, -1, 0);
"""
icxx"""
rr::AutoRemoteSyscalls remote($(current_task(current_session(timeline))));
remote.infallible_mmap_syscall($data_buffer_start+$data_buffer_size,0x1000,
    PROT_NONE,
    MAP_ANONYMOUS | MAP_PRIVATE | MAP_FIXED, -1, 0);
"""

let __current_compiler__ = TargetClang
    cxx"""
    #include <stdint.h>
    #include "/home/kfischer/rr/src/preload/instrument_interface.h"
    extern "C" {
        static inline uint64_t native_read_pmc(uint32_t counter)
        {
            uint32_t low, high;
        	asm volatile("rdpmc" : "=a" (low), "=d" (high) : "c" (counter));
        	return low | ((uint64_t)high) << 32;
        }
        extern uint64_t data_buffer_start;
        void malloc_hook() __attribute__((preserve_all)) {
            (&data_buffer_start)[++data_buffer_start] =
                rr_pmc2ticks(native_read_pmc(0));
        }

        void free_hook(uint64_t addr) __attribute__((preserve_all)) {
            (&data_buffer_start)[++data_buffer_start] =
                addr | ((uint64_t)1 << 63); // tag so we now this is a free record
        }

        void malloc_ret_hook() __attribute__((preserve_all)) {
            register uint64_t buffer asm("rax");
            (&data_buffer_start)[++data_buffer_start] = buffer;
        }
    }
    """
end

trace_func(jit, callbacks, :malloc, "malloc_hook", "malloc_ret_hook")
trace_func(jit, callbacks, :free, "free_hook")
#bp1 = Gallium.breakpoint(timeline, :malloc)
#bp2 = Gallium.breakpoint(timeline, :free)

#=
    bp = Gallium.breakpoint(timeline, :malloc)
    bp = Gallium.breakpoint(timeline, exit_hook)
=#

#=
replacement = [hook_template; zeros(UInt8,nbytes-length(hook_template))]
Gallium.store!(task, RemotePtr{UInt8}(hook_addr), replacement)
Gallium.store!(task, RemotePtr{UInt8}(buffer_start_addr),
    Gallium.Hooking.hook_tail_template([code;orig_bytes[1:nbytes]],UInt(hook_addr)+nbytes)
)
disassemble(orig_bytes[1:nbytes])

bp = Gallium.breakpoint(timeline, :malloc)
=#

#=
stacktraces = UInt64[]
using Gallium: breakpoint, conditional
cache = Gallium.Unwinder.CFICache(100_000)
bp = conditional(breakpoint(timeline, :malloc)) do loc, RC
    bt = Gallium.rec_backtrace(RC, current_vm(loc.vm), modules, true, cache) do RC
        push!(stacktraces,UInt(Gallium.ip(RC)))
        return true
    end
    push!(stacktraces, UInt(0))
    return false
end
ASTInterpreter.execute_command(nothing, nothing, Val{:mark}(), "mark")
ASTInterpreter.execute_command(nothing, nothing, Val{:timejump}(), "timejump 128000000")
=#

#=
cache = Gallium.Unwinder.CFICache(100_000)
Profile.init(n=10^9)
ASTInterpreter.execute_command(nothing, nothing, Val{:timejump}(), "timejump @1")
@profile ASTInterpreter.execute_command(nothing, nothing, Val{:timejump}(), "timejump 128000000")
IJulia=1
using ProfileView
ProfileView.svgwrite("profile.svg", C = true)

stacktraces = UInt64[]
ASTInterpreter.execute_command(nothing, nothing, Val{:timejump}(), "timejump @1")
@profile ASTInterpreter.execute_command(nothing, nothing, Val{:timejump}(), "timejump 22001343342")
=#
#=
@load "traces.jl"

symb_map = Dict(ip => Gallium.Unwinder.symbolicate(modules,ip) for ip in
    filter(x->x!=0,unique(stacktraces)))
callframes = UInt64[]
for idx in find(x->x==0, stacktraces)
    idx+2 <= length(stacktraces) || continue
    push!(callframes, stacktraces[idx+2])
end
using StatsBase
map(x->(demangle(symb_map[x[1]]),x[2]),sort(collect(countmap(callframes)), by=x->x[2]))
=#


text:
using Test, DSAModels

function read_dsa_data(filename)
    f = open(filename, "r")
    N_mach = read(f, Int64)
    M = read!(f, Vector{Float64}(undef, N_mach))
    η = read!(f, Vector{Float64}(undef, N_mach))
    η_r = read!(f, Vector{Float64}(undef, N_mach))
    close(f)

    return N_mach, M, η, η_r
end

@testset "DSAModels" begin

    @testset "Mach Models" begin

        @testset "Kang et. al. 2007" begin
            η_model = Kang07()
            N_mach, M, η, η_r = read_dsa_data("Kang07.dat")

            @testset "η" begin

                @test η_Ms(η_model, 0.5, 0.0) == 0.0

                for i = 1:N_mach
                    @test η_Ms(η_model, M[i], 0.0) ≈ η[i]
                end
            end

            @testset "η_r" begin
                @test η_Ms(η_model, 0.5, 1.0) == 0.0

                for i = 1:N_mach
                    @test η_Ms(η_model, M[i], 1.0) ≈ η_r[i]
                end
            end

        end

        @testset "Kang & Ryu 2013" begin
            η_model = KR13()
            N_mach, M, η, η_r = read_dsa_data("KR13.dat")

            @testset "η" begin
                for i = 1:N_mach
                    @test η_Ms(η_model, M[i], 0.0) ≈ η[i]
                end
            end

            @testset "η_r" begin
                for i = 1:N_mach
                    @test η_Ms(η_model, M[i], 1.0) ≈ η_r[i]
                end
            end

        end

        @testset "Caprioli & Spitkovsky 2014" begin
            η_model = CS14()
            N_mach, M, η, η_r = read_dsa_data("CS14.dat")

            @testset "η" begin
                for i = 1:N_mach
                    @test η_Ms(η_model, M[i], 0.0) ≈ η[i]
                end
            end

            @testset "η_r" begin
                for i = 1:N_mach
                    @test η_Ms(η_model, M[i], 1.0) ≈ η_r[i]
                end
            end

        end

        @testset "Ryu et. al. 2019" begin
            η_model = Ryu19()
            N_mach, M, η, η_r = read_dsa_data("Ryu19.dat")

            @testset "η" begin
                for i = 1:N_mach
                    @test η_Ms(η_model, M[i], 0.0) ≈ η[i]
                end
            end

            @testset "η_r" begin
                for i = 1:N_mach
                    @test η_Ms(η_model, M[i], 1.0) ≈ η_r[i]
                end
            end

        end

        @testset "Pfrommer et. al. 2016" begin
            η_model = P16()

            @test η_Ms(η_model, 10.0, 0.0) == 0.5
            @test η_Ms(η_model, 10.0, 1.0) == 0.5
        end
    end

    @testset "B-field Models" begin

        @testset "Electrons" begin
            # read reference data
            N_θ, θ_B, η, η_r = read_dsa_data("theta_B_e.dat")

            @testset "η" begin
                for i = 1:N_θ
                    @test ηB_acc_e(θ_B[i]) ≈ η[i]
                end
            end
            @testset "η_r" begin
                for i = 1:N_θ
                    @test ηB_reacc_e(θ_B[i]) ≈ η_r[i]
                end
            end
        end

        @testset "Protons" begin
            # read reference data
            N_θ, θ_B, η, η_r = read_dsa_data("theta_B_p.dat")

            @testset "η" begin
                for i = 1:N_θ
                    @test ηB_acc_p(θ_B[i]) ≈ η[i]
                end
            end
            @testset "η_r" begin
                for i = 1:N_θ
                    @test ηB_reacc_p(θ_B[i]) ≈ η_r[i]
                end
            end
        end
    end

end


text:
using Distributions
using Random
using LinearAlgebra
using Test
using FiniteDifferences

@testset "LKJCholesky" begin
    function test_draw(d::LKJCholesky, x; check_uplo=true)
        @test insupport(d, x)
        check_uplo && @test x.uplo == d.uplo
    end
    function test_draws(d::LKJCholesky, xs; check_uplo=true, nkstests=1)
        @test all(x -> insupport(d, x), xs)
        check_uplo && @test all(x -> x.uplo == d.uplo, xs)

        p = d.d
        dmat = LKJ(p, d.η)
        marginal = Distributions._marginal(dmat)
        ndraws = length(xs)
        zs = Array{eltype(d)}(undef, p, p, ndraws)
        for k in 1:ndraws
            zs[:, :, k] = Matrix(xs[k])
        end

        @testset "LKJCholesky marginal moments" begin
            @test mean(zs; dims=3)[:, :, 1] ≈ I atol=0.1
            @test var(zs; dims=3)[:, :, 1] ≈ var(marginal) * (ones(p, p) - I) atol=0.1
            @testset for n in 2:5
                for i in 1:p, j in 1:(i-1)
                    @test moment(zs[i, j, :], n) ≈ moment(rand(marginal, ndraws), n) atol=0.1
                end
            end
        end

        @testset "LKJCholesky marginal KS test" begin
            α = 0.01
            L = sum(1:(p - 1))
            for i in 1:p, j in 1:(i-1)
                @test pvalue_kolmogorovsmirnoff(zs[i, j, :], marginal) >= α / L / nkstests
            end
        end
    end

    # Compute logdetjac of ϕ: L → L L' where only strict lower triangle of L and L L' are unique
    function cholesky_inverse_logdetjac(L)
        size(L, 1) == 1 && return 0.0
        J = jacobian(central_fdm(5, 1), cholesky_vec_to_corr_vec, stricttril_to_vec(L))[1]
        return logabsdet(J)[1]
    end
    stricttril_to_vec(L) = [L[i, j] for i in axes(L, 1) for j in 1:(i - 1)]
    function vec_to_stricttril(l)
        n = length(l)
        p = Int((1 + sqrt(8n + 1)) / 2)
        L = similar(l, p, p)
        fill!(L, 0)
        k = 1
        for i in 1:p, j in 1:(i - 1)
            L[i, j] = l[k]
            k += 1
        end
        return L
    end
    function cholesky_vec_to_corr_vec(l)
        L = vec_to_stricttril(l)
        for i in axes(L, 1)
            w = view(L, i, 1:(i-1))
            wnorm = norm(w)
            if wnorm > 1
                w ./= wnorm
                wnorm = 1
            end
            L[i, i] = sqrt(1 - wnorm^2)
        end
        return stricttril_to_vec(L * L')
    end

    @testset "Constructors" begin
        @testset for p in (4, 5), η in (2, 3.5)
            d = LKJCholesky(p, η)
            @test d.d == p
            @test d.η == η
            @test d.uplo == 'L'
            @test d.logc0 == LKJ(p, η).logc0
        end

        @test_throws DomainError LKJCholesky(0, 2)
        @test_throws DomainError LKJCholesky(4, 0.0)
        @test_throws DomainError LKJCholesky(4, -1)

        for uplo in (:U, 'U')
            d = LKJCholesky(4, 2, uplo)
            @test d.uplo === 'U'
        end
        for uplo in (:L, 'L')
            d = LKJCholesky(4, 2, uplo)
            @test d.uplo === 'L'
        end
        @test_throws ArgumentError LKJCholesky(4, 2, :F)
        @test_throws ArgumentError LKJCholesky(4, 2, 'F')
    end

    @testset "REPL display" begin
        d = LKJCholesky(5, 1)
        @test sprint(show, d) == "$(typeof(d))(\nd: 5\nη: 1.0\nuplo: L\n)\n"
    end

    @testset "Conversion" begin
        d = LKJCholesky(5, 3.5)
        df0_1 = convert(LKJCholesky{Float32}, d)
        @test df0_1 isa LKJCholesky{Float32}
        @test df0_1.d == d.d
        @test df0_1.η ≈ d.η
        @test df0_1.uplo == d.uplo
        @test df0_1.logc0 ≈ d.logc0

        df0_2 = convert(LKJCholesky{BigFloat}, d.d, d.η, d.uplo, d.logc0)
        @test df0_2 isa LKJCholesky{BigFloat}
        @test df0_2.d == d.d
        @test df0_2.η ≈ d.η
        @test df0_2.uplo == d.uplo
        @test df0_2.logc0 ≈ d.logc0
    end

    @testset "properties" begin
        @testset for p in (4, 5), η in (2, 3.5), uplo in ('L', 'U')
            d = LKJCholesky(p, η, uplo)
            @test d.d == p
            @test size(d) == (p, p)
            @test Distributions.params(d) == (d.d, d.η, d.uplo)
            @test partype(d) <: Float64

            m = mode(d)
            @test m isa Cholesky{eltype(d)}
            @test Matrix(m) ≈ I
        end
        @test_broken partype(LKJCholesky(2, 4f0)) <: Float32

        @testset "insupport" begin
            @test insupport(LKJCholesky(40, 2, 'U'), cholesky(rand(LKJ(40, 2))))
            @test insupport(LKJCholesky(40, 2), cholesky(rand(LKJ(40, 2))))
            @test !insupport(LKJCholesky(40, 2), cholesky(rand(LKJ(41, 2))))
            z = rand(LKJ(40, 1))
            z .+= exp(Symmetric(randn(size(z)))) .* 1e-8
            x = cholesky(z)
            @test !insupport(LKJCholesky(4, 2), x)
        end
    end

    @testset "Evaluation" begin
        @testset for p in (1, 4, 10), η in (0.5, 1, 3)
            d = LKJ(p, η)
            dchol = LKJCholesky(p, η)
            z = rand(d)
            x = cholesky(z)
            x_L = typeof(x)(Matrix(x.L), 'L', x.info)
            logdetJ = sum(i -> (p - i) * log(x.UL[i, i]), 1:p)
            logdetJ_approx = cholesky_inverse_logdetjac(x.L)
            @test logdetJ ≈ logdetJ_approx

            @test logpdf(dchol, x) ≈ logpdf(d, z) + logdetJ
            @test logpdf(dchol, x_L) ≈ logpdf(dchol, x)

            @test pdf(dchol, x) ≈ exp(logpdf(dchol, x))
            @test pdf(dchol, x_L) ≈ pdf(dchol, x)

            @test loglikelihood(dchol, x) ≈ logpdf(dchol, x)
            xs = cholesky.(rand(d, 10))
            @test loglikelihood(dchol, xs) ≈ sum(logpdf(dchol, x) for x in xs)
        end
    end

    @testset "Sampling" begin
        rng = MersenneTwister(66)
        nkstests = 4 # use for appropriate Bonferroni correction for KS test

        @testset "rand" begin
            @testset for p in (2, 4, 10), η in (0.5, 1, 3), uplo in ('L', 'U')
                d = LKJCholesky(p, η, uplo)
                test_draw(d, rand(rng, d))
                test_draws(d, rand(rng, d, 10^4); nkstests=nkstests)
            end
            @test_broken rand(rng, LKJCholesky(5, Inf)) ≈ I
        end

        @testset "rand!" begin
            @testset for p in (2, 4, 10), η in (0.5, 1, 3), uplo in ('L', 'U')
                d = LKJCholesky(p, η, uplo)
                x = Cholesky(Matrix{Float64}(undef, p, p), uplo, 0)
                rand!(rng, d, x)
                test_draw(d, x)
                x = Cholesky(Matrix{Float64}(undef, p, p), uplo, 0)
                rand!(d, x)
                test_draw(d, x)

                # test that uplo of Cholesky object is respected
                x2 = Cholesky(Matrix{Float64}(undef, p, p), uplo == 'L' ? 'U' : 'L', 0)
                rand!(rng, d, x2)
                test_draw(d, x2; check_uplo = false)

                # allocating
                xs = Vector{typeof(x)}(undef, 10^4)
                rand!(rng, d, xs)
                test_draws(d, xs; nkstests=nkstests)

                F2 = cholesky(exp(Symmetric(randn(rng, p, p))))
                xs2 = [deepcopy(F2) for _ in 1:10^4]
                xs2[1] = cholesky(exp(Symmetric(randn(rng, p + 1, p + 1))))
                rand!(rng, d, xs2)
                test_draws(d, xs2; nkstests=nkstests)

                # non-allocating
                F3 = cholesky(exp(Symmetric(randn(rng, p, p))))
                xs3 = [deepcopy(F3) for _ in 1:10^4]
                rand!(rng, d, xs3)
                test_draws(d, xs3; check_uplo = uplo == 'U', nkstests=nkstests)
            end
        end
    end
end


text:
using Tokenize.Tokens: GREATER, LESS, GREATER_EQ, GREATER_THAN_OR_EQUAL_TO, LESS_EQ, LESS_THAN_OR_EQUAL_TO
using CSTParser: MacroName
begin
    struct ObsoleteVersionCheck; vers; end
    register(ObsoleteVersionCheck, Deprecation(
        "This version check is for a version of julia that is no longer supported by this package",
        "julia",
        typemin(VersionNumber), typemin(VersionNumber), typemax(VersionNumber)
    ))
    ObsoleteVersionCheck() = ObsoleteVersionCheck(Pkg.Reqs.parse(IOBuffer("julia 0.6"))["julia"])

    const comparisons = Dict(
         GREATER                  => >,
         LESS                     => <,
         GREATER_EQ               => >=,
         GREATER_THAN_OR_EQUAL_TO => ≥,
         LESS_EQ                  => <=,
         LESS_THAN_OR_EQUAL_TO    => ≤,
         #= TODO:
         EQEQ                     => ==,
         EQEQEQ                   => ===,
         IDENTICAL_TO             => ≡,
         NOT_EQ                   => !=,
         NOT_EQUAL_TO             => ≠,
         NOT_IS                   => !==
         =#)

    function detect_ver_arguments(VERSION_arg, v_arg)
        isexpr(VERSION_arg, CSTParser.IDENTIFIER) || return nothing
        VERSION_arg.val == "VERSION" || return nothing
        isexpr(v_arg, CSTParser.x_Str) || return nothing
        isexpr(children(v_arg)[1], CSTParser.IDENTIFIER) || return nothing
        isexpr(children(v_arg)[2], CSTParser.LITERAL, Tokens.STRING) || return nothing
        children(v_arg)[1].val == "v" || return nothing
        VersionNumber(children(v_arg)[2].val)
    end

    function dep_for_vers(::Type{ObsoleteVersionCheck}, vers)
        ObsoleteVersionCheck(vers["julia"])
    end

    opcode(x::CSTParser.OPERATOR) = x.kind
    iscomparison(x::CSTParser.OPERATOR) = CSTParser.precedence(x) == 6
    iscomparison(x::OverlayNode) = iscomparison(x.expr)

    function process_comparison(comparison, dep, context)
        isexpr(comparison, CSTParser.BinaryOpCall) || return nothing
        iscomparison(children(comparison)[2]) || return nothing
        comparison = comparison.expr
        opc = opcode(children(comparison)[2])
        haskey(comparisons, opc) || return nothing
        r1 = detect_ver_arguments(children(comparison)[1], children(comparison)[3])
        if r1 !== nothing
            f = comparisons[opc]
            alwaystrue = all(interval->(f(interval.lower, r1) && f(interval.upper, r1)), dep.vers.intervals)
            alwaysfalse = all(interval->(!f(interval.lower, r1) && !f(interval.upper, r1)), dep.vers.intervals)
            @assert !(alwaystrue && alwaysfalse)
            return (alwaystrue, alwaysfalse)
        end
        r2 = detect_ver_arguments(children(comparison)[3], children(comparison)[1])
        if r2 !== nothing
            f = comparisons[opc]
            alwaystrue = all(interval->(f(interval.lower, r2) && f(interval.upper, r2)), dep.vers.intervals)
            alwaysfalse = all(interval->(!f(interval.lower, r2) && !f(interval.upper, r2)), dep.vers.intervals)
            @assert !(alwaystrue && alwaysfalse)
            return (alwaystrue, alwaysfalse)
        end
        return nothing
    end

    function in_statement_position(expr)
        p = parent(expr)
        isexpr(p, CSTParser.FileH) && return true
        if isexpr(p, CSTParser.Begin) || isexpr(p, CSTParser.Block)
            # Could still be in value position if we're the last value in the
            # block and the block itself is in value position
            nextsibling(expr) !== nothing && return true
            return in_statement_position(p)
        end
        return false
    end

    match(ObsoleteVersionCheck, CSTParser.BinaryOpCall) do x
        dep, expr, resolutions, context = x
        result = process_comparison(expr, dep, context)
        result === nothing && return
        (alwaystrue, alwaysfalse) = result
        (alwaysfalse || alwaystrue) || return
        resolve_boolean(dep, resolutions, context, expr, alwaystrue)
    end
end


text:
module StorageManager

export StorageGroup

using ..GraphicsManager

API = nothing

type StorageGroup
	id::Symbol
	data::AbstractGraphicsData
	list::Dict{Symbol,StorageGroup}
	StorageGroup() = new(:_, nothing, Dict())
	StorageGroup(typ::Symbol, subtyp::Symbol,id=:_) = new(id, API.create(typ,subtyp),Dict())
end

function init()
	global API = GRAPHICSDRIVER().StorageManager
end

create(this::StorageGroup, id::Symbol, typ::Symbol, subtyp::Symbol) = haskey(this.list, id) ? this.list[id] : this.list[id]=StorageGroup(typ, subtyp, id)
link(this::StorageGroup, child::StorageGroup, id::Symbol) = (child.id=id;this.list[id]=child)

clean(this::StorageGroup) = API.clean(this.data)

getGroup(this::StorageGroup, id::Symbol) = haskey(this.list, id) ? this.list[id] : nothing

getData(this::StorageGroup, id::Symbol) = haskey(this.list, id) ? this.list[id].data : nothing
function getData(this::StorageGroup, id::Integer)
	c=collect(keys(this.list))
	v=length(c)>=id ? get(this.list,c[id],nothing) : nothing
	if v != nothing v=v.data end
	v
end

getValues(this::StorageGroup) = API.getValues(this.data)
setValues(this::StorageGroup, values::AbstractArray, elems::Integer, mode=:TRIANGLES) = API.setValues(this.data, values, elems, mode)

function link(this::StorageGroup, on=true)
	for (k,block) in this.list
		for (k,storage) in block.list
			link(block, storage, on)
		end
	end
end

link(block::StorageGroup, storage::StorageGroup, on=true) = API.bind(block.data,storage.data,on)

upload(block::StorageGroup, storage::StorageGroup, data::StorageGroup) = API.upload(block.data, storage.data, data.data)

function upload(this::StorageGroup)
	uploads=[]
	for (kb,block) in this.list
		API.prepare(block.data)
		for (ks,storage) in block.list
			API.prepare(block.data, storage.data)
			for (kd,data) in storage.list
				API.prepare(block.data, storage.data, data.data)
				push!(uploads,[block.data, storage.data, data.data])
			end
		end
		API.init(block.data)
	end

	API.upload(uploads)
end

end #StorageManager

max_stars_repo_name: Gilga/JGE

max_issues_repo_name: Gilga/JGE

max_forks_repo_name: Gilga/JGE


text:
export device, send_to_host, send_to_device

using Flux
using CuArrays
using Adapt

send_to_host(x) = send_to_device(Val(:cpu), x)
send_to_device(::Val{:cpu}, x) = x  # cpu(x) is not very efficient! So by default we do nothing here.

send_to_device(::Val{:cpu}, x::CuArray) = adapt(Array, x)
send_to_device(::Val{:gpu}, x) = Flux.fmap(a -> adapt(CuArray{Float32}, a), x)
send_to_device(::Val{:gpu}, x::SubArray{T,N,<:CircularArrayBuffer}) where {T,N} =
    CuArray{T}(x)

"""
    device(model)

Detect the suitable running device for the `model`.
Return `Val(:cpu)` by default.
"""
device(x) = device(Flux.trainable(x))
device(x::Function) = nothing
device(::CuArray) = Val(:gpu)
device(::Array) = Val(:cpu)
device(x::Tuple{}) = nothing
device(x::NamedTuple{(),Tuple{}}) = nothing

function device(x::Union{Tuple,NamedTuple})
    d1 = device(x[1])
    if isnothing(d1)
        device(Base.tail(x))
    else
        d1
    end
end


text:
mutable struct FieldHandler
    fields::Vector{Field}
    cellset::Set{Int}
    bc_values::Vector{BCValues} # for boundary conditions
end

function FieldHandler(fields::Vector{Field}, cellset)
    # TODO for now, only accept isoparamtric mapping
    bc_values = [BCValues(field.interpolation, field.interpolation) for field in fields]
    FieldHandler(fields, cellset, bc_values)
end

struct MixedDofHandler{dim,C,T} <: JuAFEM.AbstractDofHandler
    fieldhandlers::Vector{FieldHandler}
    cell_dofs::Vector{Int}
    cell_dofs_offset::Vector{Int}
    closed::ScalarWrapper{Bool}
    grid::MixedGrid{dim,C,T}
end

function MixedDofHandler(grid::MixedGrid)
    MixedDofHandler(FieldHandler[], Int[], Int[], JuAFEM.ScalarWrapper(false), grid)
end

getfieldnames(fh::FieldHandler) = [field.name for field in fh.fields]
getfielddims(fh::FieldHandler) = [field.dim for field in fh.fields]
getfieldinterpolations(fh::FieldHandler) = [field.interpolation for field in fh.fields]

"""
Returns the union of all the fields. Can be used as an iterable over all the fields in the problem.
"""
function getfieldnames(dh::MixedDofHandler)
    fieldnames = Vector{Symbol}()
    for fh in dh.fieldhandlers
        for name in getfieldnames(fh)
            push!(fieldnames, name)
        end
    end
    return unique!(fieldnames)
end

function getfielddim(dh::MixedDofHandler, name::Symbol)

    for fh in dh.fieldhandlers
        field_pos = findfirst(i->i == name, getfieldnames(fh))
        if field_pos > 0
            return fh.fields[field_pos].dim
        end
    end
    error("did not find field $name")
end

nfields(dh::MixedDofHandler) = length(getfieldnames(dh))

function Base.push!(dh::MixedDofHandler, fh::FieldHandler)
    @assert !isclosed(dh)
    push!(dh.fieldhandlers, fh)
    return dh
end

#import JuAFEM.close!
function close!(dh::MixedDofHandler{dim}) where {dim}

    @assert !JuAFEM.isclosed(dh)
    field_names = JuAFEM.getfieldnames(dh)  # all the fields in the problem
    numfields =  length(field_names)

    # Create dicts that stores created dofs
    # Each key should uniquely identify the given type
    vertexdicts = [Dict{Int, Array{Int}}() for _ in 1:numfields]
    edgedicts = [Dict{Tuple{Int,Int}, Array{Int}}() for _ in 1:numfields]
    facedicts = [Dict{Tuple{Int,Int}, Array{Int}}() for _ in 1:numfields]
    celldicts = [Dict{Int, Array{Int}}() for _ in 1:numfields]

    # Set initial values
    nextdof = 1  # next free dof to distribute
    push!(dh.cell_dofs_offset, nextdof)

    @debug "\n\nCreating dofs\n"
    for fh in dh.fieldhandlers
        # sort the cellset since we want to loop through the cells in a fixed order
        cellnumbers = sort(collect(fh.cellset))
        nextdof = _close!(
            dh,
            cellnumbers,
            JuAFEM.getfieldnames(fh),
            JuAFEM.getfielddims(fh),
            JuAFEM.getfieldinterpolations(fh),
            nextdof,
            vertexdicts,
            edgedicts,
            facedicts,
            celldicts)
    end
    dh.closed[] = true
    return dh
end

function _close!(dh::MixedDofHandler{dim}, cellnumbers, field_names, field_dims, field_interpolations, nextdof, vertexdicts, edgedicts, facedicts, celldicts) where {dim}

    ip_infos = JuAFEM.InterpolationInfo[]
    for interpolation in field_interpolations
        ip_info = JuAFEM.InterpolationInfo(interpolation)
        # these are implemented yet (or have not been tested)
        @assert(ip_info.nvertexdofs <= 1)
        @assert(ip_info.nedgedofs <= 1)
        @assert(ip_info.nfacedofs <= 1)
        @assert(ip_info.ncelldofs <= 1)  # not tested but probably works
        push!(ip_infos, ip_info)
    end

    # loop over all the cells, and distribute dofs for all the fields
    for ci in cellnumbers
        cell = dh.grid.cells[ci]
        cell_dofs = Int[]  # list of global dofs for each cell
        @debug "Creating dofs for cell #$ci"

        for fi in 1:length(field_names)
            @debug "\tfield: $(field_names[fi])"
            ip_info = ip_infos[fi]

            if ip_info.nvertexdofs > 0
                nextdof = add_vertex_dofs(cell_dofs, cell, vertexdicts[fi], field_dims[fi], ip_info.nvertexdofs, nextdof)
            end

            if ip_info.nfacedofs > 0
                nextdof = add_face_dofs(cell_dofs, cell, facedicts[fi], field_dims[fi], ip_info.nfacedofs, nextdof)
            end

            if ip_info.nedgedofs > 0
                nextdof = add_edge_dofs(cell_dofs, cell, edgedicts[fi], field_dims[fi], ip_info.nedgedofs, nextdof)
            end

            if ip_info.ncelldofs > 0
                nextdof = add_cell_dofs(cell_dofs, ci, celldicts[fi], field_dims[fi], ip_info.ncelldofs, nextdof)
            end

        end
        # after done creating dofs for the cell, push them to the global list
        push!(dh.cell_dofs, cell_dofs...)
        # push! the first index of the next cell to the offset vector
        push!(dh.cell_dofs_offset, length(dh.cell_dofs)+1)

        @debug "Dofs for cell #$ci:\n\t$cell_dofs"
    end # cell loop
    return nextdof
end

"""
Returns the next global dof number and an array of dofs.
If dofs have already been created for the object (vertex, face) then simply return those, otherwise create new dofs.
"""
function get_or_create_dofs!(nextdof, field_dim; dict, key)

    token = Base.ht_keyindex2!(dict, key)
    if token > 0  # vertex, face etc. visited before
        # reuse stored dofs (TODO unless field is discontinuous)
        @debug "\t\tkey: $key dofs: $(dict[key])  (reused dofs)"
        return nextdof, dict[key]

    else  # create new dofs
        dofs = collect(nextdof : nextdof + field_dim-1)
        @debug "\t\tkey: $key dofs: $dofs"
        Base._setindex!(dict, dofs, key, -token) #
        nextdof += field_dim
        return nextdof, dofs
    end
end

function add_vertex_dofs(cell_dofs, cell, vertexdict, field_dim, nvertexdofs, nextdof)
    for vertex in JuAFEM.vertices(cell)
        @debug "\tvertex #$vertex"
        nextdof, dofs = get_or_create_dofs!(nextdof, field_dim, dict=vertexdict, key=vertex)
        push!(cell_dofs, dofs...)
    end
    return nextdof
end

function add_face_dofs(cell_dofs, cell, facedict, field_dim, nfacedofs, nextdof)
    @assert nfacedofs == 1 "Currently only supports interpolations with nfacedofs = 1"

    for face in JuAFEM.faces(cell)
        sface = JuAFEM.sortface(face)
        @debug "\tface #$sface"
        nextdof, dofs = get_or_create_dofs!(nextdof, field_dim, dict=facedict, key=sface)
        push!(cell_dofs, dofs...)
    end
    return nextdof
end

function add_edge_dofs(cell_dofs, cell, edgedict, field_dim, nedgedofs, nextdof)
    @assert nedgedofs == 1 "Currently only supports interpolations with nedgedofs = 1"
    for edge in JuAFEM.edges(cell)
        sedge, dir = JuAFEM.sortedge(edge)
        @debug "\tedge #$sedge"
        nextdof, dofs = get_or_create_dofs!(nextdof, field_dim, dict=edgedict, key=sedge)
        push!(cell_dofs, dofs...)
    end
    return nextdof
end

function add_cell_dofs(cell_dofs, cell, celldict, field_dim, ncelldofs, nextdof)
    for celldof in 1:ncelldofs
        @debug "\tcell #$cell"
        nextdof, dofs = get_or_create_dofs!(nextdof, field_dim, dict=celldict, key=cell)
        push!(cell_dofs, dofs...)
    end
    return nextdof
end

# TODO if not too slow it can replace the "Grid-version"
function _create_sparsity_pattern(dh::MixedDofHandler, sym::Bool)
    ncells = getncells(dh.grid)
    N::Int = 0
    for element_id = 1:ncells  # TODO check for correctness
        n = ndofs_per_cell(dh, element_id)
        N += sym ? div(n*(n+1), 2) : n^2
    end
    N += ndofs(dh) # always add the diagonal elements
    I = Int[]; resize!(I, N)
    J = Int[]; resize!(J, N)

    cnt = 0
    for element_id in 1:ncells
        n = ndofs_per_cell(dh, element_id)
        global_dofs = zeros(Int, n)
        celldofs!(global_dofs, dh, element_id)
        @inbounds for j in 1:n, i in 1:n
            dofi = global_dofs[i]
            dofj = global_dofs[j]
            sym && (dofi > dofj && continue)
            cnt += 1
            if cnt > length(J)
                resize!(I, trunc(Int, length(I) * 1.5))
                resize!(J, trunc(Int, length(J) * 1.5))
            end
            I[cnt] = dofi
            J[cnt] = dofj
        end
    end
    @inbounds for d in 1:ndofs(dh)
        cnt += 1
        if cnt > length(J)
            resize!(I, trunc(Int, length(I) + ndofs(dh)))
            resize!(J, trunc(Int, length(J) + ndofs(dh)))
        end
        I[cnt] = d
        J[cnt] = d
    end
    resize!(I, cnt)
    resize!(J, cnt)
    V = zeros(length(I))
    K = sparse(I, J, V)
    return K

end

@inline create_sparsity_pattern(dh::MixedDofHandler) = _create_sparsity_pattern(dh, false)

function find_field(fh::FieldHandler, field_name::Symbol)
    j = findfirst(i->i == field_name, getfieldnames(fh))
    j == 0 && error("did not find field $field_name")
    return j
end

# Calculate the offset to the first local dof of a field
function field_offset(fh::FieldHandler, field_name::Symbol)
    offset = 0
    for i in 1:find_field(fh, field_name)-1
        offset += getnbasefunctions(getfieldinterpolations(fh)[i])::Int * getfielddims(fh)[i]
    end
    return offset
end


text:
@testset "Functions" begin
    Base.info("Testing functions...")
    inputvec = [-5.0, 0.0, 5.0]
    @test sigmoid(inputvec[2]) == 0.5
    @test sigmoid.(inputvec) == map(x -> 1.0 / (1.0 + exp(-x)), inputvec)
    @test relu(inputvec[1]) == zero(inputvec[1])
    @test relu.(inputvec) == [0.0, 0.0, inputvec[3]]

    inputvec_softmax = [0.0, 0.5, 1.0]
    @test softmax(inputvec_softmax) == [0.36787944117144233, 0.6065306597126334, 1.0] ./ 1.9744101008840758
    @test onehot([1,3,2], 1:3) == [[1,0,0] [0,0,1] [0,1,0]]
    @test crossentropyerror([1,0,0], [1,0,0]) == -log(1.0 + 1.0e-7)
    @test crossentropyerror([0.6,0.4,0], [0,0,1]) == -log(1.0e-7)
    Base.info("Done testing functions.")
end


text:
withenv("FORCE_COLOR" => true) do
    include("test_custompass.jl")
    include("test_bracketmatch.jl")
    include("test_highlighter.jl")
    include("test_rainbowbrackets.jl")
end


text:
module VortexHelperLilly
using ..Ahorn, Maple
using Cairo

@mapdef Entity "VortexHelper/Lilly" Lilly(
                                    x::Integer, y::Integer,
                                    width::Integer = 24, height::Integer = 24,
                                    maxLength::Int = 64)

const placements = Ahorn.PlacementDict(
    "Lilly (Vortex Helper)" => Ahorn.EntityPlacement(
        Lilly,
        "rectangle",
        Dict{String, Any}()
    )
)

Ahorn.resizable(entity::Lilly) = false, true
Ahorn.minimumSize(entity::Lilly) = 24, 24
Ahorn.nodeLimits(entity::Lilly) = 2, 2

const block = "objects/VortexHelper/squareBumperNew/block00"
const lillyFace = "objects/VortexHelper/squareBumperNew/face12" # chose specific frame
const armend = "objects/VortexHelper/squareBumperNew/armend"
const arm = "objects/VortexHelper/squareBumperNew/arm00"
const faceColor = (0, 208, 255, 255) ./ 255

function Ahorn.selection(entity::Lilly)
    x, y = Ahorn.position(entity)
    height = Int(get(entity.data, "height", 24))
    return Ahorn.Rectangle(x, y, 24, height)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Lilly, room::Maple.Room)
    maxLength = Int(get(entity.data, "maxLength", 64))
    height = Int(get(entity.data, "height", 24))

    save(ctx)
	set_antialias(ctx, 1)
	set_line_width(ctx, 1)
    set_dash(ctx, [0.6, 0.2])
    Ahorn.drawRectangle(ctx, -maxLength + 1, 1, 2 * maxLength + 23, height - 7, (0.0, 0.0, 0.0, 0.0), (1.0, 1.0, 1.0, 0.3))
    restore(ctx)

    # arms
    x = maxLength + 16
    while x > 8
        Ahorn.drawImage(ctx, arm, x, 0)
        Ahorn.drawImage(ctx, arm, 16-x, 0)
        x -= 8
    end

    tilesHeight = tilesHeight = div(height, 8)
    for i in 0 : tilesHeight - 1
        ty = (i == 0 ? 0 : (i == tilesHeight - 1 ? 16 : (i == tilesHeight - 2 ? 8 : 24)))
        # block
        Ahorn.drawImage(ctx, block, 0, i * 8, 0, ty, 8, 8)
        Ahorn.drawImage(ctx, block, 8, i * 8, 8, ty, 8, 8)
        Ahorn.drawImage(ctx, block, 16, i * 8, 16, ty, 8, 8)

        # arm ends
        Ahorn.drawImage(ctx, armend, maxLength + 16, i * 8, 0, ty, 8, 8)
        Ahorn.drawImage(ctx, armend, maxLength + 16, i * 8, 0, ty, 8, 8)
        Ahorn.drawImage(ctx, armend, maxLength + 16, i * 8, 0, ty, 8, 8)
        save(ctx)
        scale(ctx, -1, 1)
        Ahorn.drawImage(ctx, armend, maxLength - 8, i * 8, 0, ty, 8, 8)
        Ahorn.drawImage(ctx, armend, maxLength - 8, i * 8, 0, ty, 8, 8)
        Ahorn.drawImage(ctx, armend, maxLength - 8, i * 8, 0, ty, 8, 8)
        restore(ctx)
    end

    Ahorn.drawImage(ctx, lillyFace, 4, 5 + (height - 24) / 2, tint = faceColor)
end

end

max_stars_repo_name: catapillie/VortexHelper

max_issues_repo_name: catapillie/VortexHelper

max_forks_repo_name: catapillie/VortexHelper