text:
#####
##### Adapting structures to be able to pass them to GPU CUDA kernels.
#####

# Adapt an offset CuArray to work nicely with CUDA kernels.
Adapt.adapt_structure(to, x::OffsetArray) = OffsetArray(adapt(to, parent(x)), x.offsets)

# Need to adapt SubArray indices as well.
# See: https://github.com/JuliaGPU/Adapt.jl/issues/16
#Adapt.adapt_structure(to, A::SubArray{<:Any,<:Any,AT}) where {AT} =
#    SubArray(adapt(to, parent(A)), adapt.(Ref(to), parentindices(A)))

#####
##### Utilities that make it easier to juggle around Arrays and CuArrays.
#####

         array_type(::CPU) = Array
@hascuda array_type(::GPU) = CuArray

#####
##### Convenient macros
#####

macro loop_xyz(i, j, k, grid, expr)
    return esc(
        quote
            @loop for $k in (1:$grid.Nz; (blockIdx().z - 1) * blockDim().z + threadIdx().z)
                @loop for $j in (1:$grid.Ny; (blockIdx().y - 1) * blockDim().y + threadIdx().y)
                    @loop for $i in (1:$grid.Nx; (blockIdx().x - 1) * blockDim().x + threadIdx().x)
                        $expr
                    end
                end
            end
        end)
end

macro loop_xy(i, j, grid, expr)
    return esc(
        quote
            @loop for $j in (1:$grid.Ny; (blockIdx().y - 1) * blockDim().y + threadIdx().y)
                @loop for $i in (1:$grid.Nx; (blockIdx().x - 1) * blockDim().x + threadIdx().x)
                    $expr
                end
            end
        end)
end

macro loop_xz(i, k, grid, expr)
    return esc(
        quote
            @loop for $k in (1:$grid.Nz; (blockIdx().y - 1) * blockDim().y + threadIdx().y)
                @loop for $i in (1:$grid.Nx; (blockIdx().x - 1) * blockDim().x + threadIdx().x)
                    $expr
                end
            end
        end)
end

#####
##### Convenient definitions
#####

"""
    second

A `Float64` constant equal to 1.0. Useful for increasing the clarity of scripts, e.g. `Δt = 1second`.
"""
const second = 1.0

"""
    minute

A `Float64` constant equal to 60`second`. Useful for increasing the clarity of scripts, e.g. `Δt = 15minute`.
"""
const minute = 60second

"""
    hour

A `Float64` constant equal to 60`minute`. Useful for increasing the clarity of scripts, e.g. `Δt = 3hour`.
"""
const hour   = 60minute

"""
    day

A `Float64` constant equal to 24`hour`. Useful for increasing the clarity of scripts, e.g. `Δt = 0.5day`.
"""
const day    = 24hour

"""
    meter

A `Float64` constant equal to 1.0. Useful for increasing the clarity of scripts, e.g. `Lx = 100meter`.
"""
const meter = 1.0

"""
    kilometer

A `Float64` constant equal to 1000`meter`. Useful for increasing the clarity of scripts, e.g. `Lx = 250kilometer`.
"""
const kilometer = 1000meter

"""
    KiB

A `Float64` constant equal to 1024.0. Useful for increasing the clarity of scripts, e.g. `max_filesize = 250KiB`.
"""
const KiB = 1024.0

"""
    MiB

A `Float64` constant equal to 1024`KiB`. Useful for increasing the clarity of scripts, e.g. `max_filesize = 100MiB`.
"""
const MiB = 1024KiB

"""
    GiB

A `Float64` constant equal to 1024`MiB`. Useful for increasing the clarity of scripts, e.g. `max_filesize = 50GiB`.
"""
const GiB = 1024MiB

"""
    TiB

A `Float64` constant equal to 1024`GiB`. Useful for increasing the clarity of scripts, e.g. `max_filesize = 2TiB`.
"""
const TiB = 1024GiB

#####
##### Pretty printing
#####

"""
    prettytime(t)

Convert a floating point value `t` representing an amount of time in seconds to a more
human-friendly formatted string with three decimal places. Depending on the value of `t`
the string will be formatted to show `t` in nanoseconds (ns), microseconds (μs),
milliseconds (ms), seconds (s), minutes (min), hours (hr), or days (day).
"""
function prettytime(t)
    # Modified from: https://github.com/JuliaCI/BenchmarkTools.jl/blob/master/src/trials.jl
    iszero(t) && return "0.000 s"

    if t < 1e-6
        value, units = t * 1e9, "ns"
    elseif t < 1e-3
        value, units = t * 1e6, "μs"
    elseif t < 1
        value, units = t * 1e3, "ms"
    elseif t < minute
        value, units = t, "s"
    elseif t < hour
        value, units = t / minute, "min"
    elseif t < day
        value, units = t / hour, "hr"
    else
        value, units = t / day, "day"
    end

    return @sprintf("%.3f", value) * " " * units
end

"""
    pretty_filesize(s, suffix="B")

Convert a floating point value `s` representing a file size to a more human-friendly
formatted string with one decimal places with a `suffix` defaulting to "B". Depending on
the value of `s` the string will be formatted to show `s` using an SI prefix from bytes,
kiB (1024 bytes), MiB (1024² bytes), and so on up to YiB (1024⁸ bytes).
"""
function pretty_filesize(s, suffix="B")
    # Modified from: https://stackoverflow.com/a/1094933
    for unit in ["", "Ki", "Mi", "Gi", "Ti", "Pi", "Ei", "Zi"]
        abs(s) < 1024 && return @sprintf("%3.1f %s%s", s, unit, suffix)
        s /= 1024
    end
    return @sprintf("%.1f %s%s", s, "Yi", suffix)
end

#####
##### Creating fields by dispatching on architecture
#####

function OffsetArray(underlying_data, grid::AbstractGrid)
    # Starting and ending indices for the offset array.
    i1, i2 = 1 - grid.Hx, grid.Nx + grid.Hx
    j1, j2 = 1 - grid.Hy, grid.Ny + grid.Hy
    k1, k2 = 1 - grid.Hz, grid.Nz + grid.Hz

    return OffsetArray(underlying_data, i1:i2, j1:j2, k1:k2)
end

function Base.zeros(T, ::CPU, grid)
    underlying_data = zeros(T, grid.Tx, grid.Ty, grid.Tz)
    return OffsetArray(underlying_data, grid)
end

function Base.zeros(T, ::GPU, grid)
    underlying_data = CuArray{T}(undef, grid.Tx, grid.Ty, grid.Tz)
    underlying_data .= 0  # Gotta do this otherwise you might end up with a few NaN values!
    return OffsetArray(underlying_data, grid)
end

Base.zeros(T, ::CPU, grid, Nx, Ny, Nz) = zeros(T, Nx, Ny, Nz)
Base.zeros(T, ::GPU, grid, Nx, Ny, Nz) = zeros(T, Nx, Ny, Nz) |> CuArray

# Default to type of Grid
Base.zeros(arch, grid::AbstractGrid{T}) where T = zeros(T, arch, grid)
Base.zeros(arch, grid::AbstractGrid{T}, Nx, Ny, Nz) where T = zeros(T, arch, grid, Nx, Ny, Nz)

#####
##### Courant–Friedrichs–Lewy (CFL) condition number calculation
#####

# Note: these functions will have to be refactored to work on non-uniform grids.

"Returns the time-scale for advection on a regular grid across a single grid cell."
function cell_advection_timescale(u, v, w, grid)
    umax = maximum(abs, u)
    vmax = maximum(abs, v)
    wmax = maximum(abs, w)

    Δx = grid.Δx
    Δy = grid.Δy
    Δz = grid.Δz

    return min(Δx/umax, Δy/vmax, Δz/wmax)
end

cell_advection_timescale(model) =
    cell_advection_timescale(model.velocities.u.data.parent,
                             model.velocities.v.data.parent,
                             model.velocities.w.data.parent,
                             model.grid)

#####
##### Adaptive time stepping
#####

"""
    TimeStepWizard{T}

    TimeStepWizard(cfl=0.1, max_change=2.0, min_change=0.5, max_Δt=Inf)

A type for calculating adaptive time steps based on capping the CFL number at `cfl`.

On calling `update_Δt!(wizard, model)`, the `TimeStepWizard` computes a time-step such that
``cfl = max(u/Δx, v/Δy, w/Δz) Δt``, where ``max(u/Δx, v/Δy, w/Δz)`` is the maximum ratio
between model velocity and along-velocity grid spacing anywhere on the model grid. The new
`Δt` is constrained to change by a multiplicative factor no more than `max_change` or no
less than `min_change` from the previous `Δt`, and to be no greater in absolute magnitude
than `max_Δt`.
"""
Base.@kwdef mutable struct TimeStepWizard{T}
              cfl :: T = 0.1
    cfl_diffusion :: T = 2e-2
       max_change :: T = 2.0
       min_change :: T = 0.5
           max_Δt :: T = Inf
               Δt :: T = 0.01
end

"""
    update_Δt!(wizard, model)

Compute `wizard.Δt` given the velocities and diffusivities of `model`, and the parameters
of `wizard`.
"""
function update_Δt!(wizard, model)
    Δt = wizard.cfl * cell_advection_timescale(model)

    # Put the kibosh on if needed
    Δt = min(wizard.max_change * wizard.Δt, Δt)
    Δt = max(wizard.min_change * wizard.Δt, Δt)
    Δt = min(wizard.max_Δt, Δt)

    wizard.Δt = Δt

    return nothing
end

#####
##### Some utilities for tupling
#####

tupleit(t::Tuple) = t
tupleit(a::AbstractArray) = Tuple(a)
tupleit(nt) = tuple(nt)

parenttuple(obj) = Tuple(f.data.parent for f in obj)

@inline datatuple(obj::Nothing) = nothing
@inline datatuple(obj::AbstractArray) = obj
@inline datatuple(obj::Tuple) = Tuple(datatuple(o) for o in obj)
@inline datatuple(obj::NamedTuple) = NamedTuple{propertynames(obj)}(datatuple(o) for o in obj)
@inline datatuples(objs...) = (datatuple(obj) for obj in objs)

#####
##### Dynamic launch configuration
#####

function launch_config(grid, dims)
    return function (kernel)
        fun = kernel.fun
        config = launch_configuration(fun)

        # Adapt the suggested config from 1D to the requested grid dimensions
        if dims == :xyz
            t = floor(Int, cbrt(config.threads))
            threads = [t, t, t]
            blocks  = ceil.(Int, [grid.Nx, grid.Ny, grid.Nz] ./ t)
        elseif dims == :xy
            t = floor(Int, sqrt(config.threads))
            threads = [t, t]
            blocks  = ceil.(Int, [grid.Nx, grid.Ny] ./ t)
        elseif dims == :xz
            t = floor(Int, sqrt(config.threads))
            threads = [t, t]
            blocks  = ceil.(Int, [grid.Nx, grid.Nz] ./ t)
        else
            error("Unsupported launch configuration: $dims")
        end

        return (threads=Tuple(threads), blocks=Tuple(blocks))
    end
end

#####
##### Utilities shared between diagnostics and output writers
#####

defaultname(::AbstractDiagnostic, nelems) = Symbol(:diag, nelems+1)
defaultname(::AbstractOutputWriter, nelems) = Symbol(:writer, nelems+1)

const DiagOrWriterDict = OrderedDict{S, <:Union{AbstractDiagnostic, AbstractOutputWriter}} where S

function push!(container::DiagOrWriterDict, elem)
    name = defaultname(elem, length(container))
    container[name] = elem
    return nothing
end

getindex(container::DiagOrWriterDict, inds::Integer...) = getindex(container.vals, inds...)
setindex!(container::DiagOrWriterDict, newvals, inds::Integer...) = setindex!(container.vals, newvals, inds...)

function push!(container::DiagOrWriterDict, elems...)
    for elem in elems
        push!(container, elem)
    end
    return nothing
end

"""
    validate_interval(frequency, interval)

Ensure that frequency and interval are not both `nothing`.
"""
function validate_interval(frequency, interval)
    isnothing(frequency) && isnothing(interval) && @error "Must specify a frequency or interval!"
    return nothing
end

has_interval(obj) = :interval in propertynames(obj) && obj.interval != nothing
has_frequency(obj) = :frequency in propertynames(obj) && obj.frequency != nothing

function interval_is_ripe(clock, obj)
    if has_interval(obj) && clock.time >= obj.previous + obj.interval
        obj.previous = clock.time - rem(clock.time, obj.interval)
        return true
    else
        return false
    end
end

frequency_is_ripe(clock, obj) = has_frequency(obj) && clock.iteration % obj.frequency == 0

function time_to_run(clock, output_writer)

    interval_is_ripe(clock, output_writer) && return true
    frequency_is_ripe(clock, output_writer) && return true

    # If the output writer does not specify an interval or frequency,
    # it is unable to write output and we throw an error as a convenience.
    has_interval(output_writer) || has_frequency(output_writer) ||
        error("$(typeof(output_writer)) must have a frequency or interval specified!")

    return false
end

#####
##### Utils for models
#####

function ordered_dict_show(dict, padchar)
    if length(dict) == 0
        return string(typeof(dict), " with no entries")
    elseif length(dict) == 1
        return string(typeof(dict), " with 1 entry:", '\n',
                      padchar, "   └── ", dict.keys[1], " => ", typeof(dict.vals[1]))
    else
        return string(typeof(dict), " with $(length(dict)) entries:", '\n',
                      Tuple(string(padchar,
                                   "   ├── ", name, " => ", typeof(dict[name]), '\n')
                            for name in dict.keys[1:end-1]
                           )...,
                           padchar, "   └── ", dict.keys[end], " => ", typeof(dict.vals[end])
                      )
    end
end

"""
    with_tracers(tracers, initial_tuple, tracer_default)

Create a tuple corresponding to the solution variables `u`, `v`, `w`,
and `tracers`. `initial_tuple` is a `NamedTuple` that at least has
fields `u`, `v`, and `w`, and may have some fields corresponding to
the names in `tracers`. `tracer_default` is a function that produces
a default tuple value for each tracer if not included in `initial_tuple`.
"""
function with_tracers(tracers, initial_tuple::NamedTuple, tracer_default; with_velocities=false)
    solution_values = [] # Array{Any, 1}
    solution_names = []

    if with_velocities
        push!(solution_values, initial_tuple.u)
        push!(solution_values, initial_tuple.v)
        push!(solution_values, initial_tuple.w)

        append!(solution_names, [:u, :v, :w])
    end

    for name in tracers
        tracer_elem = name ∈ propertynames(initial_tuple) ?
                        getproperty(initial_tuple, name) :
                        tracer_default(tracers, initial_tuple)

        push!(solution_values, tracer_elem)
    end

    append!(solution_names, tracers)

    return NamedTuple{Tuple(solution_names)}(Tuple(solution_values))
end


text:
#production and acquisition

function ICCost!(options::zefiroOptions,sizes::zefiroSizes,data::zefiroData,output::zefiroOutput)
    key = (options.postprocess ? "postprocessing" : "evaluating")
    println(key*" port cost")
    port!(options,sizes,data,output)
    println(key*" installationOfTheComponents cost")
    installationOfTheComponents!(options,sizes,data,output)
    println(key*" commissioning cost")
    commissioning!(options,sizes,data,output)
    println(key*" insurance cost")
    insurance!(options,sizes,data,output)

    nothing
end

function port!(options::zefiroOptions,sizes::zefiroSizes,data::zefiroData,output::zefiroOutput)

    if options.port == 0
        if options.postprocess
            nothing
        else
            for i in 1:sizes.winds,j in 1:sizes.turbines
                output.CAPEX[i,j] += (data.Cport[i] + data.N1d[j]*data.Lr[i])*data.turbineQuantity[j]
                output.port[i,j] += (data.Cport[i] + data.N1d[j]*data.Lr[i])*data.turbineQuantity[j]
            end
        end
    else
        error("port mode not found.")
    end
    nothing

end

function installationOfTheComponents!(options::zefiroOptions,sizes::zefiroSizes,data::zefiroData,output::zefiroOutput)

    if options.installationOfTheComponents == 0
        if options.postprocess
            nothing
        else
            output.CAPEX .+= transpose(data.Ccomp).*transpose(data.turbineQuantity)
            output.installationOfTheComponents .+= transpose(data.Ccomp).*transpose(data.turbineQuantity)
        end
    else
        error("installationOfTheComponents mode not found.")
    end
    nothing

end

function commissioning!(options::zefiroOptions,sizes::zefiroSizes,data::zefiroData,output::zefiroOutput)

    if options.commissioning == 0
        if options.postprocess
            nothing
        else
            output.CAPEX .+= transpose(data.Lcom.*data.Ncom).*transpose(data.turbineQuantity)
            output.commissioning .+= transpose(data.Lcom.*data.Ncom).*transpose(data.turbineQuantity)
        end
    else
        error("commissioning mode not found.")
    end
    nothing

end

function insurance!(options::zefiroOptions,sizes::zefiroSizes,data::zefiroData,output::zefiroOutput)

    if options.insurance == 0
        if options.postprocess
            nothing
        else
            for i in 1:sizes.winds,j in 1:sizes.turbines
                output.CAPEX[i,j] += data.Cins[i]*data.IC[j]
                output.insurance[i,j] += data.Cins[i]*data.IC[j]
            end
        end
    else
        error("insurance mode not found.")
    end
    nothing

end

ext: jl

max_stars_repo_name: Thiago-NovaesB/zefiro

max_issues_repo_name: Thiago-NovaesB/zefiro

max_forks_repo_name: Thiago-NovaesB/zefiro


text:
function DistanceResto(id_aliment,id_work)

    simple_graph = loadgraph("SaintVicNoChildrenLeftBehind.graphml", GraphIO.GraphML.GraphMLFormat()) ;

    chemin = zeros(Int32,(length(id_work),length(id_aliment))) ;

    @time @threads for i = 1:length(id_work)
        for j=1:length(id_aliment)
            ch = a_star(simple_graph,id_work[i],id_aliment[j]) ;
            chemin[i,j] = length(ch) ;
        end

        println("J\'ai fini la",i,"eme entreprise")

    end

    return(chemin)

end

DistanceResto


text:
module Domains

export Tdjsets, find!, union!, setallroot!, Tdomains

mutable struct Tdjsets
    sets::Vector{UInt32}
    setsz::Vector{UInt32}
end

# constructor
function Tdjsets(N)
    sets = Array{UInt32,1}(0x00000001:UInt32(N))
    setsz = ones(UInt32, size(sets))
    Tdjsets(sets, setsz)
end

function find!( djsets::Tdjsets, vid::UInt32 )
    # find root id or domain id
    rid = vid
    while rid != djsets.sets[rid]
        rid = djsets.sets[rid]
    end

    # path compression
    # current id
    cid = vid
    while rid != cid
        # parent id
        pid = djsets.sets[cid]
        djsets.sets[cid] = rid
        cid = pid
    end
    return rid
end

# import Base.union! for extention
import Base.union!
function union!( djsets::Tdjsets, sid1::UInt32, sid2::UInt32 )
    # find root id
    rid1 = find!(djsets, sid1)
    rid2 = find!(djsets, sid2)

    if rid1 == rid2
        # already in the same domain
        return
    end

    # reduce set number
    if djsets.setsz[ rid1 ] >= djsets.setsz[ rid2 ]
        # assign sid1 as the parent of sid2
        djsets.sets[ rid2 ] = rid1
        djsets.setsz[ rid1 ] += djsets.setsz[ rid2 ]
    else
        djsets.sets[ rid1 ] = rid2
        djsets.setsz[ rid2 ] += djsets.setsz[ rid1 ]
    end
end

function setallroot!( djsets::Tdjsets )
    # label all the voxels to root id
    for vid in 0x00000001:UInt32( length( djsets.sets ) )
        # with patch compress
        # all the voxels will be labeled as root id
        rid = find!(djsets, vid)
        djsets.setsz[vid] = djsets.setsz[rid]
    end
    return djsets.sets
end

# size of each label in a domain
# key is the label id in ground truth
# value is the number of voxels in that label
const Tdlsz = Dict{UInt32,UInt32}

"""
compare two domain label sizes to get the number of same voxel pair and different voxel pair
In default, it is foreground restricted.
"""
function get_pair_num(dlsz1::Tdlsz, dlsz2::Tdlsz, is_fr::Bool = true)
    n_same_pair = UInt32(0)
    n_diff_pair = UInt32(0)
    for (lid1, sz1) in dlsz1
        for (lid2, sz2) in dlsz2
            if is_fr && lid2==UInt32(0)
                continue
            end
            if lid1 == lid2
                # have common segment id, merge together
                n_same_pair += sz1 * sz2
            else
                # do not have common id, create new one
                n_diff_pair += sz1 * sz2
            end
        end
    end
    return n_same_pair, n_diff_pair
end

# union domain label size dict 2 to 1, only 1 was changed
function union!( dlsz1::Tdlsz, dlsz2::Tdlsz )
    for (lid1, sz1) in dlsz1
        for (lid2, sz2) in dlsz2
            if lid1 == lid2
                # have common segment id, merge together
                dlsz1[lid1] += sz2
            else
                # do not have common id, create new one
                dlsz1[lid2] = sz2
            end
        end
    end
end

# list of dictionary, each represents the label sizes
const Tdlszes = Array{Tdlsz,1}

mutable struct Tdomains
    # domain label sizes
    dlszes::Tdlszes
    # disjoint sets
    djsets::Tdjsets
end

# constructor function
function Tdomains(seg::Array{T,3}) where T
    # number of voxels
    N = T( length(seg) )
    # initialize the disjoint sets
    djsets = Tdjsets( N )

    # initialize the dms as an empty vector/list/1D array
    dlszes = Tdlszes([])
    seg_flat = seg[:]
    for vid in T(1):N
        # initial manual labeled segment id
        push!(dlszes, Tdlsz( seg_flat[vid] => 0x00000001) )
    end
    Tdomains(dlszes, djsets)
end

# find the corresponding domain of a voxel
function find!(dms::Tdomains, vid::UInt32)
    rid = find!(dms.djsets, vid)
    dlsz = dms.dlszes[ rid ]
    return rid, dlsz
end

function union!(dms::Tdomains, rid1::UInt32, dlsz1::Tdlsz, rid2::UInt32, dlsz2::Tdlsz)
    # alread in one domain
    if rid1 == rid2
        return
    end

    # attach the small one to the big one to make the tree as flat as possible
    union!( dms.djsets, rid1, rid2 )

    # merge small one to big one
    if dms.djsets.setsz[ rid1 ] >= dms.djsets.setsz[ rid2 ]
        union!(dms.dlszes[rid1], dms.dlszes[rid2])
    else
        union!(dms.dlszes[rid2], dms.dlszes[rid1])
    end
end

# union the two domains of two voxel ids
function union!(dms::Tdomains, vid1::UInt32, vid2::UInt32)
    # domain id and domain
    rid1, dlsz1 = find!(dms, vid1)
    rid2, dlsz2 = find!(dms, vid2)

    union!(dms, rid1, dlsz1, rid2, dlsz2)
end

end # end of module


text:
using EffectiveWaves, Test

#tests the values of the integrand against results from the Mathematica file integrate_hankels.nb, which used no approximations.

# When these tests were done there was an error in the integration domain of B. Likely these will need to be redone.
@testset "Compare integrand values against brute force Mathematica result" begin

    # using OffsetArrays
    import EffectiveWaves: BS_matrices

    function intergrand_kernels(X, a12k; M = M, θin=θin, num_coefs = 2000)

        q = min(Int(floor(a12k/(X[2]-X[1]))),J)
        B, S = BS_matrices(X, a12k; θin=θin, M=M, num_coefs = num_coefs)

        function intergrand(l,j,m,n)
            P = S[j-l,n-m]
            Q = (abs(j-l) <= q) ? (B[j-l,n-m] - S[j-l,n-m]) : zero(Complex{Float64})
            P + Q
        end

        [intergrand(l,j,m,n) for l=0:J, m=-M:M, j=0:J, n=-M:M]
    end

    h=0.02; max_x=2.0;
    J = Int(round(max_x/h))
    x = (0:J).*h
    i1 = 1 + Int(round(1/h)) # for x1 = 1
    i2 = 1 + Int(round(2/h)) # for x1 = 2
    a12k = 1.0

    # # both commented tests below should also pass, but am commenting to save time
    # θin=0.; M=0; n=0;
    #
    # # math was produced by Mathematica
    # math00 = [172.31259187840652, 143.34097100449316, 172.31259187840536]
    # intergrand_quad = intergrand_kernels(x, a12k; θin=θin, M=M, num_coefs = 2000);
    # julia00 = [sum(abs.(intergrand_quad[i,M+1,:,M+1+n])) for i in [1,i1,i2]];
    # @test maximum(abs.(1.0 .- julia00./math00)) < 1e-6
    #
    # math0p50 = [192.61985415285892, 155.98244729650008, 192.61985415285886]
    # θin=0.5;
    # intergrand_quad = intergrand_kernels(x, a12k; θin=θin, M=M, num_coefs = 3000);
    # julia0p50 = [sum(abs.(intergrand_quad[i,M+1,:,M+1+n])) for i in [1,i1,i2]];
    # @test maximum(abs.(1.0 .- julia0p50./math0p50)) < 3e-6

    math0p32 = [241.5144625759003, 211.13367934660897, 181.54133990599098]
    θin=0.3; M=2; n=2;
    intergrand_quad = intergrand_kernels(x, a12k; θin=θin, M=M, num_coefs = 2000);
    julia0p32 = [sum(abs.(intergrand_quad[i,M+1,:,M+1+n])) for i in [1,i1,i2]];
    @test maximum(abs.(1.0 .- julia0p32./math0p32)) < 4e-6

end


text:
using Flux
using Flux: @functor

using LinearAlgebra: tril!

abstract type AbstractAttention end

using ZygoteRules: AContext
import ZygoteRules: _pullback

function create_atten_mask(T::Type, mask::AbstractArray, future::Bool=false)
  #size(mask) == (q, k, n, b)

  ql, kl = size(mask)
  mask = copy(mask)

  maskval = convert(T, -1e9)
  !future && batched_tril!(mask, -1)
  mask .= (1 .- mask) .* maskval
  return mask
end

Flux.@nograd create_atten_mask

struct MultiheadAttention <: AbstractAttention
    head::Int
    future::Bool
    iqproj::Dense
    ikproj::Dense
    ivproj::Dense
    oproj::Dense
    drop::Dropout
end

Flux.functor(mh::MultiheadAttention) = (mh.iqproj, mh.ikproj, mh.ivproj, mh.oproj), m -> MultiheadAttention(mh.head, mh.future, m..., mh.drop)

"""
    MultiheadAttention(head::Int, is::Int, hs::Int, os::Int;
                       future::Bool=true, pdrop = 0.1)

Multihead dot product Attention Layer, `head` is the number of head,
`is` is the input size, `hs` is the hidden size of input projection layer of each head,
`os` is the output size. When `future` is `false`, the k-th token can't see tokens at > k.
`pdrop` is the dropout rate.
"""
MultiheadAttention(head::Int,
                   is::Int,
                   hs::Int,
                   os::Int;
                   future::Bool=true, pdrop = 0.1) = MultiheadAttention(head,
                                                                        future,
                                                                        Dense(is, hs*head),
                                                                        Dense(is, hs*head),
                                                                        Dense(is, hs*head),
                                                                        Dense(hs*head, os),
                                                                        Dropout(pdrop),
                                                                        )

function Base.show(io::IO, mh::MultiheadAttention)
    hs = div(size(mh.iqproj.W)[1], mh.head)
    is = size(mh.iqproj.W)[end]
    os = size(mh.oproj.W)[1]

    print(io, "MultiheadAttention(")
    print(io, "head=$(mh.head), ")
    print(io, "head_size=$(hs), ")
    print(io, "$(is)=>$(os)")

    if Flux.istraining()
        print(io, ", dropout=$(mh.drop.p))")
    else
        print(io, ")")
    end
end

function (mh::MultiheadAttention)(query::A1,
                                  key::A2,
                                  value::A3;
                                  mask=nothing) where {T,
                                                       A1 <: Abstract3DTensor{T},
                                                       A2 <: Abstract3DTensor{T},
                                                       A3 <: Abstract3DTensor{T}}
    qs = size(query)
    ks = size(key)
    vs = size(value)

    #size(ipq) == (h, q_seq_len, batch)
    ipq = @toNd mh.iqproj(query)
    ipk = @toNd mh.ikproj(key)
    ipv = @toNd mh.ivproj(value)

    h = size(ipq, 1)
    hs = div(h, mh.head)

    #size(ipq) == (hs, q_seq_len, head, batch)
    ipq = permutedims(reshape(ipq, hs, mh.head, qs[2], qs[3]), [1, 3, 2, 4])
    ipk = permutedims(reshape(ipk, hs, mh.head, ks[2], ks[3]), [1, 3, 2, 4])
    ipv = permutedims(reshape(ipv, hs, mh.head, vs[2], vs[3]), [1, 3, 2, 4])

    #size(ipq) == (hs, q_seq_len, head * batch)
    ipq = reshape(ipq, hs, qs[2], :)
    ipk = reshape(ipk, hs, ks[2], :)
    ipv = reshape(ipv, hs, vs[2], :)

    atten = attention(ipq,ipk,ipv,
                      mask;
                      future=mh.future,
                      dropout=mh.drop)

    atten = permutedims(reshape(atten, hs, qs[2], mh.head, qs[3]), [1, 3, 2, 4]) #size(atten) == (hs, head, ql, b)
    atten = reshape(atten, h, qs[2], qs[3]) #size(atten) == (h, ql, b)

    out = @toNd mh.oproj(atten)
    out #size(out) == (h, q_seq_len, batch)
end

function (mh::MultiheadAttention)(query::A1,
                                  key::A2,
                                  value::A3;
                                  mask=nothing) where {T,
                                                           A1 <: AbstractMatrix{T},
                                                           A2 <: AbstractMatrix{T},
                                                           A3 <: AbstractMatrix{T}}

    # size(query) == (dims, seq_len)
    ipq = mh.iqproj(query)
    ipk = mh.ikproj(key)
    ipv = mh.ivproj(value)

    h = size(ipq)[1] #h == hs * head
    hs = div(h, mh.head)

    #size(hq) == (hs, seq_len, head)
    hq = permutedims(reshape(ipq, hs, mh.head, :), [1, 3, 2])
    hk = permutedims(reshape(ipk, hs, mh.head, :), [1, 3, 2])
    hv = permutedims(reshape(ipv, hs, mh.head, :), [1, 3, 2])

    atten = attention(hq, hk, hv,
                      mask;
                      future=mh.future,
                      dropout=mh.drop)

    # size(atten) == (head*hs, seq_len)
    atten = reshape(permutedims(atten, [1, 3, 2]), h, :)

    mh.oproj(atten)
end

# unused function
# only for understand how attention works
# function attention(query::AbstractMatrix{T},
#                    key::AbstractMatrix{T},
#                    value::AbstractMatrix{T};
#                    mask=nothing, future::Bool = false,
#                    dropout=nothing) where T
#     # size(query) == (dims, {q,k}_seq_len) == size(key) == size(value)
#     # size(score) == (k_seq_len, q_seq_len)
#     dk = size(key)[1]
#     score = transpose(key) * query
#     score = score ./ convert(T, sqrt(dk))

#     if mask !== nothing
#         @. mask = (1 - mask) * convert(T, -1e9)
#         score = score .+ mask
#     end

#     if !future
#         fmask = tril!(fill!(similar(score), convert(T, -1e9)), -1)
#         score = score .+ fmask
#     end

#     score = softmax(score)
#     dropout !== nothing && (score = dropout(score))
#     value * score #size(return) == (dims, q_seq_len)
# end

function attention(query::A1,
                           key::A2,
                           value::A3,
                           mask; future::Bool = false,
                           dropout=nothing) where {T,
                                                       A1 <: Abstract3DTensor{T},
                                                       A2 <: Abstract3DTensor{T},
                                                       A3 <: Abstract3DTensor{T}}
    #size(query) == (dims, {q,k}_seq_len, batch) == size(key) == size(value)
    #size(score) == (k_seq_len, q_seq_len, batch)
    dk = size(key, 1)
    score = batchedmul(key, query; transA = true)
    score = score ./ convert(T, sqrt(dk))

    s = size(score)
    if mask !== nothing
        ms = size(mask)

        if s[end] == ms[end]
            score = score .+ create_atten_mask(T, mask, future)
        else
            bxn = s[end]
            b = ms[end]
            score = reshape(reshape(score, s[1:end-1]..., :, b) .+ reshape(create_atten_mask(T, mask, future), ms[1:end-1]..., 1, b), s)
        end
    end

    score = @toNd softmax(score) #reshape(softmax(reshape(score, s[1], :)) , s)
    dropout !== nothing && (score = dropout(score))
    batchedmul(value, score) #size(return) == (dims, q_seq_len, batch)
end


text:
export PooledString
"""
    PooledString

A singleton type that can be used for signaling that a column of a csv file should be pooled,
with the output array type being a `PooledArray`.
"""
struct PooledString <: AbstractString end

# PointerString is an internal-only type for efficiently tracking string data + length
# all strings indexed from a column/row will always be a full String
# specifically, it allows avoiding materializing full Strings for pooled string columns while parsing
# and allows a fastpath for materializing a full String when no escaping is needed
struct PointerString <: AbstractString
    ptr::Ptr{UInt8}
    len::Int
end

function Base.hash(s::PointerString, h::UInt)
    h += Base.memhash_seed
    ccall(Base.memhash, UInt, (Ptr{UInt8}, Csize_t, UInt32), s.ptr, s.len, h % UInt32) + h
end

import Base: ==
function ==(x::String, y::PointerString)
    sizeof(x) == y.len && ccall(:memcmp, Cint, (Ptr{UInt8}, Ptr{UInt8}, Csize_t), pointer(x), y.ptr, y.len) == 0
end
function ==(x::PointerString, y::PointerString)
    x.len == y.len && ccall(:memcmp, Cint, (Ptr{UInt8}, Ptr{UInt8}, Csize_t), x.ptr, y.ptr, y.len) == 0
end
==(y::PointerString, x::String) = x == y

Base.ncodeunits(s::PointerString) = s.len
@inline function Base.codeunit(s::PointerString, i::Integer)
    @boundscheck checkbounds(s, i)
    GC.@preserve s unsafe_load(s.ptr + i - 1)
end

Base.String(x::PointerString) = _unsafe_string(x.ptr, x.len)

# column bit flags; useful so we don't have to pass a bunch of arguments/state around manually

# whether the user provided the type or not
const USER       = 0b00000001
user(flag) = flag & USER > 0

# whether any missing values have been found in this column so far
const ANYMISSING = 0b00000010
anymissing(flag) = flag & ANYMISSING > 0

# whether a column type has been detected yet
const TYPEDETECTED = 0b00000100
typedetected(flag) = flag & TYPEDETECTED > 0

# whether a column will be "dropped" from the select/drop keyword arguments
const WILLDROP = 0b00001000
willdrop(flag) = flag & WILLDROP > 0

# whether strings should be lazy; results in LazyStringVectors
# this setting isn't per column, but we store it on the column bit flags anyway for convenience
const LAZYSTRINGS = 0b00010000
lazystrings(flag) = flag & LAZYSTRINGS > 0

const MAYBEPOOLED = 0b00100000
maybepooled(flag) = flag & MAYBEPOOLED > 0
# ~95% z-score, 10% MoE
const POOLSAMPLESIZE = 100

flag(T, lazystrings) = (T === Union{} ? 0x00 : ((USER | TYPEDETECTED) | (hasmissingtype(T) ? ANYMISSING : 0x00))) | (lazystrings ? LAZYSTRINGS : 0x00)

# we define our own bit flag on a Parsers.ReturnCode to signal if a column needs to promote to string
const PROMOTE_TO_STRING = 0b0100000000000000 % Int16
promote_to_string(code) = code & PROMOTE_TO_STRING > 0

hasmissingtype(T) = T === Missing || T !== Core.Compiler.typesubtract(T, Missing)

@inline function promote_types(@nospecialize(T), @nospecialize(S))
    if T === Union{} || S === Union{} || T === Missing || S === Missing || T === S || nonmissingtype(T) === nonmissingtype(S)
        return Union{T, S}
    elseif T === Int64
        return S === Float64 ? S : S === Union{Float64, Missing} ? S : hasmissingtype(S) ? Union{String, Missing} : String
    elseif T === Union{Int64, Missing}
        return S === Float64 || S === Union{Float64, Missing} ? Union{Float64, Missing} : Union{String, Missing}
    elseif T === Float64
        return S === Int64 ? T : S === Union{Int64, Missing} ? Union{Float64, Missing} : hasmissingtype(S) ? Union{String, Missing} : String
    elseif T === Union{Float64, Missing}
        return S === Int64 || S === Union{Int64, Missing} ? Union{Float64, Missing} : Union{String, Missing}
    elseif hasmissingtype(T) || hasmissingtype(S)
        return Union{String, Missing}
    else
        return String
    end
end

## lazy strings
# bit patterns for missing value, int value, escaped string, position and len in lazy string parsing
const PosLen = UInt64

# primitive type PosLen 64 end
# PosLen(x::UInt64) = Core.bitcast(PosLen, x)
# UInt64(x::PosLen) = Core.bitcast(UInt64, x)

const MISSING_BIT = 0x8000000000000000
missingvalue(x) = (UInt64(x) & MISSING_BIT) == MISSING_BIT

const ESCAPE_BIT = 0x4000000000000000
escapedvalue(x) = (UInt64(x) & ESCAPE_BIT) == ESCAPE_BIT

getpos(x) = (UInt64(x) & 0x3ffffffffff00000) >> 20
getlen(x) = UInt64(x) & 0x00000000000fffff

_unsafe_string(p, len) = ccall(:jl_pchar_to_string, Ref{String}, (Ptr{UInt8}, Int), p, len)

@inline function str(buf, e, poslen)
    missingvalue(poslen) && return missing
    escapedvalue(poslen) && return unescape(PointerString(pointer(buf, getpos(poslen)), getlen(poslen)), e)
    pos, len = getpos(poslen), getlen(poslen)
    return _unsafe_string(pointer(buf, getpos(poslen)), getlen(poslen))
end

@inline function strnomiss(buf, e, poslen)
    escapedvalue(poslen) && return unescape(PointerString(pointer(buf, getpos(poslen)), getlen(poslen)), e)
    pos, len = getpos(poslen), getlen(poslen)
    return _unsafe_string(pointer(buf, getpos(poslen)), getlen(poslen))
end

struct LazyStringVector{T, A <: AbstractVector{PosLen}} <: AbstractVector{T}
    buffer::Vector{UInt8}
    e::UInt8
    poslens::A
end

LazyStringVector{T}(buffer, e, poslens::A) where {T, A} = LazyStringVector{T, A}(buffer, e, poslens)

Base.IndexStyle(::Type{LazyStringVector}) = Base.IndexLinear()
Base.size(x::LazyStringVector) = (length(x.poslens),)

Base.@propagate_inbounds function Base.getindex(x::LazyStringVector{Union{String, Missing}}, i::Int)
    @boundscheck checkbounds(x, i)
    @inbounds s = str(x.buffer, x.e, x.poslens[i])
    return s
end

Base.@propagate_inbounds function Base.getindex(x::LazyStringVector{String}, i::Int)
    @boundscheck checkbounds(x, i)
    @inbounds s = strnomiss(x.buffer, x.e, x.poslens[i])
    return s
end

# optimize iterate for ChainedVector
@inline function Base.iterate(x::LazyStringVector{T}) where {T}
    st = iterate(x.poslens)
    st === nothing && return nothing
    @inbounds s = T === String ? strnomiss(x.buffer, x.e, st[1]) : str(x.buffer, x.e, st[1])
    return s, st[2]
end

@inline function Base.iterate(x::LazyStringVector{T}, state) where {T}
    st = iterate(x.poslens, state)
    st === nothing && return nothing
    @inbounds s = T === String ? strnomiss(x.buffer, x.e, st[1]) : str(x.buffer, x.e, st[1])
    return s, st[2]
end

## column array allocating
# we don't want to use SentinelVector for small integer types due to the higher risk of
# sentinel value collision, so we just use Vector{Union{T, Missing}} and convert to Vector{T} if no missings were found
const SmallIntegers = Union{Int8, UInt8, Int16, UInt16, Int32, UInt32}

# allocate columns for a full file
function allocate(rowsguess, ncols, types, flags, refs)
    columns = Vector{AbstractVector}(undef, ncols)
    for i = 1:ncols
        @inbounds columns[i] = allocate(lazystrings(flags[i]) && (types[i] === String || types[i] === Union{String, Missing}) ? PosLen : types[i], rowsguess)
        if types[i] === PooledString || types[i] ===  Union{PooledString, Missing} || types[i] === CategoricalValue{String, UInt32} || types[i] ===  Union{CategoricalValue{String, UInt32}, Missing}
            refs[i] = RefPool()
        end
    end
    return columns
end

# MissingVector is an efficient representation in SentinelArrays.jl package
allocate(::Type{Union{}}, len) = MissingVector(len)
allocate(::Type{Missing}, len) = MissingVector(len)
function allocate(::Type{PosLen}, len)
    A = Vector{PosLen}(undef, len)
    memset!(pointer(A), typemax(UInt8), sizeof(A))
    return A
end
allocate(::Type{String}, len) = SentinelVector{String}(undef, len)
allocate(::Type{Union{String, Missing}}, len) = SentinelVector{String}(undef, len)
allocate(::Type{PooledString}, len) = Vector{UInt32}(undef, len)
allocate(::Type{Union{PooledString, Missing}}, len) = Vector{UInt32}(undef, len)
allocate(::Type{CategoricalValue{String, UInt32}}, len) = Vector{UInt32}(undef, len)
allocate(::Type{Union{CategoricalValue{String, UInt32}, Missing}}, len) = Vector{UInt32}(undef, len)
allocate(::Type{Bool}, len) = Vector{Union{Missing, Bool}}(undef, len)
allocate(::Type{Union{Missing, Bool}}, len) = Vector{Union{Missing, Bool}}(undef, len)
allocate(::Type{T}, len) where {T <: SmallIntegers} = Vector{Union{Missing, T}}(undef, len)
allocate(::Type{Union{Missing, T}}, len) where {T <: SmallIntegers} = Vector{Union{Missing, T}}(undef, len)
allocate(T, len) = SentinelVector{nonmissingtype(T)}(undef, len)

reallocate!(A, len) = resize!(A, len)
# when reallocating, we just need to make sure the missing bit is set for lazy string PosLen
function reallocate!(A::Vector{PosLen}, len)
    oldlen = length(A)
    resize!(A, len)
    memset!(pointer(A, oldlen + 1), typemax(UInt8), (len - oldlen) * 8)
    return
end

const SVec{T} = SentinelVector{T, T, Missing, Vector{T}}
const SVec2{T} = SentinelVector{T, typeof(undef), Missing, Vector{T}}

ts(T, S) = Core.Compiler.typesubtract(T, S)
# when users pass non-standard types, we need to keep track of them in a Tuple{...} to generate efficient custom parsing kernel codes
function nonstandardtype(T)
    S = ts(ts(ts(ts(ts(ts(ts(ts(ts(T, Int64), Float64), String), PooledString), Bool), Date), DateTime), Time), Missing)
    if S === Union{}
        return S
    elseif S <: SmallIntegers
        return Tuple{Vector{Union{Missing, S}}, S}
    elseif isbitstype(S)
        return Tuple{SVec{S}, S}
    else
        return Tuple{SVec2{S}, S}
    end
end

# one-liner suggested from ScottPJones
consumeBOM(buf, pos) = (length(buf) >= 3 && buf[pos] == 0xef && buf[pos + 1] == 0xbb && buf[pos + 2] == 0xbf) ? pos + 3 : pos

# whatever input is given, turn it into an AbstractVector{UInt8} we can parse with
function getsource(x)
    if x isa AbstractVector{UInt8}
        return x, 1, length(x)
    elseif x isa Base.GenericIOBuffer
        return x.data, x.ptr, x.size
    elseif x isa Cmd || x isa IO
        Base.depwarn("`CSV.File` or `CSV.Rows` with `$(typeof(x))` object is deprecated; pass a filename, `IOBuffer`, or byte buffer directly (via `read(x)`)", :getsource)
        buf = Base.read(x)
        return buf, 1, length(buf)
    else
        try
            buf = Mmap.mmap(string(x))
            return buf, 1, length(buf)
        catch e
            # if we can't mmap, try just `read`ing the whole thing into a byte vector
            buf = read(x)
            return buf, 1, length(buf)
        end
    end
end

getname(buf::Vector{UInt8}) = "<raw buffer>"
getname(cmd::Cmd) = string(cmd)
getname(str) = string(str)
getname(io::I) where {I <: IO} = string("<", I, ">")

# normalizing column name utilities
const RESERVED = Set(["local", "global", "export", "let",
    "for", "struct", "while", "const", "continue", "import",
    "function", "if", "else", "try", "begin", "break", "catch",
    "return", "using", "baremodule", "macro", "finally",
    "module", "elseif", "end", "quote", "do"])

normalizename(name::Symbol) = name
function normalizename(name::String)::Symbol
    uname = strip(Unicode.normalize(name))
    id = Base.isidentifier(uname) ? uname : map(c->Base.is_id_char(c) ? c : '_', uname)
    cleansed = string((isempty(id) || !Base.is_id_start_char(id[1]) || id in RESERVED) ? "_" : "", id)
    return Symbol(replace(cleansed, r"(_)\1+"=>"_"))
end

function makeunique(names)
    set = Set(names)
    length(set) == length(names) && return Symbol[Symbol(x) for x in names]
    nms = Symbol[]
    for nm in names
        if nm in nms
            k = 1
            newnm = Symbol("$(nm)_$k")
            while newnm in set || newnm in nms
                k += 1
                newnm = Symbol("$(nm)_$k")
            end
            nm = newnm
        end
        push!(nms, nm)
    end
    return nms
end

initialtypes(T, x::AbstractDict{String}, names) = Type[haskey(x, string(nm)) ? x[string(nm)] : T for nm in names]
initialtypes(T, x::AbstractDict{Symbol}, names) = Type[haskey(x, nm) ? x[nm] : T for nm in names]
initialtypes(T, x::AbstractDict{Int}, names)    = Type[haskey(x, i) ? x[i] : T for i = 1:length(names)]

initialflags(T, x::AbstractDict{String}, names, lazystrings) = UInt8[haskey(x, string(nm)) ? flag(x[string(nm)], lazystrings) : T for nm in names]
initialflags(T, x::AbstractDict{Symbol}, names, lazystrings) = UInt8[haskey(x, nm) ? flag(x[nm], lazystrings) : T for nm in names]
initialflags(T, x::AbstractDict{Int}, names, lazystrings)    = UInt8[haskey(x, i) ? flag(x[i], lazystrings) : T for i = 1:length(names)]

# given a DateFormat, is it meant for parsing Date, DateTime, or Time?
function timetype(df::Dates.DateFormat)
    date = false
    time = false
    for token in df.tokens
        T = typeof(token)
        if T in (Dates.DatePart{'H'}, Dates.DatePart{'I'}, Dates.DatePart{'M'}, Dates.DatePart{'S'}, Dates.DatePart{'s'})
            time = true
        elseif T in (Dates.DatePart{'y'}, Dates.DatePart{'Y'}, Dates.DatePart{'m'}, Dates.DatePart{'d'}, Dates.DatePart{'u'}, Dates.DatePart{'U'})
            date = true
        end
    end
    return ifelse(date & time, DateTime, ifelse(time, Time, Date))
end

# if a cell value of a csv file has escape characters, we need to unescape it
function unescape(s, e)
    n = ncodeunits(s)
    buf = Base.StringVector(n)
    len = 1
    i = 1
    @inbounds begin
        while i <= n
            b = codeunit(s, i)
            if b == e
                i += 1
                b = codeunit(s, i)
            end
            @inbounds buf[len] = b
            len += 1
            i += 1
        end
    end
    resize!(buf, len - 1)
    return String(buf)
end

"""
    CSV.detect(str::String)

Use the same logic used by `CSV.File` to detect column types, to parse a value from a plain string.
This can be useful in conjunction with the `CSV.Rows` type, which returns each cell of a file as a String.
The order of types attempted is: `Int64`, `Float64`, `Date`, `DateTime`, `Bool`, and if all fail, the input String is returned.
No errors are thrown.
For advanced usage, you can pass your own `Parsers.Options` type as a keyword argument `option=ops` for sentinel value detection.
"""
function detect end

detect(str::String; options=Parsers.OPTIONS) = something(detect(codeunits(str), 1, sizeof(str), options), str)

function detect(buf, pos, len, options)
    int, code, vpos, vlen, tlen = Parsers.xparse(Int64, buf, pos, len, options)
    if Parsers.sentinel(code) && code > 0
        return missing
    end
    if Parsers.ok(code) && vpos + vlen - 1 == len
        return int
    end
    float, code, vpos, vlen, tlen = Parsers.xparse(Float64, buf, pos, len, options)
    if Parsers.ok(code) && vpos + vlen - 1 == len
        return float
    end
    if options.dateformat === nothing
        try
            date, code, vpos, vlen, tlen = Parsers.xparse(Date, buf, pos, len, options)
            if Parsers.ok(code) && vpos + vlen - 1 == len
                return date
            end
        catch e
        end
        try
            datetime, code, vpos, vlen, tlen = Parsers.xparse(DateTime, buf, pos, len, options)
            if Parsers.ok(code) && vpos + vlen - 1 == len
                return datetime
            end
        catch e
        end
    else
        try
            # use user-provided dateformat
            DT = timetype(options.dateformat)
            dt, code, vpos, vlen, tlen = Parsers.xparse(DT, buf, pos, len, options)
            if Parsers.ok(code) && vpos + vlen - 1 == len
                return dt
            end
        catch e
        end
    end
    bool, code, vpos, vlen, tlen = Parsers.xparse(Bool, buf, pos, len, options)
    if Parsers.ok(code) && vpos + vlen - 1 == len
        return bool
    end
    return nothing
end

# a ReversedBuf takes a byte vector and indexes backwards;
# used for the footerskip keyword argument, which starts at the bottom of the file
# and skips lines backwards
struct ReversedBuf <: AbstractVector{UInt8}
    buf::Vector{UInt8}
end

Base.size(a::ReversedBuf) = size(a.buf)
Base.IndexStyle(::Type{ReversedBuf}) = Base.IndexLinear()
Base.getindex(a::ReversedBuf, i::Int) = a.buf[end + 1 - i]

memset!(ptr, value, num) = ccall(:memset, Ptr{Cvoid}, (Ptr{Cvoid}, Cint, Csize_t), ptr, value, num)

# a RefPool holds our refs as a Dict, along with a lastref field which is incremented when a new ref is found while parsing pooled columns
mutable struct RefPool
    refs::Dict{Union{String, Missing}, UInt32}
    lastref::UInt32
end

RefPool() = RefPool(Dict{Union{String, Missing}, UInt32}(), 0)


text:
using Random
abstract type Layer end
#This type isn't Optimized.
abstract type NParam <: Layer end
#This type is Optimized.
abstract type Param <: Layer end
export Param, NParam

include("Pooling.jl")
include("Conv.jl")

"""
    NetWork(layers...)

Connect multiple layers, and build a NeuralNetwork.
NetWork also supports index.
You can also add layers later using the add_layer!() Function.

# Example
```jldoctest
julia> N = NetWork(Dense(10=>5, relu), Dense(5=>1, relu))

julia> N[1]

Dense(IO:10=>5, σ:relu)
```
"""
struct NetWork
    net::Dict{Int, Layer}
    function NetWork()
        new(Dict{Int, Layer}())
    end
end

function NetWork(layers...)
    N = NetWork()
    add_layer!(N, layers...)
    return N
end

function Base.getindex(N::NetWork, i::Int)
    return N.net[i]
end

function (N::NetWork)(x)
    layers = N.net
    for i in 1 : length(layers)
        x = layers[i](x)
    end
    return x
end

function Base.show(io::IO, N::NetWork)
    layers = N.net
    for i in 1 : length(layers)
        print("Layer$i : ")
        if i != length(layers)
            println(layers[i])
        else
            print(layers[i])
        end
    end
end

"""
    Dense(in=>out, σ; set_w = "Xavier", set_b = zeros)
Crate a traditinal `Dense` layer, whose forward propagation is given by:
    y = σ.(W * x .+ b)
The input of `x` should be a Vactor of length `in`, (Sorry for you can't learn using batch. I'll implement)

# Example
```jldoctest
julia> D = Dense(5=>2, relu)
Dense(IO:5=>2, σ:relu)

julia> D(rand(Float64, 5)) |> size
(2,)
```
"""
struct Dense{W, B, F} <: Param
    w::W
    b::B
    σ::F
end

function Dense(io::Pair{<:Integer, <:Integer}, σ; set_w = "Xavier", set_b = zeros)
    w, b = dense_w(io..., set_w), set_b(io[2])
    Dense(w, b, σ)
end

trainable(D::Dense) = D.w, D.b

function Base.show(io::IO, D::Dense)
    o, i = size(D.w)
    σ = D.σ
    print("Dense(IO:$i => $o, σ:$σ)")
end

function (D::Dense)(X::AbstractVecOrMat)
    W, b, σ = D.w, D.b, D.σ
    σ.(muladd(W, X, b))
end

(D::Dense)(x::AbstractArray) = reshape(D(reshape(x, size(x, 1), :)), :, size(x)[2:end]...)

"""
    Flatten()
This layer change the dimentions Image to Vector.

# Example
```jldoctest
julia> F = Flatten()
Flatten(())

julia> F(rand(10, 10, 2, 5)) |> size
(1000, )
```
"""
mutable struct Flatten <: NParam
    csize::Tuple
    function Flatten()
        return new(())
    end
end

function (F::Flatten)(x::Array)
    return reshape(x, :, size(x)[end])
end

"""
    Dropout(p)
This layer dropout the input data.

# Example
```jldoctest
julia> D = Dropout(0.25)
Dropout(0.25)

julia> D(rand(10))
10-element Array{Float64,1}:
 0.0
 0.3955865029078952
 0.8157710047424143
 1.0129613533211907
 0.8060508293474877
 1.1067504108970596
 0.1461289547292684
 0.0
 0.04581776023870532
 1.2794087133638332
```
"""
struct Dropout <: NParam
    p::Float64
    function Dropout(p)
        if 0<=p<=1
            new(p)
        else
            throw(ArhumentError("p must be between 0 and 1!"))
        end
    end
end

function Base.show(io::IO, D::Dropout)
    print("Dropout(")
    print(string(D.p))
    print(")")
end

_dropout_kernel(y::T, p, q) where {T} = (y > p) ? T(1 / q) : T(0)

function dropout_kernel(x, p)
    y = rand!(similar(x))
    y .= _dropout_kernel.(y, p, 1-p)
    return y
end

function (D::Dropout)(x)
    y = dropout_kernel(x, D.p)
    return x .* y
end

@adjoint function (D::Dropout)(x)
    y = dropout_kernel(x, D.p)
    return x.*y, Δ -> (Δ.*y, nothing)
end

"""
     add_layer!(model, layers...)
This function add layers to model. You can add layers when you create a NeuralNetwork with `NetWork()`, and you can also use this function to add layers later.

# Example
```jldoctest
julia> N = NetWork()

julia> add_layer!(N, Dense(10=>5, relu), Dense(5=>1, relu))

julia> N
Layer1 : Dense(IO:10 => 5, σ:relu)
Layer2 : Dense(IO:5 => 1, σ:relu)
```
"""
function add_layer!(model::NetWork, obj::T) where T <: Layer
    l = length(model.net)
    model.net[l+1] = obj
end

function add_layer!(model::NetWork, layers...)
    for layer in layers
        add_layer!(model, layer)
    end
end


text:
# ---
# title: 184. Department Highest Salary
# id: problem184
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories:
# link: <https://leetcode.com/problems/department-highest-salary/description/>
# hidden: true
# ---
#
# The `Employee` table holds all employees. Every employee has an Id, a salary,
# and there is also a column for the department Id.
#
#
#
#     +----+-------+--------+--------------+
#     | Id | Name  | Salary | DepartmentId |
#     +----+-------+--------+--------------+
#     | 1  | Joe   | 70000  | 1            |
#     | 2  | Jim   | 90000  | 1            |
#     | 3  | Henry | 80000  | 2            |
#     | 4  | Sam   | 60000  | 2            |
#     | 5  | Max   | 90000  | 1            |
#     +----+-------+--------+--------------+
#
#
# The `Department` table holds all departments of the company.
#
#
#
#     +----+----------+
#     | Id | Name     |
#     +----+----------+
#     | 1  | IT       |
#     | 2  | Sales    |
#     +----+----------+
#
#
# Write a SQL query to find employees who have the highest salary in each of the
# departments. For the above tables, your SQL query should return the following
# rows (order of rows does not matter).
#
#
#
#     +------------+----------+--------+
#     | Department | Employee | Salary |
#     +------------+----------+--------+
#     | IT         | Max      | 90000  |
#     | IT         | Jim      | 90000  |
#     | Sales      | Henry    | 80000  |
#     +------------+----------+--------+
#
#
# Explanation:
#
# Max and Jim both have the highest salary in the IT department and Henry has
# the highest salary in the Sales department.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
baremodule SimpleSegmentedStacks

function call end

module Internal

using ..SimpleSegmentedStacks: SimpleSegmentedStacks

include("utils.jl")
include("internal.jl")

end  # module Internal

end  # baremodule SimpleSegmentedStacks


text:
include("helpers/listing.jl")

hello_str = "\"Hello, World\""

isinstalled(pkg) = try Pkg.installed(pkg) != nothing catch e false end

md_str = """
md\"\"\"
## This is a heading

*italic*, **bold** and `monospace`.

- one
- two
- three

Code:

```julia
function foo()
    42
end
foo()
```
\"\"\" """

function main(window)
    push!(window.assets, "widgets")
    push!(window.assets, "codemirror")
    push!(window.assets, "tex")

    vbox(
        title(3, "Content"),

        vskip(1em),
        h2("String"),
        "A string is rendered as-is to plain text",
        vskip(1em),
        listing(hello_str),

        vskip(1em),
        h2("Markdown"),
        "Markdown strings get rendered appropriately.",
        vskip(1em),
        listing(md_str),

        vskip(1em),
        h2("Code"),
        md"""To show code with syntax highlighting, you can use the `codemirror` function. Codemirror requires the `"codemirror"` asset.

        Use `push!(window.assets, "codemirror")` to do this.""",
        vskip(1em),
        listing("""
            push!(window.assets, "codemirror")

            codemirror(\"\"\"
                function foo()
                42
                end
                \"\"\"
            )
            """,
            codemirror("""
                function foo()
                    42
                end
                """
            )
        ),
        md"""Note that you also need the `"codemirror"` asset to be loaded even if you are writing code inside `md`.""",

        vskip(1em),
        h2("LaTeX"),
        md"LaTeX strings can be rendered with the `tex` function.",
        vskip(1em),
        listing("""
        tex("(a+b)^2=a^2+b^2+2ab")"""),

        md"[SymPy.jl](https://github.com/jverzani/SymPy.jl) can be used to generate mathematical expressions. Escher automatically typesets them to LaTeX.",
        vskip(1em),
        listing("""
        using SymPy
        Escher.external_setup() # sets up rendering of SymPy symbols

        x = Sym("x")
        SymPy.diff(sin(x^2), x, 5)""",
        isinstalled("SymPy") ? nothing : "SymPy.jl is not installed."),

        h2("Images"),
        md"External images can be included using the `image` function.",
        vskip(1em),
        listing("""
        image(
            "https://upload.wikimedia.org/wikipedia/en/thumb/2/24/Lenna.png/220px-Lenna.png",
            alt="Lenna",
        )
        """),

        md"Images can also be read in from the [`Images.jl`](https://github.com/timholy/Images.jl/) package.",
        vskip(1em),
        listing("""
        using Images
        Escher.external_setup() # sets up rendering of Image objects
        grayim(rand(100,100))
        """,
        isinstalled("Images") ? nothing : "Images.jl is not installed."),

    ) |> pad(2em)
end


text:
using contact_walk
using Test

@testset "contact_walk.jl" begin
    # Write your tests here.
end


text:
using FilePaths
using Test

@testset "FilePaths" begin

include("compat.jl")
include("glob.jl")
include("test_uri.jl")

end


text:
# This dictionary maps easy to remember names to Youtube video IDs
# after adding an ID here, you can use the {{youtube <shortname>}}
# syntax in your markdown files to embed the video into the page!
videos = Dict(
              "test" => ("7827457","BV1Gs411Y7eS","12840293"),
              "test2" => ("8350625","BV18s41187Kw","13743870"),
)


text:
using GetindexArrays
using CoordinateTransformations, Rotations, StaticArrays
using OffsetArrays
using AxisArrays: AxisArrays, AxisArray
using ImageCore, ImageAxes
using Test

if Base.VERSION < v"1.1"
    using Compat   # oneunit(::CartesianIndex)
end

@testset "GetindexArrays.jl" begin
    # A purely computational array
    A = GetindexArray((Base.OneTo(5), Base.OneTo(7))) do y, x
        asin((x-1)^2/72 + (y-1)^2/32)
    end
    @test @inferred(A[1, 1]) == 0
    @test @inferred(A[5, 7]) ≈ π/2
    @test_throws BoundsError A[0, 1]
    @test_throws BoundsError A[6, 1]

    P = [1 2; 3 4]
    @test P[1, 1] == 1
    @test P[1, 2] == 2
    @test P[2, 1] == 3
    @test P[2, 2] == 4
    A = GetindexArray(axes(P)) do idx...
        # reverse both indexes
        P[(size(P) .+ 1 .- idx)...]
    end
    @test A[1, 1] == 4
    @test A[1, 2] == 3
    @test A[2, 1] == 2
    @test A[2, 2] == 1
    @test_throws BoundsError A[3, 2]
    @test_throws BoundsError A[2, 3]
    @test @inferred(axes(A)) == axes(P)

    # Same as above, but inferrably
    makeA(P) = GetindexArray(axes(P)) do idx...
        # reverse both indexes
        P[(size(P) .+ 1 .- idx)...]
    end
    A = @inferred(makeA(P))
    @test @inferred(A[1, 1]) == 4

    # Test that construction is inferred (@inferred can't handle `do` syntax)
    # on an example that mimics PermutedDimsArray
    P = rand(3, 2)
    A = @inferred(GetindexArray{eltype(P)}((i1, i2)->P[i2, i1], reverse(axes(P))))
    @test @inferred(size(A)) == (2, 3)
    @test @inferred(axes(A)) == (1:2, 1:3)
    @test A[1, 1] == P[1, 1]
    @test A[2, 1] == P[1, 2]
    @test_throws BoundsError A[3, 1]
    @test A[1, 2] == P[2, 1]
    @test A[2, 2] == P[2, 2]
    @test A[1, 3] == P[3, 1]
    @test A[2, 3] == P[3, 2]
    @test_throws BoundsError A[1, 4]

    P = rand(3, 5, 4)
    A = GetindexArray((axes(P, 1), axes(P, 3))) do i1, i2
        sum(view(P, i1, :, i2))
    end
    @test size(A) == (3, 4)
    for i in eachindex(A)
        @test A[i] == sum(P[i[1], :, i[2]])
    end

    # A much more complex example using traits, applying a coordinate transformation,
    # mapping values, and local averaging. Also enforces inferrability.

    function special_getindex(A::AbstractArray{T,N}, tforms, i::Vararg{Int,N}) where {T,N}
        # Get the spatial slice and transformation for this time point
        tax, ti = timeaxis(A), i[timedim(A)]
        Aspatial = view(A, tax(ti))
        tform = tforms[ti]
        # Compute the transformed center position
        si = ImageAxes.filter_space_axes(AxisArrays.axes(A), i)
        tformi = tform(SVector(si))
        ci = CartesianIndex(Tuple(round.(Int, tformi)))
        # Average over a spatial 3x3x... region, considering just those points that are in-bounds
        sinds = CartesianIndices(axes(Aspatial))
        s, n = float(zero(eltype(A))), 0
        @inbounds for j in max(first(sinds), ci-oneunit(ci)):min(last(sinds), ci+oneunit(ci))
            s += Aspatial[j]      # cumulative sum
            n += 1                # number of in-bounds voxels
        end
        # square-root transform and return the result
        # Both `s` and `n` will be zero if there were no in-bounds points, in which case we return NaN (which is good)
        return sqrt(s/n)
    end

    # Because tforms is captured, construct everything inside a function
    function make_all(P, tforms)
        # forced specialization
        capture_getindex(i::Vararg{Int,N}) where {N} = special_getindex(P, tforms, i...)
        A = @inferred(GetindexArray(capture_getindex, axes(P)))
        return A
    end
    tforms = [LinearMap(RotMatrix(θ)) for θ = 0:π/16:2π]
    P = zeros(8, 13, length(tforms))
    P[3:6, 4:9, :] .= 1
    P = AxisArray(OffsetArray(P, -3:4, -6:6, 1:length(tforms)), :y, :x, :time)
    A = make_all(P, tforms)
    @test @inferred(A[0, 0, lastindex(tforms)]) === 1.0
end


text:
#!/usr/bin/env julia

import Base.show

struct ParsedLine
  min::Int
  max::Int
  c::Char
  str::String
end

function show(io::IO, parsed_line::ParsedLine)
  print(io, parsed_line.min, " ", parsed_line.max, " ", parsed_line.c, " ", parsed_line.str)
end

function debug_print(parsed_line::ParsedLine)
  println(parsed_line.c, " ", parsed_line.str[parsed_line.min], " ", parsed_line.str[parsed_line.max])
end

function parse_line(line::String)
  m = match(r"(\d+)-(\d+) ([a-z]): (\w+)", line)

  ParsedLine(parse(Int, m.captures[1]), parse(Int, m.captures[2]), m.captures[3][1], m.captures[4])
end

function check_line(parsed_line::ParsedLine)
  parsed_line.min <= count(c -> c == parsed_line.c, parsed_line.str) <= parsed_line.max
end

function check_line2(parsed_line::ParsedLine)
  min_has = parsed_line.str[parsed_line.min] == parsed_line.c
  max_has = parsed_line.str[parsed_line.max] == parsed_line.c

  min_has ⊻ max_has
end

function main()
  res = 0
  for line in eachline("./input")
    parsed_line = parse_line(line)
    res += check_line2(parsed_line) == true
    debug_print(parsed_line)
    println(check_line2(parsed_line))
  end
  println(res)
end

main()


text:
@testset "2_encoders" begin
  @testset "batch_encoder" begin
    parms = EncryptionParameters(SchemeType.bfv)
    poly_modulus_degree = 8192
    set_poly_modulus_degree!(parms, poly_modulus_degree)
    set_coeff_modulus!(parms, coeff_modulus_bfv_default(poly_modulus_degree))
    set_plain_modulus!(parms, plain_modulus_batching(poly_modulus_degree, 20))
    context = SEALContext(parms)

    @testset "first_context_data" begin
      @test_nowarn first_context_data(context)
    end
    context_data = first_context_data(context)

    @testset "qualifiers" begin
      @test_nowarn qualifiers(context_data)
    end
    epq = qualifiers(context_data)

    @testset "using_batches" begin
      @test using_batching(epq) == true
    end

    keygen = KeyGenerator(context)
    public_key_ = PublicKey()
    create_public_key!(public_key_, keygen)
    secret_key_ = secret_key(keygen)
    relin_keys_ = RelinKeys()
    create_relin_keys!(relin_keys_, keygen)
    encryptor = Encryptor(context, public_key_)
    evaluator = Evaluator(context)
    decryptor = Decryptor(context, secret_key_)

    @testset "BatchEncoder" begin
      @test_nowarn BatchEncoder(context)
    end
    batch_encoder = BatchEncoder(context)

    @testset "slot_count" begin
      @test slot_count(batch_encoder) == 8192
    end
    slot_count_ = slot_count(batch_encoder)
    row_size = div(slot_count_, 2)

    pod_matrix = zeros(UInt64, slot_count_)
    pod_matrix[1] = 0
    pod_matrix[2] = 1
    pod_matrix[3] = 2
    pod_matrix[4] = 3
    pod_matrix[row_size + 1] = 4
    pod_matrix[row_size + 2] = 5
    pod_matrix[row_size + 3] = 6
    pod_matrix[row_size + 4] = 7

    plain_matrix = Plaintext()
    @testset "encode!" begin
      @test_nowarn encode!(plain_matrix, pod_matrix, batch_encoder)
    end

    pod_result = similar(pod_matrix)
    @testset "decode!" begin
      @test_nowarn decode!(pod_result, plain_matrix, batch_encoder)
    end

    # Extra: encode/decode for Int64 data type
    pod_matrix_int64 = Int64.(pod_matrix)
    plain_matrix_int64 = Plaintext()
    @testset "encode!" begin
      @test_nowarn encode!(plain_matrix_int64, pod_matrix_int64, batch_encoder)
    end
    pod_result_int64 = similar(pod_matrix_int64)
    @testset "decode!" begin
      @test_nowarn decode!(pod_result_int64, plain_matrix_int64, batch_encoder)
    end

    encrypted_matrix = Ciphertext()
    @testset "encrypt!" begin
      @test_nowarn encrypt!(encrypted_matrix, plain_matrix, encryptor)
    end

    @testset "noise budget 1" begin
      @test invariant_noise_budget(encrypted_matrix, decryptor) in (145, 146, 147)
    end

    pod_matrix2 = ones(UInt64, slot_count_)
    pod_matrix2[2:2:slot_count_] .= 2
    plain_matrix2 = Plaintext()

    @testset "encode!" begin
      @test_nowarn encode!(plain_matrix2, pod_matrix2, batch_encoder)
    end

    @testset "sum, square, and relinearize" begin
      @test_nowarn add_plain_inplace!(encrypted_matrix, plain_matrix2, evaluator)
      @test_nowarn square_inplace!(encrypted_matrix, evaluator)
      @test_nowarn relinearize_inplace!(encrypted_matrix, relin_keys_, evaluator)
    end

    @testset "noise budget 2" begin
      @test invariant_noise_budget(encrypted_matrix, decryptor) in (113, 114, 115)
    end

    plain_result = Plaintext()
    @testset "decrypt! and decode!" begin
      @test_nowarn decrypt!(plain_result, encrypted_matrix, decryptor)
      @test_nowarn decode!(pod_result, plain_result, batch_encoder)
    end
  end
end


text:
function generaSupBruta(ps_predio, conjuntoLadosCalle, anchoEspacioPublico)

    ps_bruto = polyShape.polyExpandSides_v2(ps_predio, anchoEspacioPublico/2, conjuntoLadosCalle)
    ps_publico = polyShape.polyExpandSides_v2(ps_predio, anchoEspacioPublico, conjuntoLadosCalle)

    return ps_bruto, ps_publico
end


text:
hexMap = Dict{Char, String}(
    '0' => "0000",
    '1' => "0001",
    '2' => "0010",
    '3' => "0011",
    '4' => "0100",
    '5' => "0101",
    '6' => "0110",
    '7' => "0111",
    '8' => "1000",
    '9' => "1001",
    'A' => "1010",
    'B' => "1011",
    'C' => "1100",
    'D' => "1101",
    'E' => "1110",
    'F' => "1111",
    '\n' => ""
)

function hexToBin(str)
    return join([hexMap[c] for c in str])
end

function binToDec(str)
    len = length(str)
    return sum([parse(Int, c) * 2^(len - i) for (i, c) in enumerate(str)])
end

function parsePacket(bits, ind=1)
    version = bits[ind:ind+2]
    typeId = bits[ind+3:ind+5]

    ind += 6

    if typeId == "100"
        # Literal Value
        valStr = ""

        for index in ind:5:length(bits)
            valStr *= bits[index+1:index+4]
            ind += 5
            if bits[index] == '0'
                dec = binToDec(valStr)
                return (version = binToDec(version), typeId = binToDec(typeId), literal = dec, index = ind, innerPackets = [])
            end
        end
    else
        # Operator
        lengthTypeId = bits[ind]
        ind += 1
        if lengthTypeId == '0'
            subPacketslengthBits = binToDec(bits[ind:ind+14]) # Next 15 bits
            ind += 15
            subPackets = []
            oldInd = ind
            while true
                packet = parsePacket(bits, ind)
                push!(subPackets, packet)
                ind = packet[:index]

                if ind - oldInd == subPacketslengthBits
                    break
                end
            end
            return (version = binToDec(version), typeId = binToDec(typeId), literal = nothing, index = ind, innerPackets = subPackets)
        else
            nSubPackets = binToDec(bits[ind:ind+10]) # Next 11 bits
            ind += 11
            subPackets = []
            for _ in 1:nSubPackets
                packet = parsePacket(bits, ind)
                ind = packet[:index]
                push!(subPackets, packet)
            end
            return (version = binToDec(version), typeId = binToDec(typeId), literal = nothing, index = ind, innerPackets = subPackets)
        end
    end
end

packet = parsePacket(hexToBin(read("input.txt", String)))

#! Part 1

function sumVersions(packet)
    version = packet[:version]
    for pckt in packet[:innerPackets]
        version += sumVersions(pckt)
    end
    return version
end

sumVersions(packet)

#! Part 2
function evalPacket(packet)
    typeId = packet[:typeId]
    if typeId == 0
        # Sum
        return sum([evalPacket(p) for p in packet[:innerPackets]])
    elseif typeId == 1
        # Product
        return prod([evalPacket(p) for p in packet[:innerPackets]])
    elseif typeId == 2
        # Minimum
        return minimum([evalPacket(p) for p in packet[:innerPackets]])
    elseif typeId == 3
        # Maximum
        return maximum([evalPacket(p) for p in packet[:innerPackets]])
    elseif typeId == 4
        # Literal
        return packet[:literal]
    elseif typeId == 5
        # Greater than packet
        p1 = evalPacket(packet[:innerPackets][1])
        p2 = evalPacket(packet[:innerPackets][2])
        return p1 > p2 ? 1 : 0
    elseif typeId == 6
        # Less than packet
        p1 = evalPacket(packet[:innerPackets][1])
        p2 = evalPacket(packet[:innerPackets][2])
        return p1 < p2 ? 1 : 0
    elseif typeId == 7
        # Equality packet
        p1 = evalPacket(packet[:innerPackets][1])
        p2 = evalPacket(packet[:innerPackets][2])
        return p1 == p2 ? 1 : 0
    end
end

evalPacket(packet)

max_stars_repo_path: Day 16/main.jl

max_issues_repo_path: Day 16/main.jl

max_forks_repo_path: Day 16/main.jl


text:
# verbose flag
const VERBOSE = Ref(false)

"""
    set_verbose!(flag::Bool)

For debugging - set flag to print macro expansions
"""
function set_verbose!(verbose::Bool)
    VERBOSE[] = verbose
end


text:
using GeoMakie, Makie, Proj4

lons = LinRange(0.5, 359.5, 360)
lats = LinRange(-89.5, 89.5, 180)

field = [exp(cosd(l)) + 3(y/90) for l in lons, y in lats]

cf = circshift(field, 180) # shift the field to the correct position

source = Projection("+proj=lonlat +lon_0=180 +pm=180")
dest = Projection("+proj=moll +lon_0=0")

xs, ys = xygrid(lons, lats)
Proj4.transform!(source, dest, vec(xs), vec(ys))

scene = surface(xs, ys; color = cf, shading = false, show_axis = false)

geoaxis!(scene, -180, 180, -90, 90; crs = (src = source, dest = dest,))

coastlines!(scene, 1; crs = (src = source, dest = dest,))


text:
"""
    _foldl_blockarray(rf, acc, coll::BlockArrays.BlockArray)
"""
function _foldl_blockarray(rf::RF, acc, coll) where {RF}
    # `isempty` is required for `map(first, blockaxes(coll))` below
    isempty(coll) && return complete(rf, acc)
    return complete(rf, @return_if_reduced _foldl_blockarray(
        rf, acc, coll,
        (),  # (outer) `Block`s
        (),  # (outer) offset indices
        Val(ndims(coll) - 1),  # stop recursion when `== Val(length(block))`
    ))
end

# Looping over outer/right dimensions; i.e., one of `b`, `c`, `d` in
# `coll[a, b, c, d]` but not `a`.
@inline function _foldl_blockarray(rf::RF, acc, coll, block, offset, ValN₋₁) where {RF}
    blockaxes = BlockArrays.blockaxes

    d = ndims(coll) - length(block)
    sample_block_template = map(first, blockaxes(coll))
    for b in blockaxes(coll, d)
        sample_block = @set sample_block_template[d] = b
        for k in axes(coll[sample_block...], d)
            acc = @return_if_reduced _foldl_blockarray(
                rf, acc, coll, (b, block...), (k, offset...), ValN₋₁,
            )
        end
    end
    return acc
end

# Finally, looping over the inner/left -most dimension:
@inline function _foldl_blockarray(
    rf::RF,
    acc,
    coll,
    block::NTuple{N₋₁, Any},
    offset::NTuple{N₋₁, Int},
    ::Val{N₋₁},
) where {RF,N₋₁}

    blockaxes = BlockArrays.blockaxes

    @inbounds for b in blockaxes(coll, 1)
        array = coll[b, block...]
        @simd_if rf for k in axes(array, 1)
            acc = @next(rf, acc, array[k, offset...])
        end
    end
    return acc
end

# TODO: write reduce for BlockArrays which can be done in the "natural" order


text:
@testset "empty" begin
    @rules Demo [] begin end
    @test Demo(a) == a
    @test Demo(f(a)) == f(a)
end

@testset "free, simple" begin
    @rules Demo [x] begin
        f(x) := g(x)
    end
    @test Demo(f(a)) == g(a)
    @test Demo(h(f(a))) == h(g(a))
    @test Demo(g(a)) == g(a)
    @test Demo(h(a)) == h(a)
end

@testset "free, complex" begin
    @testset "orthogonal" begin
        @rules Demo [x, y] begin
            f(x)       := g(x)
            h(x, x, y) := p(x, y)
        end

        @test Demo(f(a)) == g(a)
        @test Demo(h(a, a, b)) == p(a, b)
        @test Demo(h(a, b, b)) == h(a, b, b)
    end

    @testset "overlapping" begin
        @rules Demo [x] begin
            f(x)    := g(x)
            f(g(x)) := h(x)
        end

        @test Demo(f(a)) == g(a)
        @test Demo(f(g(a))) ∈ [g(g(a)), h(a)]
        @test Demo(g(a)) == g(a)
    end

    @testset "identical" begin
        @rules Demo [x] begin
            f(x) := g(x)
            f(x) := h(x)
        end

        @test Demo(f(a)) ∈ [g(a), h(a)]
        @test Demo(f(f(a))) ∈ [g(g(a)), g(h(a)), h(g(a)), h(h(a))]
    end
end

@testset "commutative, simple" begin
    @rules Demo [x] begin
        p(x, x) := x
    end

    @test Demo(p(a, a)) == a
    @test Demo(p(b, b)) == b
    @test Demo(p(a, b)) == p(a, b)
end

@testset "commutative, complex" begin
    @testset "orthogonal" begin
        @rules Demo [x] begin
            p(a, x)   := x
            p(b, b) := c
        end

        @test Demo(p(a, b)) == b
        @test Demo(f(p(a, b))) == f(b)
    end

    @testset "overlapping" begin
        @rules Demo [x] begin
            p(a, x) := a
            p(x, b) := b
        end

        @test Demo(p(a, c)) == a
        @test Demo(p(b, c)) == b
        @test Demo(p(a, b)) ∈ [a, b]
    end

    @testset "identical" begin
        @rules Demo [x, y] begin
            p(x, y) := q(x, x)
            q(x, x) := a
            q(x, x) := b
        end

        @test Demo(p(c, c)) ∈ [a, b]
    end
end

@testset "sample cases" begin
    @testset "boolean logic" begin
        @theory! begin
            F   => FreeTheory()
            T   => FreeTheory()
            and => CTheory()
            or  => CTheory()
            not => FreeTheory()
        end

        @rules Prop [x, y] begin
            and(x, F) := F
            and(x, T) := x
            and(x, x) := x

            or(x, F) := x
            or(x, T) := T
            or(x, x) := x

            not(and(x, y)) := or(not(x), not(y))
            not(or(x, y))  := and(not(x), not(y))
            not(not(x))    := x
        end

        @test @rewrite(Prop, and(A, not(not(A)))) == @term(A)
        @test @rewrite(Prop, and(or(F, F), B)) == @term(F)
    end

    @testset "natural arithmetic" begin
        @theory! begin
            z => FreeTheory()
            s => FreeTheory()
            add => FreeTheory()
            mul => FreeTheory()
        end

        @rules Arithmetic [x, y] begin
            add(x, z)    := x
            add(x, s(y)) := s(add(x, y))
            mul(x, z)    := z
            mul(x, s(y)) := add(x, mul(x, y))
        end

        _nat(n) = n == 0 ? @term(z) : @term(s($(_nat(n - 1))))

        @testset "$x * $y" for x ∈ 0:5, y ∈ 0:5
            nx, ny = _nat(x), _nat(y)
            @test @rewrite(Arithmetic, mul($nx, $ny)) == _nat(x * y)
        end
    end

    @testset "commutative addition" begin
        @theory! begin
            z => FreeTheory()
            s => FreeTheory()
            (++) => CTheory()
        end

        @rules Addition [x, y] begin
            x ++ z    := x
            s(x) ++ y := s(x ++ y)
        end

        _nat(n) = n == 0 ? @term(z) : @term(s($(_nat(n - 1))))

        @test @term(a ++ b) == @term(b ++ a)
        @test @rewrite(Addition, a ++ s(b)) == @term(s(a ++ b))
        @test @rewrite(Addition, s(a) ++ b) == @term(s(a ++ b))
        @test @rewrite(Addition, s(a) ++ s(b)) == @term(s(s(a ++ b)))
        @testset "$x ++ $y" for x ∈ 0:5, y ∈ 0:5
            nx, ny = _nat(x), _nat(y)
            @test @rewrite(Addition, $nx ++ $ny) == _nat(x + y)
        end
    end

    @testset "list reverse" begin
        @theory! begin
            nil  => FreeTheory()
            cons => FreeTheory()
            rev     => FreeTheory()
            rev_aux => FreeTheory()
        end

        @rules Reverse [x, y, z] begin
            rev_aux(nil, x)        := x
            rev_aux(cons(x, y), z) := rev_aux(y, cons(x, z))
            rev(x)                 := rev_aux(x, nil)
        end

        function _from(arr)
            l = @term(nil)
            for x ∈ reverse(arr)
                l = @term(cons($x, $l))
            end
            return l
        end

        @testset "reverse length $n" for n ∈ 0:10
            local list
            arr = [@term($(Symbol(i))) for i ∈ 1:n]
            @test @rewrite(Reverse, rev($(_from(arr)))) == _from(reverse(arr))
        end
    end

    @testset "symbolic differentiation" begin
        @theory! begin
            ∂ => FreeTheory()
            z => FreeTheory()
            s => FreeTheory()
            (+) => ACTheory()
            (*) => ACTheory()
            sin => FreeTheory()
            cos => FreeTheory()
        end

        @rules Diff [x, y, t] begin
            x + z := x
            x + s(y) := s(x + y)

            x * z := z
            x * s(y) := x + x * y

            ∂(x, x) := s(z)

            ∂(z, t) := z
            ∂(s(x), t) := ∂(x, t)
            ∂(x + y, t) := ∂(x,t) + ∂(y,t)
            ∂(x * y, t) := ∂(x,t)*y + x*∂(y,t)
            ∂(sin(x), t) := cos(x) * ∂(x, t)
       end

       @test @rewrite(Diff, ∂(sin(s(s(z())) * x()) + x + x + x, x())) ==
             @rewrite(Diff, s(s(z))*cos(s(s(z)) * x) + s(s(s(z))))
    end
end


text:
module AgentBasedModels

using DataFrames: AbstractAggregate, getiterator
using CUDA: findfirst, atomic_add!
using Base: add_with_overflow
using Random
using LinearAlgebra
using Distributions
using CUDA
using DataFrames
using CSV
import MacroTools: postwalk, @capture, inexpr, prettify, gensym_ids, flatten, unblock
#using WriteVTK

export MediumFlat, δMedium_
export configurator_

#Constants
include("./constants/constants.jl")

#Agent
export @agent, Agent
export BoundaryFlat, Periodic, Bounded, Free
export PeriodicBoundaryCondition, DirichletBoundaryCondition, DirichletBoundaryCondition_NewmannBoundaryCondition, NewmannBoundaryCondition, NewmannBoundaryCondition_DirichletBoundaryCondition
include("./agent/boundary/boundaryAbstract.jl")
include("./agent/boundary/boundaryFlatStructures.jl")
include("./agent/agentStructure.jl")
include("./agent/constructAgent.jl")
include("./model/structProgram.jl") #Structure
include("./agent/boundary/boundaryFlatFunctions.jl")

#Model
export Model, compile
include("./model/model.jl")
include("./model/agentCode/basic/addParameters.jl")
include("./model/agentCode/basic/checkBounds.jl")
include("./model/agentCode/basic/cleanInteraction.jl")
include("./model/agentCode/basic/eventRemoveAgent.jl")
include("./model/agentCode/basic/eventAddAgent.jl")
include("./model/agentCode/basic/updateGlobal.jl")
include("./model/agentCode/basic/updateLocal.jl")
include("./model/agentCode/basic/updateLocalInteraction.jl")
include("./model/compile.jl")

#Community
export Community, CommunityInTime, saveCSV, loadCommunityFromCSV, loadCommunityInTimeFromCSV
export initialiseCommunityCompactHexagonal, initialiseCommunityCompactCubic
include("./community/community.jl")
include("./community/baseModuleExtensions.jl")
include("./community/constructors/latices/compactHexagonal.jl")
include("./community/constructors/latices/cubic.jl")
include("./community/constructors/extrude.jl")
include("./community/constructors/initialisers.jl")
include("./community/IO/save.jl")
include("./community/IO/load.jl")

#Random
include("./model/random/distribution.jl")
include("./model/random/randomAdapt.jl")

#Neighbors
include("./model/agentCode/neighbors/neighborsFull.jl")
include("./model/agentCode/neighbors/neighborsGrid.jl")
include("./model/agentCode/neighbors/neighbors.jl")

#Integrators
include("./model/agentCode/integrator/euler.jl")
include("./model/agentCode/integrator/heun.jl")
include("./model/agentCode/integrator/integrators.jl")

#Integrators Medium
include("./model/agentCode/integratorMedium/ftcs.jl")
include("./model/agentCode/integratorMedium/lax.jl")
include("./model/agentCode/integratorMedium/leapfrog.jl")
include("./model/agentCode/integratorMedium/integratorsMedium.jl")
include("./model/agentCode/basic/updateMediumInteraction.jl")

#Saving
include("./model/agentCode/saving/saveRAM.jl")
include("./model/agentCode/saving/saveCSV.jl")
include("./model/agentCode/saving/saving.jl")

#Cuda specific functions
include("./model/cuda/cudaAdapt.jl")
include("./model/cuda/cudaConfigurator.jl")

#Auxiliar function
include("./auxiliar/checkDeclared.jl")
include("./auxiliar/clean.jl")
include("./auxiliar/substitution.jl")
include("./auxiliar/vectorize.jl")
include("./auxiliar/arguments.jl")
include("./auxiliar/wrapping.jl")
include("./auxiliar/emptyquote.jl")
include("./auxiliar/updates.jl")
include("./auxiliar/substitution2.jl")
include("./auxiliar/extract.jl")

end


text:
# machine code generation

function machine(cap::VersionNumber, triple::String)
    InitializeNVPTXTarget()
    InitializeNVPTXTargetInfo()
    t = Target(triple)

    InitializeNVPTXTargetMC()
    cpu = "sm_$(cap.major)$(cap.minor)"
    if cuda_driver_version >= v"9.0" && v"6.0" in ptx_support
        # in the case of CUDA 9, we use sync intrinsics from PTX ISA 6.0+
        feat = "+ptx60"
    else
        feat = ""
    end
    tm = TargetMachine(t, triple, cpu, feat)
    asm_verbosity!(tm, true)

    return tm
end

function mcgen(ctx::CompilerContext, mod::LLVM.Module, f::LLVM.Function)
    tm = machine(ctx.cap, triple(mod))

    InitializeNVPTXAsmPrinter()
    return String(emit(tm, mod, LLVM.API.LLVMAssemblyFile))
end


text:
using JLD
# @load "BF-Journal_result_6scen_delayOff.jld"
@load "BF-Journal_result_6scen_delayOn.jld"

status
obj_value
simulation_time
ESlocation
SoC
u
z
pg
qg
pls
qls
IC
OC
EC
OCgen
OCloadshed
OCstorage
x
delta
mu
sigma_d
sigma_l

N_storage = size(u,1)
N_bus = size(u,2)
N_time = size(u,3)
N_scen = size(u,4)

ESlocation_normal = Array{Int64}(N_storage)
temp_z = z
temp_z[find(temp_z.<0.5)]=0
temp_z[find(temp_z.>0.5)]=1
temp_z_idx = ind2sub(temp_z,find(temp_z.==1))
for k = 1:N_storage
	ESlocation_normal[k] = temp_z_idx[2][k] - 1
end

ESlocation = Array{Int64}(N_storage, N_time, N_scen)
for s = 1:N_scen
	for k = 1:N_storage
		temp = u[k,:,:,s]
		temp[find(temp.<0.5)]=0
		temp[find(temp.>0.5)]=1
		temp_idx = ind2sub(temp,find(temp.==1))
		temp_location = Array{Int64}(N_time,1)
		temp_location[1:end] = -1
		for i = 1:length(temp_idx[1])
			temp_location[temp_idx[2][i]]  = temp_idx[1][i] - 1
		end
		ESlocation[k,:,s] = temp_location
	end
end

firstbus = Array{Int64}(N_storage, N_time, N_scen)
lastbus = Array{Int64}(N_storage, N_time, N_scen)
for k in 1:N_storage
	for t in 1:N_time
		for s in 1:N_scen
			firstbus[k,t,s] = findfirst(delta[k,:,t,s],1) - 1
			lastbus[k,t,s] = findlast(delta[k,:,t,s],1) - 1
		end
	end
end

sigma_d[b,t,s];
sigma_l[l,t,s];
ind2sub(sigma_d,find(x->x==0,sigma_d))
ind2sub(sigma_d[:,:,2],find(x->x==0,sigma_d[:,:,2]))
ind2sub(sigma_l,find(x->x==0,sigma_l))

writecsv("result_sigma_d.csv", sigma_d)


text:
module Mongo

include("BSON.jl")
using Mongo.BSON

export UPSERT, MULTI,
       find, find_one, count, update, insert, remove, run_command

const MONGO_LIB = "libmongoc"
const MONGO_OK = 0
const MONGO_ERROR = -1

# Update flags
const UPSERT = 1
const MULTI  = 2

include("mongo_client.jl")
include("mongo_cursor.jl")

function find(client::MongoClient, namespace::String, query, fields, limit::Int, skip::Int)
    MongoCursor(client, namespace, BSONObject(query), BSONObject(fields), limit, skip)
end
find(client::MongoClient, namespace::String, query, fields) = find(client, namespace, query, fields, 0, 0)
find(client::MongoClient, namespace::String, query) = find(client, namespace, query, BSONObject(), 0, 0)
find(client::MongoClient, namespace::String) = find(client, namespace, BSONObject(), BSONObject(), 0, 0)

function find_one(client::MongoClient, namespace::String, query, fields)
    bson = BSONObject()
    errno = ccall((:mongo_find_one, MONGO_LIB), Int32, (Ptr{Void},Ptr{Uint8},Ptr{Void},Ptr{Void},Ptr{Void}),
                client._mongo, bytestring(namespace), BSONObject(query)._bson, BSONObject(fields)._bson, bson._bson)
    errno == MONGO_ERROR ? nothing : dict(bson)
end
find_one(client::MongoClient, namespace::String, query) = find_one(client, namespace, query, BSONObject())
find_one(client::MongoClient, namespace::String) = find_one(client, namespace, BSONObject(), BSONObject())

function count(client::MongoClient, namespace::String, query)
    db, collection = split(namespace, '.')
    c = ccall((:mongo_count, MONGO_LIB), Float64, (Ptr{Void}, Ptr{Uint8}, Ptr{Uint8}, Ptr{Void}),
                client._mongo, bytestring(db), bytestring(collection), BSONObject(query)._bson)
    if c == MONGO_ERROR
        error("Unable to perform Mongo count – mongo_count()")
    end
    convert(Int, c)
end
count(client::MongoClient, namespace::String) = count(client, namespace, BSONObject())

function update(client::MongoClient, namespace::String, query, op, flags::Int)
    errno = ccall((:mongo_update, MONGO_LIB), Int32, (Ptr{Void}, Ptr{Uint8}, Ptr{Void}, Ptr{Void}, Int32, Ptr{Void}),
                client._mongo, bytestring(namespace), BSONObject(query)._bson, BSONObject(op)._bson, convert(Int32, flags), C_NULL)
    if errno == MONGO_ERROR
        error("Unable to update document(s) – mongo_update()")
    end
end
update(client::MongoClient, namespace::String, query, op) = update(client, namespace, query, op, 0)

function insert(client::MongoClient, namespace::String, bson)
    errno = ccall((:mongo_insert, MONGO_LIB), Int32, (Ptr{Void}, Ptr{Uint8}, Ptr{Void}, Ptr{Void}),
                client._mongo, bytestring(namespace), BSONObject(bson)._bson, C_NULL)
    if errno == MONGO_ERROR
        error("Unable to insert document(s) – mongo_insert()")
    end
end

function remove(client::MongoClient, namespace::String, query)
    errno = ccall((:mongo_remove, MONGO_LIB), Int32, (Ptr{Void}, Ptr{Uint8}, Ptr{Void}, Ptr{Void}),
        client._mongo, bytestring(namespace), BSONObject(query)._bson, C_NULL)
    if errno == MONGO_ERROR
        error("Unable to remove document(s) – mongo_remove()")
    end
end

function run_command(client::MongoClient, dbname::String, command)
    bson = BSONObject()
    errno = ccall((:mongo_run_command, MONGO_LIB), Int32,
                  (Ptr{Void}, Ptr{Uint8}, Ptr{Void}, Ptr{Void}),
                  client._mongo, bytestring(dbname), BSONObject(command)._bson, bson._bson)
    errno == MONGO_ERROR ? nothing : dict(bson)
end

end # module Mongo


text:
using ITensors
using ITensorMPOCompression

function make_2D_Heisenberg_AutoMPO(sites,Nx::Int64,Ny::Int64,hz::Float64=0.0,J::Float64=1.0)
    N=length(sites)
    @assert N==Nx*Ny
    ampo = OpSum()
    for j=1:N
        add!(ampo, hz   ,"Sz", j)
    end
    lattice = square_lattice(Nx, Ny; yperiodic = false)

    # Define the Heisenberg spin Hamiltonian on this lattice
    ampo = OpSum()
    for b in lattice
        ampo .+= 0.5*J, "S+", b.s1, "S-", b.s2
        ampo .+= 0.5*J, "S-", b.s1, "S+", b.s2
        ampo .+=     J, "Sz", b.s1, "Sz", b.s2
    end
    H = MPO(ampo,sites)
    fix_autoMPO!(H) #swap row[2]<->row[Dw] and col[2]<->col[Dw]
    return H
end
using Printf
Base.show(io::IO, f::Float64) = @printf(io, "%1.3f", f)
Nx=10
Ny=6
N=Nx*Ny; #10 sites
sites = siteinds("S=1/2",N);
H=make_2D_Heisenberg_AutoMPO(sites,Nx,Ny)
Dw_auto=get_Dw(H)
ss=truncate!(H)
Dw_trunc=get_Dw(H)
if Dw_auto==Dw_trunc
    println("You can't beat autoMPO!!!!")
else
    println("Wa-hoo truncate did something useful :)")
end


text:
# Julia 1.0: Type definition not allowed inside a local scope, therefore we must define our
# TestLayer outside of the @testset
abstract type TestLayer{Next <: HTTP.Layer} <: HTTP.Layer{Next} end

@testset "AWS Signature Version 4" begin
    r = @SymDict(
        creds = AWSCredentials("AKIDEXAMPLE","wJalrXUtnFEMI/K7MDENG+bPxRfiCYEXAMPLEKEY"),
        region = "us-east-1",
        verb = "POST",
        service = "iam",
        url = "http://iam.amazonaws.com/",
        content = "Action=ListUsers&Version=2010-05-08",
        headers = Dict(
            "Content-Type" => "application/x-www-form-urlencoded; charset=utf-8",
            "Host" => "iam.amazonaws.com"
        )
    )

    AWSCore.sign!(r, DateTime("2011-09-09T23:36:00"))

    h = r[:headers]
    out = join(["$k: $(h[k])\n" for k in sort(collect(keys(h)))])

    expected = (
        "Authorization: AWS4-HMAC-SHA256 " *
        "Credential=AKIDEXAMPLE/20110909/us-east-1/iam/aws4_request, " *
        "SignedHeaders=content-md5;content-type;host;" *
        "x-amz-content-sha256;x-amz-date, " *
        "Signature=1a6db936024345449ef4507f890c5161bbfa2ff2490866653bb8b58b7ba1554a\n" *
        "Content-MD5: r2d9jRneykOuUqFWSFXKCg==\n" *
        "Content-Type: application/x-www-form-urlencoded; charset=utf-8\n" *
        "Host: iam.amazonaws.com\n" *
        "x-amz-content-sha256: b6359072c78d70ebee1e81adcbab4f01bf2c23245fa365ef83fe8f1f955085e2\n" *
        "x-amz-date: 20110909T233600Z\n"
    )

    @test out == expected
end

@testset "HTTP Request - AWS4AuthLayer" begin
    test_access_key = "TEST_ACCESS_KEY"
    test_secret_key = "TEST_SECRET_KEY"

    function HTTP.request(::Type{TestLayer{Next}}, io::IO, req, body; kw...) where Next
        @test kw[:aws_access_key_id] == test_access_key
        @test kw[:aws_secret_access_key] == test_secret_key
        return HTTP.request(Next, io, req, body; kw...)
    end

    function _create_stack()
        custom_stack = insert(stack(), StreamLayer, TestLayer)
        custom_stack = insert(custom_stack, RetryLayer, SignatureV4.AWS4AuthLayer)
        result = HTTP.request(custom_stack, "GET", "http://httpbin.org/ip")
        @test result.status == 200
    end

    @testset "Environment Variables" begin
        withenv(
            "AWS_ACCESS_KEY_ID" => test_access_key,
            "AWS_SECRET_ACCESS_KEY" => test_secret_key,
        ) do
            _create_stack()
        end
    end

    @testset "Credentials File - Default" begin
        mktemp() do creds_file, creds_io
            write(creds_io, """
                [default]
                aws_access_key_id=$(test_access_key)
                aws_secret_access_key=$(test_secret_key)
            """)
            close(creds_io)

            withenv(
                "AWS_ACCESS_KEY_ID" => nothing,
                "AWS_SECRET_ACCESS_KEY" => nothing,
                "AWS_SHARED_CREDENTIALS_FILE" => creds_file,
            ) do
                _create_stack()
            end
        end
    end

    @testset "Credentials File - Specified Profile" begin
        aws_profile = "test"

        mktemp() do creds_file, creds_io
            write(creds_io, """
                [$(aws_profile)]
                aws_access_key_id=$(test_access_key)
                aws_secret_access_key=$(test_secret_key)
            """)
            close(creds_io)

            withenv(
                "AWS_PROFILE" => aws_profile,
                "AWS_ACCESS_KEY_ID" => nothing,
                "AWS_SECRET_ACCESS_KEY" => nothing,
                "AWS_SHARED_CREDENTIALS_FILE" => creds_file,
            ) do
                _create_stack()
            end
        end
    end

    @testset "Configuration File" begin
        mktemp() do config_file, config_io
            write(config_io, """
                [default]
                aws_access_key_id=$(test_access_key)
                aws_secret_access_key=$(test_secret_key)
            """)
            close(config_io)

            withenv(
                "AWS_ACCESS_KEY_ID" => nothing,
                "AWS_SECRET_ACCESS_KEY" => nothing,
                "AWS_CONFIG_FILE" => config_file,
            ) do
                cred_patch = @patch dot_aws_credentials_file() = ""

                apply([cred_patch]) do
                    _create_stack()
                end
            end
        end
    end
end


text:
#!/usr/bin/env julia

using VibronicToolkit

using DelimitedFiles

using ArgParse

s = ArgParseSettings()
s.autofix_names = true
@add_arg_table s begin
    "--conf"
        metavar = "FILE"
        help = "path to config file"
        required = true
    "--max-iter"
        metavar = "ITER"
        help = "maximum number of iterations"
        arg_type = Int
        default = 100
    "--out-conf"
        metavar = "FILE"
        help = "path to system configuration output file"
        required = true
    "--out-vs"
        metavar = "FILE"
        help = "path to transformation vector output file"
end
c = parse_args(ARGS, s, as_symbols=true)

sys = read(c[:conf], DenseSystem)
max_iter = c[:max_iter]
out_conf = c[:out_conf]
out_vs = c[:out_vs]

decomp = IterativeDecomposition(sys, max_iter)

write(out_conf, DiagonalSystem(decomp))
if out_vs !== nothing
    writedlm(out_vs, permutedims(decomp.vs))
end