text:
export initAddEcho
export addEcho

"""
    arcrc(T₆₀::Number, Lᵣ::NTuple{3,Number})

                  β * V
`T₆₀ = ———————————————————————`\n
          - S * ln(1-α) + m₄ * V
where β = 0.160465, m₄ = 0.0104.
V is volume unit m³, S is area unit m²
"""
function arcrc(T₆₀::Number, Lᵣ::NTuple{3,Number})
    # acoustical reflaction coefficient and counts
    (T₆₀>3.00) ? (T₆₀=3.00) : (T₆₀=T₆₀) # max time
    (T₆₀<0.05) ? (T₆₀=0.05) : (T₆₀=T₆₀) # min time
    x,y,z = Lᵣ
    V = x * y * z                 # volume
    S = (x*y + x*z + y*z)*2       # area
    β  = 0.1611                   # 24log(10)/344
    m₄ = 0.0104                   # 4m value
    L  = 4 * V / S                # mean free propagatation length
    α = 1 - exp((m₄ - β/T₆₀)*V/S) # acoustical absorption coefficient
    C = floor(Int, 344 * T₆₀ / L) # reflaction counts
    γ = sqrt(1-α)                 # reflaction coefficient
    return (-γ, C÷6)
end

"""
    y = directconv(x, h)
convolves vectors x and h. The resulting vector is length length(x)+length(h)-1 which
has a O(Lx*Lh) complexity.
"""
function directconv(x, h)
    Lx = veclen(x)
    Lh = veclen(h)
    L  = Lx + Lh - 1
    y  = zeros(eltype(x), L)
    c  = h[end:-1:1]
    # boundary left y[1:Lh-1]
    for i=1:Lh-1
        y[i] = sum(x[1:i] .* h[i:-1:1])
    end
    # middle part y[Lh:Lx]
    for i=Lh:Lx
        y[i] = sum(x[i-Lh+1:i] .* c)
    end
    # right part y[Lx+1:Ly]
    for i=0:Lh-2
        y[L-i] = sum(x[Lx-i:Lx] .* h[Lh:-1:Lh-i])
    end
    return y
end

"""
    y = conv(x, h)
convolves signals x and h. The resulting signal has length (Lx + Lh - 1) which
has nearly O(Lx*Log(Lx)) complexity when Lx>>Lh. This function uses fft to reduce complexity.
"""
function conv(x, h)
    Lx = veclen(x)
    Lh = veclen(h)
    L  = Lx + Lh - 1
    xₑ = zeros(eltype(x), L)
    hₑ = zeros(eltype(h), L)
    xₑ[1:Lx] = x
    hₑ[1:Lh] = h
    return real(ifft( fft(xₑ) .* fft(hₑ) ))
end

"""
    y = rir(fs::Number, T60::Number, room::NTuple{3,Number}, src::NTuple{3,Number}, mic::NTuple{3,Number})
Generate room impulse response.
"""
function rir(fs::Number, T60::Number, room::NTuple{3,Number}, src::NTuple{3,Number}, mic::NTuple{3,Number})
    γ, C = arcrc(T60, room)
    N = 2*C + 1
    L = floor(Int, fs*T60)
    T = eltype(T60)
    x = zeros(T, N)
    y = zeros(T, N)
    z = zeros(T, N)
    h = zeros(T, L)
    xr,yr,zr = room  # room's size in meters
    xs,ys,zs = src   # source's coordinates in meters
    xm,ym,zm = mic   # micphone's coordinates in meters
    for i=-C:1:C
        if i&1==1 # odd number
            x[i+C+1] = xs + i*xr - xm;
            y[i+C+1] = ys + i*yr - ym;
            z[i+C+1] = zs + i*zr - zm;
        else # even number
            x[i+C+1] = -xs + (i+1)*xr - xm;
            y[i+C+1] = -ys + (i+1)*yr - ym;
            z[i+C+1] = -zs + (i+1)*zr - zm;
        end
    end
    maxi = 0
    for i=-C:1:C
        for j=-C:1:C
            for k=-C:1:C
                d = sqrt( x[i+C+1]^2 + y[j+C+1]^2 + z[k+C+1]^2 )
                n = floor(Int, d / 344.0 * fs)
                if n>L; continue; end
                if maxi<n; maxi=n;end
                m = abs(i) + abs(j) + abs(k)
                h[n] += γ^m / d
            end
        end
    end
    return h[1:maxi]
end

"""
    y = addEcho(wav::Array, fs::Number, T60::Number,
                room::NTuple{3,Number},
                src::NTuple{3,Number},
                mic::NTuple{3,Number})
Generate impulse response function online, and add reverberation effect to wav.
## Arguments
- `wav`: sound samples
- `fs`: sound's sampling rate
- `T60`: effective reverberation time
- `room`: room's size in meters, e.g., (3, 4, 2.5)
- `src`: source's coordinates in meters, e.g., (1, 2, 1.2)
- `mic`: micphone's coordinates in meters, e.g., (1.6, 2.6, 1.0)
"""
function addEcho(wav::Array, fs::Number, T60::Number,
                 room::NTuple{3,Number},
                 src::NTuple{3,Number},
                 mic::NTuple{3,Number})
    return conv(wav, rir(fs, T60, room, src, mic))
end

function initAddEcho(path::String, period::Int)
    counter = 1
    h = nothing
    FILES = readdir(path)
    function addecho(speech::Array)
        if counter == 1
            file = rand(FILES,1)[1]
            @assert endswith(file, "wav") "$path should only keep *.wav files"
            h, fs = wavread(joinpath(path, file))
        end # read another Room Impulse Response every period
        (counter == period) ? (counter=1) : (counter+=1)
        return conv(speech, h)
    end
    return addecho
end

"""
    initAddEcho(fs::Number, T₆₀Span::NTuple{2,Number}, roomSpan::NTuple{6,Number}) -> addecho(wav::Array)
init reverberation effect function.
+ `fs` sampling rate
+ `T₆₀Span` T₆₀ range e.g. (0.05, 0.5)
+ `roomSpan`room size e.g. (MinL, MaxL, MinW, MaxW, MinH, MaxH)
"""
function initAddEcho(fs::Number, T₆₀Span::NTuple{2,Number}, roomSpan::NTuple{6,Number})
    MinT₆₀, MaxT₆₀ = T₆₀Span
    MinL, MaxL, MinW, MaxW, MinH, MaxH = roomSpan
    @assert MinT₆₀ <= MaxT₆₀
    @assert MinL <= MaxL
    @assert MinW <= MaxW
    @assert MinH <= MaxH
    function addecho(wav::Array)
        T₆₀ = rand()*(MaxT₆₀ - MinT₆₀) + MinT₆₀
        Lx  = rand()*(MaxL - MinL) + MinL
        Ly  = rand()*(MaxW - MinW) + MinW
        Lz  = rand()*(MaxH - MinH) + MinH
        room = (Lx, Ly, Lz)
        src  = (rand()*Lx, rand()*Ly, rand()*Lz)
        mic  = (rand()*Lx, rand()*Ly, rand()*Lz)
        return conv(wav, rir(fs, T₆₀, room, src, mic))
    end
    return addecho
end

max_stars_repo_name: sonosole/SpeechAugment

max_issues_repo_name: sonosole/SpeechAugment

max_forks_repo_name: sonosole/SpeechAugment


text:
mutable struct IgHolding
    symbol::String
    shares::Float64
    date::Date
    purchase_price::Float64
end

mutable struct IgPortfolio
    cash::Float64
    holdings::Vector{IgHolding}
end

struct IgUserError
    message::String
end

struct IgSystemError
    message::String
end

abstract type AbstractPortfolioOutputView end
struct PrettyView <: AbstractPortfolioOutputView end
struct SimpleView <: AbstractPortfolioOutputView end

# JSON3 struct bindings
StructTypes.StructType(::Type{IgHolding}) = StructTypes.Struct()
StructTypes.StructType(::Type{IgPortfolio}) = StructTypes.Struct()

ext: jl


text:
using TuringMachine

input_file = "example_input_1.txt"
program_file = "example_program_1.txt"

input = load_input(input_file)
program, init_state, halt_state = load_program(program_file)

state, tape_left, tape_right = set_up(init_state, input)

# run simulation step-by-step until we reach a halting state of 'accept' or 'reject'
while true
    try
        global state, tape_left, tape_right = simulate(state, program, tape_left, tape_right)
        println(state)
        state ==("qAccept") && break
        state ==("qReject") && break
    catch err
        println(err) && break
    end
end


text:
longest(a::String, b::String) = length(a) ≥ length(b) ? a : b

"""
julia> lcsrecursive("thisisatest", "testing123testing")
"tsitest"
"""
# Recursive
function lcsrecursive(xstr::String, ystr::String)
    if length(xstr) == 0 || length(ystr) == 0
        return ""
    end

    x, xs, y, ys = xstr[1], xstr[2:end], ystr[1], ystr[2:end]
    if x == y
        return string(x, lcsrecursive(xs, ys))
    else
        return longest(lcsrecursive(xstr, ys), lcsrecursive(xs, ystr))
    end
end

# Dynamic
function lcsdynamic(a::String, b::String)
    lengths = zeros(Int, length(a) + 1, length(b) + 1)

    # row 0 and column 0 are initialized to 0 already
    for (i, x) in enumerate(a), (j, y) in enumerate(b)
        if x == y
            lengths[i+1, j+1] = lengths[i, j] + 1
        else
            lengths[i+1, j+1] = max(lengths[i+1, j], lengths[i, j+1])
        end
    end

    # read the substring out from the matrix
    result = ""
    x, y = length(a) + 1, length(b) + 1
    while x > 1 && y > 1
        if lengths[x, y] == lengths[x-1, y]
            x -= 1
        elseif lengths[x, y] == lengths[x, y-1]
            y -= 1
        else
            @assert a[x-1] == b[y-1]
            result = string(a[x-1], result)
            x -= 1
            y -= 1
        end
    end

    return result
end

@show lcsrecursive("thisisatest", "testing123testing")
@time lcsrecursive("thisisatest", "testing123testing")
@show lcsdynamic("thisisatest", "testing123testing")
@time lcsdynamic("thisisatest", "testing123testing")

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
"""
    PolarQuadrature{N,T<:Real}

Holds information of the polar quadrature, such as the polar angles `θs` and the polar
weights `ωₚ`.
"""
struct PolarQuadrature{N,T<:Real}
    sinθs::Vector{T} # TODO: SVector as type!
    θs::Vector{T}
    ωₚ::Vector{T}
end

npolar(::PolarQuadrature{N}) where {N} = N
npolar2(::PolarQuadrature{N}) where {N} = div(N, 2)

for polar_quad in (:TabuchiYamamoto, :GaussLegendre, :EqualWeight, :EqualAngle, :Leonard)
    @eval begin
        $polar_quad(N::Int, T::Type{<:Real}=Float64) = $polar_quad(Val(N), T)
    end
end

const _TL_NΘ_ = (2, 4, 6)
const _GL_NΘ_ = (2, 3, 4, 6, 8, 10, 12)

function TabuchiYamamoto(n_polar::Val{N}, T::Type{<:Real}=Float64) where {N}
    N in _TL_NΘ_ || throw(DomainError(N, "number of polar angles for TabuchiYamamoto must" *
                                         " be in $(string(_TL_NΘ_))."))

    sinθs, θs, ωₚ = ntuple(_ -> Vector{T}(undef, N), 3)

    # Nuclear Engineering Handbook, pg. 1154
    set_tabuchiyamamoto_data!(n_polar, sinθs, ωₚ)

    # compute remaining data from loaded data
    for i in 1:div(N, 2)
        θs[i] = asin(sinθs[i])

        # suplementaries
        j = N - i + 1
        sinθs[j] = sinθs[i]
        θs[j] = π - θs[i]
        ωₚ[j] = ωₚ[i]
    end

    return PolarQuadrature{N,T}(sinθs, θs, ωₚ)
end

function set_tabuchiyamamoto_data!(::Val{2}, sinθs, ωₚ)
    sinθs[1] = 0.798184
    ωₚ[1] = 1 / 2
    return nothing
end

function set_tabuchiyamamoto_data!(::Val{4}, sinθs, ωₚ)
    sinθs[1] = 0.363900
    sinθs[2] = 0.899900
    ωₚ[1] = 0.212854 / 2
    ωₚ[2] = 0.787146 / 2
    return nothing
end

function set_tabuchiyamamoto_data!(::Val{6}, sinθs, ωₚ)
    sinθs[1] = 0.166648
    sinθs[2] = 0.537707
    sinθs[3] = 0.932954
    ωₚ[1] = 0.046233 / 2
    ωₚ[2] = 0.283619 / 2
    ωₚ[3] = 0.670148 / 2
    return nothing
end

function GaussLegendre(n_polar::Val{N}, T::Type{<:Real}=Float64) where {N}
    N in _GL_NΘ_ || throw(DomainError(N, "number of polar angles for GaussLegendre must be in $(string(_GL_NΘ_))."))

    sinθs, θs, ωₚ = ntuple(_ -> Vector{T}(undef, N), 3)

    # Nuclear Engineering Handbook, pg. 1153
    set_gausslegendre_data!(n_polar, θs, ωₚ)

    # compute remaining data from loaded data
    for i in 1:div(N, 2)
        sinθs[i] = sin(θs[i])

        # suplementaries
        j = N - i + 1
        sinθs[j] = sinθs[i]
        θs[j] = π - θs[i]
        ωₚ[j] = ωₚ[i]
    end

    return PolarQuadrature{N,T}(sinθs, θs, ωₚ)
end

function set_gausslegendre_data!(::Val{2}, θs, ωₚ)
    θs[1] = acos(0.5773502691)
    ωₚ[1] = 1 / 2
    return nothing
end

function set_gausslegendre_data!(::Val{4}, θs, ωₚ)
    θs[1] = acos(0.3399810435)
    θs[2] = acos(0.8611363115)
    ωₚ[1] = 0.6521451549 / 2
    ωₚ[2] = 0.3478548451 / 2
    return nothing
end

function set_gausslegendre_data!(::Val{6}, θs, ωₚ)
    θs[1] = acos(0.2386191860)
    θs[2] = acos(0.6612093864)
    θs[3] = acos(0.9324695142)
    ωₚ[1] = 0.4679139346 / 2
    ωₚ[2] = 0.3607615730 / 2
    ωₚ[3] = 0.1713244924 / 2
    return nothing
end

function set_gausslegendre_data!(::Val{8}, θs, ωₚ)
    θs[1] = acos(0.1834346424)
    θs[2] = acos(0.5255324099)
    θs[3] = acos(0.7966664774)
    θs[4] = acos(0.9602898564)
    ωₚ[1] = 0.3626837834 / 2
    ωₚ[2] = 0.3137066459 / 2
    ωₚ[3] = 0.2223810344 / 2
    ωₚ[4] = 0.1012285363 / 2
    return nothing
end

function set_gausslegendre_data!(::Val{10}, θs, ωₚ)
    θs[1] = acos(0.1488743387)
    θs[2] = acos(0.4333953941)
    θs[3] = acos(0.6794095682)
    θs[4] = acos(0.8650633666)
    θs[5] = acos(0.9739065285)
    ωₚ[1] = 0.2955242247 / 2
    ωₚ[2] = 0.2692667193 / 2
    ωₚ[3] = 0.2190863625 / 2
    ωₚ[4] = 0.1494513492 / 2
    ωₚ[5] = 0.0666713443 / 2
    return nothing
end

function set_gausslegendre_data!(::Val{12}, θs, ωₚ)
    θs[1] = acos(0.1252334085)
    θs[2] = acos(0.3678314989)
    θs[3] = acos(0.5873179542)
    θs[4] = acos(0.7699026741)
    θs[5] = acos(0.9041172563)
    θs[6] = acos(0.9815606342)
    ωₚ[1] = 0.2491470458 / 2
    ωₚ[2] = 0.2334925365 / 2
    ωₚ[3] = 0.2031674267 / 2
    ωₚ[4] = 0.1600783286 / 2
    ωₚ[5] = 0.1069393260 / 2
    ωₚ[6] = 0.0471753364 / 2
    return nothing
end

function Leonard(n_polar::Val{N}, T::Type{<:Real}=Float64) where {N}
    N in _TL_NΘ_ || throw(DomainError(N, "number of polar angles for Leonard must be in $(string(_TL_NΘ_))."))

    sinθs, θs, ωₚ = ntuple(_ -> Vector{T}(undef, N), 3)

    # TODO: decide which source we should use since Tabuchi paper and OpenMoC have differ.
    # On the other hand, the Nuclear Engineering Handbook has Herbert optimized values.
    # The following values are from Tabuchi paper
    set_leonard_data!(n_polar, θs, ωₚ)

    # compute remaining data from loaded data
    for i in 1:div(N, 2)
        θs[i] = asin(sinθs[i])

        # suplementaries
        j = N - i + 1
        sinθs[j] = sinθs[i]
        θs[j] = π - θs[i]
        ωₚ[j] = ωₚ[i]
    end

    return PolarQuadrature{N,T}(sinθs, θs, ωₚ)
end

function set_leonard_data!(::Val{2}, sinθs, ωₚ)
    sinθs[1] = 0.752244
    ωₚ[1] = 1 / 2
    return nothing
end

function set_leonard_data!(::Val{4}, sinθs, ωₚ)
    sinθs[1] = 0.273658
    sinθs[2] = 0.865714
    ωₚ[1] = 0.139473 / 2
    ωₚ[2] = 0.860527 / 2
    return nothing
end

function set_leonard_data!(::Val{6}, sinθs, ωₚ)
    sinθs[1] = 0.103840
    sinθs[2] = 0.430723
    sinθs[3] = 0.905435
    ωₚ[1] = 0.020530 / 2
    ωₚ[2] = 0.219161 / 2
    ωₚ[3] = 0.760309 / 2
    return nothing
end


text:
using ShaleDrillingLikelihood: zero_out_small_probs!

import ShaleDrillingLikelihood: logprice, logrigrate

using Statistics: var, cov

export GridTransition

# ----------------------------
# Make Grid + transition
# ----------------------------

mysteprange(a,b,n) = range(a,b;length=n)

function PriceGrid(z, delta, len)
    x = first.(z)
    prng = range(minimum(x)-delta, maximum(x)+delta; length=len)
    sigsq = var(diff(x))
    println("Var(logprice) = $sigsq and sd(logprice) = $(sqrt(sigsq))")
    P = tauchen_1d(prng, identity, sigsq)
    return prng, P
end

function PriceCostGrid(z, delta, len)
    p = logprice.(z)
    c = logrigrate.(z)
    pc = hcat(p,c)

    # ranges
    mins = vec(minimum(pc; dims=1))
    maxs = vec(maximum(pc; dims=1))
    a = mins .- delta
    b = maxs .+ delta
    pcrng = mysteprange.(a, b, len)
    pcprod = product(pcrng...)

    # sigma
    dpc = diff(pc; dims=1)

    Sigma = cov(dpc)
    stdvs = sqrt.(diag(Sigma))
    sigcor = Sigma[2,1] / prod(stdvs)
    println("Sigma(logprice, logrigrate) = $Sigma.\nStdv(logprice, logrigrate) = $(round.(stdvs;digits=5)) and correlation = $(round(sigcor;digits=4))")

    P = tauchen_2d(pcprod, identity, Sigma)
    return pcrng, P
end

function YearGrid(z, args...)
    y = ShaleDrillingLikelihood.year.(z)
    miny, maxy = extrema(y)
    yrng = UnitRange(miny, maxy)
    nyears = length(yrng)
    nyearsm1 = nyears-1
    @assert nyearsm1 == maxy-miny

    l = zeros(nyearsm1)
    m = vcat(fill(3/4,nyearsm1), 1.0)
    r = fill(1/4, nyearsm1)

    P = Tridiagonal(l,m,r)
    return yrng, P
end

function GridTransition(rwrd::DrillReward{R,C}, z::Vector, delta::Number, len::Number; minp=minp_default()) where {
        R, C<:Union{DrillingCost_constant, DrillingCost_dgt1}
    }
    rng, P = PriceGrid(z,delta,len)
    zero_out_small_probs!(P, minp)
    return tuple(rng), sparse(P)
end

# # rigrate + price
# function GridTransition(rwrd::DrillReward{R,C}, z, delta, len; minp=minp_default()) where {
#         R, C<:Union{DrillingCost_constant, DrillingCost_dgt1}
#     }
#     rng, P = PriceCostGrid(z, delta, len)
#     zero_out_small_probs!(P, minp)
#     return rng, sparse(P)
# end

function GridTransition(rwrd::DrillReward{R,C}, z::Vector, delta::Number, len::Number; minp=minp_default()) where {
        R, C<:AbstractDrillingCost_TimeFE
    }
    p = tuple.(logprice.(z))
    prng, pP = PriceGrid(p, delta, len)
    yrng, yP = YearGrid(z)
    rng = (prng, yrng)
    P = kron(yP, pP)
    zero_out_small_probs!(P, minp)
    return rng, sparse(P)
end

function GridTransition(rwrd::DrillReward{R,C}, z::Vector, delta::Number, len::Number; minp=minp_default()) where {
        R, C<:Union{DrillingCost_TimeFE_rigrate, DrillingCost_TimeFE_rig_costdiffs}
    }
    pcrng, pcP = PriceCostGrid(z,delta,len)
    yrng, yP = YearGrid(z)
    rng = (pcrng..., yrng)
    P = kron(yP, pcP)
    zero_out_small_probs!(P, minp)
    return rng, sparse(P)
end

GridTransition(rwrd::DrillReward, etv::ExogTimeVars, args...; kwargs...) = GridTransition(rwrd, _timevars(etv), args...; kwargs...)
GridTransition(d::DataDrillPrimitive, args...; kwargs...) = GridTransition(reward(d), zchars(d), args...; kwargs...)
GridTransition(d::DataDrill,          args...; kwargs...) = zspace(_model(d)), ztransition(_model(d))


text:
@testset "callback" begin
    # Exploitability Callback
    game = CFR.Kuhn()
    sol = CFRSolver(game)
    cb = CFR.ExploitabilityCallback(sol, 100)
    train!(sol, 100_000, cb = cb)

    @test RecipesBase.apply_recipe(Dict{Symbol,Any}(), cb) ≠ nothing
    @test length(cb.hist.y) == length(cb.hist.x) == 1_000
    @test 0.0 < last(cb.hist.y) < 1e-2

    game = CFR.IIEMatrixGame([(randn(), randn()) for i in 1:5, j in 1:5])
    sol = CFRSolver(game)
    cb = CFR.ExploitabilityCallback(sol, 100)
    train!(sol, 100_000, cb = cb)

    @test RecipesBase.apply_recipe(Dict{Symbol,Any}(), cb) ≠ nothing
    @test RecipesBase.apply_recipe(Dict{Symbol,Any}(), cb.hist) ≠ nothing
    @test 0.0 < last(cb.hist.y) < 1e-2
    @test length(cb.hist.y) == length(cb.hist.x) == 1_000

    ## Throttle
    sol = ESCFRSolver(game)
    io = IOBuffer()
    f = () -> println(io, length(sol.I))
    save_freq = 10
    train_iter = 101
    cb = CFR.Throttle(f,save_freq)
    train!(sol, train_iter, cb = cb)
    str = String(take!(io))
    @test length(split(str, "\n")) == div(train_iter, save_freq) + 2

    ## chain
    sol = ESCFRSolver(game)
    io = IOBuffer()
    g = () -> println(io, length(sol.I))
    save_freq = 10
    train_iter = 1001
    ecb = CFR.ExploitabilityCallback(sol, 100)
    cb = CFR.CallbackChain(ecb,  CFR.Throttle(g, save_freq))
    train!(sol, train_iter, cb = cb)

    str = String(take!(io))
    @test first(ecb.hist.y) > last(ecb.hist.y) > 0.0
    @test length(ecb.hist.y) == length(ecb.hist.x) == 11
    @test length(split(str, "\n")) == div(train_iter, save_freq) + 2
end

max_stars_repo_name: WhiffleFish/HelloCFR

max_issues_repo_name: WhiffleFish/HelloCFR

max_forks_repo_name: WhiffleFish/HelloCFR


text:
using VoronoiCells
using GeometryBasics
import Random
using Test

@testset "Rectangles" begin
    @testset "Construct Rectangles" begin
        @testset "Construct with with GeometryBasics integer points" begin
            rect = Rectangle(Point2(0, 0), Point2(1, 1))
            @test isa(rect, Rectangle{GeometryBasics.Point{2, Float64}})
        end

        @testset "Construct with with GeometryBasics floating points" begin
            rect = Rectangle(Point2(0.0, 0.0), Point2(1.0, 1.0))
            @test isa(rect, Rectangle{Point{2, Float64}})
        end

        @testset "Construct with VoronoiDelaunay points" begin
            rect = Rectangle(VoronoiDelaunay.Point2D(1.1, 1.1), VoronoiDelaunay.Point2D(1.6, 1.7))
            @test isa(rect, Rectangle{VoronoiDelaunay.Point2D})
        end

        @testset "Invalid Rectangles" begin
            p1 = GeometryBasics.Point2(0.0, 0.0)
            p2 = GeometryBasics.Point2(0.0, 1.0)

            @test_throws ArgumentError Rectangle(p1, p2)
        end
    end

    @testset "Is point in rectangle?" begin
        @testset "Point is in rectangle" begin
            point = GeometryBasics.Point2(1, 1)
            rect = Rectangle(GeometryBasics.Point2(0, 0), GeometryBasics.Point2(2, 2))

            @test VoronoiCells.isinside(point, rect)
        end

        @testset "Point is not in rectangle" begin
            point = GeometryBasics.Point2(-1, 1)
            rect = Rectangle(GeometryBasics.Point2(0, 0), GeometryBasics.Point2(2, 2))

            @test !VoronoiCells.isinside(point, rect)
        end
    end

    @testset "Map points between recangles" begin
        Random.seed!(1337)
        points = [GeometryBasics.Point2(rand(), rand()) for _ in 1:3]

        from_rect = Rectangle(GeometryBasics.Point2(0, 0), GeometryBasics.Point2(1, 1))
        to_rect = Rectangle(GeometryBasics.Point2(0, 0), GeometryBasics.Point2(2, 2))

        transformed_points = VoronoiCells.map_rectangle(points, from_rect, to_rect)
        double_transformed_points = VoronoiCells.map_rectangle(transformed_points, to_rect, from_rect)

        @test all(points .== double_transformed_points)
    end

    @testset "Error mapping points outside rectangle" begin
        points = [GeometryBasics.Point2(-1.0, 0.0)]

        from_rect = Rectangle(GeometryBasics.Point2(0, 0), GeometryBasics.Point2(1, 1))
        to_rect = Rectangle(GeometryBasics.Point2(0, 0), GeometryBasics.Point2(2, 2))

        @test_throws ErrorException VoronoiCells.map_rectangle(points, from_rect, to_rect)
    end

    @testset "Find points nearest to each rectangle corner" begin
        rect = Rectangle(GeometryBasics.Point2(0, 0), GeometryBasics.Point2(1, 1))

        @testset "Unique nearest neighbors" begin
            rng = Random.MersenneTwister(1)
            points = [Point2(rand(rng), rand(rng)) for _ in 1:5]

            corner_neighbors = VoronoiCells.corner_nearest_neighbor(points, rect)

            expected_neighbors = Dict(
                1 => [4], 2 => [5], 3 => [2], 4 => [3]
            )

            @test corner_neighbors == expected_neighbors
        end

        @testset "Multiple corners with the same nearest neighbors" begin
            points = [Point2(0.75, 0.5), Point2(0.25, 0.75), Point2(0.25, 0.25)]
            corner_neighbors = VoronoiCells.corner_nearest_neighbor(points, rect)

            expected_neighbors = Dict(
                1 => [1], 2 => [2], 3 => [3], 4 => [1]
            )

            @test corner_neighbors == expected_neighbors
        end

        @testset "Corners with multiple nearest neighbors" begin
            points = [Point2(0.75, 0.5), Point2(0.5, 0.75), Point2(0.5, 0.5)]
            corner_neighbors = VoronoiCells.corner_nearest_neighbor(points, rect)

            expected_neighbors = Dict(
                1 => [1, 2], 2 => [2], 3 => [3], 4 => [1]
            )

            @test corner_neighbors == expected_neighbors
        end
    end
end


text:
module HalfIntegers

#   This file is part of HalfIntegers.jl. It is licensed under the MIT license
#   HalfIntegers Copyright (C) 2021 Michael Reed

export Half, twice

struct Half{T} <: Real
    n::T
    Half(n::T) where T<:Integer = new{T}(n)
    Half{T}(n::T) where T<:Integer = new{T}(n)
end

Base.show(io::IO,n::Half) = print(io,Float64(n))

# access / conversion

twice(n::Half) = n.n
Base.Float64(n::Half) = twice(n)/2
Base.convert(::Type{Half{T}},n::T) where T<:Integer = Half{T}(2n)
Base.convert(::Type{Half},n::T) where T<:Integer = Half{T}(2n)

# algebra

for op ∈ (:+,:-,:*,:/)
    @eval begin
        Base.$op(a::Half,b::Half) = Half($op(twice(a),twice(b)))
        Base.$op(a::Half,b::T) where T<:Integer = Half($op(twice(a),2b))
        Base.$op(a::T,b::Half) where T<:Integer = Half($op(2a,twice(b)))
        Base.$op(a::Half,b::T) where T<:Number = $op(Float64(a),b)
        Base.$op(a::T,b::Half) where T<:Number = $op(a,Float64(b))
    end
end

for op ∈ (:(==),:isapprox)
    @eval begin
        Base.$op(a::Half,b::Half) = $op(twice(a),twice(b))
        Base.$op(a::Half,b::T) where T<:Integer = $op(twice(a),2b)
        Base.$op(a::T,b::Half) where T<:Integer = $op(2a,twice(b))
        Base.$op(a::Half,b::T) where T<:Number = $op(Float64(a),b)
        Base.$op(a::T,b::Half) where T<:Number = $op(a,Float64(b))
    end
end

end # module


text:
unary_operators = [:(+), :(-), :(!)]

numeric_unary_operators = [:(+), :(-)]

logical_unary_operators = [:(!)]

elementary_functions = [:abs, :sign, :acos, :acosh, :asin,
                        :asinh, :atan, :atanh, :sin, :sinh,
                        :cos, :cosh, :tan, :tanh, :ceil, :floor,
                        :round, :trunc, :exp, :exp2, :expm1, :log, :log10, :log1p,
                        :log2, :logb, :sqrt, :gamma, :lgamma, :digamma,
                        :erf, :erfc, :square]

two_argument_elementary_functions = [:round, :ceil, :floor, :trunc]

special_comparison_operators = [:isless]

comparison_operators = [:(==), :(.==), :(!=), :(.!=),
                        :(>), :(.>), :(>=), :(.>=), :(<), :(.<),
                        :(<=), :(.<=)]

scalar_comparison_operators = [:(==), :(!=), :(>), :(>=),
                               :(<), :(<=)]

array_comparison_operators = [:(.==), :(.!=), :(.>), :(.>=), :(.<), :(.<=)]

vectorized_comparison_operators = [(:(.==), :(==)), (:(.!=), :(!=)),
                                   (:(.>), :(>)), (:(.>=), :(>=)),
                                   (:(.<), :(<)), (:(.<=), :(<=))]

binary_operators = [:(+), :(.+), :(-), :(.-), :(*), :(.*), :(/), :(./),
                    :(.^), :(div), :(mod), :(fld), :(rem)]

induced_binary_operators = [:(^)]

arithmetic_operators = [:(+), :(.+), :(-), :(.-), :(*), :(.*), :(/), :(./),
                        :(.^), :(div), :(mod), :(fld), :(rem)]

induced_arithmetic_operators = [:(^)]

biscalar_operators = [:(max), :(min)]

scalar_arithmetic_operators = [:(+), :(-), :(*), :(/),
                               :(div), :(mod), :(fld), :(rem)]

induced_scalar_arithmetic_operators = [:(^)]

array_arithmetic_operators = [:(+), :(.+), :(-), :(.-), :(.*), :(./), :(.^)]

bit_operators = [:(&), :(|), :($)]

unary_vector_operators = [:min, :max, :prod, :sum, :mean, :median, :std,
                          :var, :mad, :norm, :skewness, :kurtosis]

# TODO: dist, iqr, rle, inverse_rle

pairwise_vector_operators = [:diff, :reldiff, :percent_change]

cumulative_vector_operators = [:cumprod, :cumsum, :cumsum_kbn, :cummin, :cummax]

ffts = [:fft]

binary_vector_operators = [:dot, :cor_pearson, :cov_pearson,
                           :cor_spearman, :cov_spearman]

rowwise_operators = [:rowmins, :rowmaxs, :rowprods, :rowsums,
                     :rowmeans, :rowmedians, :rowstds, :rowvars,
                     :rowffts, :rownorms]

columnar_operators = [:colmins, :colmaxs, :colprods, :colsums,
                      :colmeans, :colmedians, :colstds, :colvars,
                      :colffts, :colnorms]

boolean_operators = [:any, :all]

for f in unary_operators
    @eval begin
        function ($f)(d::NAtype)
            return NA
        end
        function ($f){T}(dv::DataVector{T})
            res = deepcopy(dv)
            for i in 1:length(dv)
                res[i] = ($f)(dv[i])
            end
            return res
        end
        function ($f){T}(dm::DataMatrix{T})
            res = deepcopy(dm)
            for i in 1:length(dm)
                res[i] = ($f)(dm[i])
            end
            return res
        end
        function ($f)(df::DataFrame)
            res = deepcopy(df)
            n, p = nrow(df), ncol(df)
            for j in 1:p
                if typeof(df[j]).parameters[1] <: Number
                    for i in 1:n
                        res[i, j] = ($f)(df[i, j])
                    end
                else
                    for i in 1:n
                        res[i, j] = NA
                    end
                end
            end
            return res
        end
    end
end

# Treat ctranspose and * in a special way for now
function ctranspose(d::DataArray)
    return DataArray(d.data', d.na')
end

# TODO: Check there are no better algorithms
function (*){S <: Real, T <: Real}(a::DataVector{S}, b::DataMatrix{T})
    if size(b, 1) != 1
        error("DataVector and matrix sizes must match")
    end
    n, p = length(a), size(b, 2)
    res = datazeros(n, p)
    for i in 1:n
        for j in 1:p
            res[i, j] = a[i] * b[j]
        end
    end
    return res
end

function (*){S <: Real, T <: Real}(a::Vector{S}, b::DataMatrix{T})
    if size(b, 1) != 1
        error("Vector and matrix sizes must match")
    end
    n, p = length(a), size(b, 2)
    res = datazeros(n, p)
    for i in 1:n
        for j in 1:p
            res[i, j] = a[i] * b[j]
        end
    end
    return res
end

# TODO: Check there are no better algorithms
function (*){S <: Real, T <: Real}(a::DataMatrix{S}, b::DataVector{T})
    if size(a, 2) != length(b)
        error("The number of columns of the DataMatrix must match the length of the DataVector")
    end
    n, p = size(a, 1), length(b)
    res = datazeros(n)
    for i in 1:n
        res[i] = 0.0
        for j in 1:p
            res[i] += a[i, j] * b[j]
        end
    end
    return res
end

# TODO: Check there are no better algorithms
function (*){S <: Real, T <: Real}(a::DataMatrix{S}, b::Vector{T})
    if size(a, 2) != length(b)
        error("The number of columns of the DataMatrix must match the length of the Vector")
    end
    n, p = size(a, 1), length(b)
    res = datazeros(n)
    for i in 1:n
        res[i] = 0.0
        for j in 1:p
            res[i] += a[i, j] * b[j]
        end
    end
    return res
end

# Propagates NA's
# For a dissenting view,
# http://radfordneal.wordpress.com/2011/05/21/slowing-down-matrix-multiplication-in-r/
# But we're getting 10x R while maintaining NA's
function (*){S <: Real, T <: Real}(a::DataMatrix{S}, b::DataMatrix{T})
    n1, p1 = size(a)
    n2, p2 = size(b)
    if p1 != n2
        error("DataMatrix sizes must align for matrix multiplication")
    end
    res = DataArray(a.data * b.data, falses(n1, p2))
    # Propagation can be made more efficient by storing record of corrupt
    # rows and columns, then doing fast edits.
    corrupt_rows = falses(n1)
    corrupt_cols = falses(p2)
    for i in 1:n1
        for j in 1:p1
            if a.na[i, j]
                # Propagate NA's
                # Corrupt all rows based on i
                corrupt_rows[i] = true
            end
        end
    end
    for i in 1:n2
        for j in 1:p2
            if b.na[i, j]
                # Propagate NA's
                # Corrupt all columns based on j
                corrupt_cols[j] = true
            end
        end
    end
    for i in 1:n1
        if corrupt_rows[i]
            res.na[i, :] = true
        end
    end
    for j in 1:p2
        if corrupt_cols[j]
            res.na[:, j] = true
        end
    end
    return res
end

function (*){S <: Real, T <: Real}(a::DataMatrix{S}, b::Matrix{T})
    n1, p1 = size(a)
    n2, p2 = size(b)
    if p1 != n2
        error("DataMatrix and Matrix sizes must align for matrix multiplication")
    end
    res = DataArray(a.data * b, falses(n1, p2))
    for i in 1:n1
        for j in 1:p1
            if a.na[i, j]
                # Propagate NA's
                # Corrupt all rows based on i
                res.na[i, :] = true
            end
        end
    end
    return res
end

function (*){S <: Real, T <: Real}(a::Matrix{S}, b::DataMatrix{T})
    n1, p1 = size(a)
    n2, p2 = size(b)
    if p1 != n2
        error("Matrix and DataMatrix sizes must align for matrix multiplication")
    end
    res = DataArray(a * b.data, falses(n1, p2))
    for i in 1:n2
        for j in 1:p2
            if b.na[i, j]
                # Propagate NA's
                # Corrupt all columns based on j
                res.na[:, j] = true
            end
        end
    end
    return res
end

for f in elementary_functions
    @eval begin
        function ($f)(d::NAtype)
            return NA
        end
        function ($f){T}(dv::DataVector{T})
            n = length(dv)
            res = DataArray(Array(T, n), falses(n))
            for i = 1:n
                res[i] = ($f)(dv[i])
            end
            return res
        end
        function ($f){T}(adv::AbstractDataVector{T})
            res = deepcopy(adv)
            for i = 1:length(adv)
                if isna(adv[i])
                    res[i] = NA
                else
                    res[i] = ($f)(adv[i])
                end
            end
            return res
        end
        function ($f){T}(dm::DataMatrix{T})
            res = DataArray(Array(T, size(dm)), falses(size(dm)))
            for i = 1:length(dm)
                res[i] = ($f)(dm[i])
            end
            return res
        end
        function ($f)(df::DataFrame)
            n, p = nrow(df), ncol(df)
            res = DataFrame(coltypes(df), colnames(df), n)
            for j in 1:p
                if eltype(df[j]) <: Number
                    for i in 1:n
                        res[i, j] = ($f)(df[i, j])
                    end
                else
                    for i in 1:n
                        res[i, j] = NA
                    end
                end
            end
            return res
        end
    end
end

for f in two_argument_elementary_functions
    @eval begin
        function ($f){T}(dv::DataVector{T}, arg2)
            n = length(dv)
            res = DataArray(Array(T, n), falses(n))
            for i = 1:n
                res[i] = ($f)(dv[i], arg2)
            end
            return res
        end
    end
end

for f in special_comparison_operators
    @eval begin
        function ($f)(d::NAtype, e::NAtype)
            return false
        end
        function ($f){T <: Union(String, Number)}(d::NAtype, x::T)
            return true
        end
        function ($f){T <: Union(String, Number)}(x::T, d::NAtype)
            return false
        end
    end
end

for f in comparison_operators
    @eval begin
        function ($f)(d::NAtype, e::NAtype)
            return NA
        end
        function ($f){T <: Union(String, Number)}(d::NAtype, x::T)
            return NA
        end
        function ($f){T <: Union(String, Number)}(x::T, d::NAtype)
            return NA
        end
    end
end

#
# Bit operators
#

(&)(a::NAtype, b::NAtype) = NA
(&)(a::NAtype, b::Bool) = b ? NA : false
(&)(a::Bool, b::NAtype) = a ? NA : false
(&)(a::NAtype, b::Real) = NA
(&)(a::Real, b::NAtype) = NA

(|)(a::NAtype, b::NAtype) = NA
(|)(a::NAtype, b::Bool) = b ? true : NA
(|)(a::Bool, b::NAtype) = a ? true : NA
(|)(a::NAtype, b::Real) = NA
(|)(a::Real, b::NAtype) = NA

($)(a::NAtype, b::NAtype) = NA
($)(a::NAtype, b::Bool) = NA
($)(a::Bool, b::NAtype) = NA
($)(a::NAtype, b::Real) = NA
($)(a::Real, b::NAtype) = NA

for f in bit_operators
    @eval begin
        function ($f)(a::Array, b::AbstractDataArray)
            res = similar(b, size(b))
            for i in 1:length(res)
                res[i] = ($f)(a[i], b[i])
            end
            return res
        end
        function ($f)(a::AbstractDataArray, b::Array)
            res = similar(a, size(a))
            for i in 1:length(res)
                res[i] = ($f)(a[i], b[i])
            end
            return res
        end
        function ($f)(a::AbstractDataArray, b::AbstractDataArray)
            res = similar(a, size(a))
            for i in 1:length(res)
                res[i] = ($f)(a[i], b[i])
            end
            return res
        end
    end
end

for (f, scalarf) in vectorized_comparison_operators
    @eval begin
        function ($f){S, T <: Union(String, Number)}(a::DataVector{S}, v::T)
            res = DataArray(Array(Bool, length(a)), BitArray(length(a)))
            for i in 1:length(a)
                if isna(a[i])
                    res[i] = NA
                else
                    res[i] = ($scalarf)(a[i], v)
                end
            end
            return res
        end
        function ($f){S <: Union(String, Number), T}(v::S, a::DataVector{T})
            res = DataArray(Array(Bool, length(a)), BitArray(length(a)))
            for i in 1:length(a)
                if isna(a[i])
                    res[i] = NA
                else
                    res[i] = ($scalarf)(v, a[i])
                end
            end
            return res
        end
        function ($f){S, T <: Union(String, Number)}(a::PooledDataVector{S}, v::T)
            res = PooledDataArray(Array(Bool, length(a)), BitArray(length(a)))
            for i in 1:length(a)
                if isna(a[i])
                    res[i] = NA
                else
                    res[i] = ($scalarf)(a[i], v)
                end
            end
            return res
        end
        function ($f){S <: Union(String, Number), T}(v::S, a::PooledDataVector{T})
            res = PooledDataArray(Array(Bool, length(a)), BitArray(length(a)))
            for i in 1:length(a)
                if isna(a[i])
                    res[i] = NA
                else
                    res[i] = ($scalarf)(v, a[i])
                end
            end
            return res
        end
        function ($f){T}(a::AbstractDataVector{T}, v::NAtype)
            res = DataArray(Array(Bool, length(a)), BitArray(length(a)))
            for i in 1:length(a)
                res[i] = NA
            end
            res
        end
        function ($f){T}(v::NAtype, a::AbstractDataVector{T})
            res = DataArray(Array(Bool, length(a)), BitArray(length(a)))
            for i in 1:length(a)
                res[i] = NA
            end
            res
        end
        function ($f){T}(a::PooledDataVector{T}, v::NAtype)
            res = PooledDataArray(Array(Bool, length(a)), BitArray(length(a)))
            for i in 1:length(a)
                res[i] = NA
            end
            res
        end
        function ($f){T}(v::NAtype, a::PooledDataVector{T})
            res = PooledDataArray(Array(Bool, length(a)), BitArray(length(a)))
            for i in 1:length(a)
                res[i] = NA
            end
            res
        end
        function ($f){T <: Number}(a::DataFrame, v::T)
            res = DataFrame(Array(Bool, size(a)))
            n, p = nrow(a), ncol(a)
            for j in 1:p
                if typeof(a[j]).parameters[1] <: Number
                    for i in 1:n
                        res[i, j] = isna(a[i, j]) ? NA : ($scalarf)(a[i, j], v)
                    end
                else
                    for i in 1:n
                        res[i, j] = NA
                    end
                end
            end
            return res
        end
        function ($f){T <: Number}(v::T, a::DataFrame)
            res = DataFrame(Array(Bool, size(a)))
            n, p = nrow(a), ncol(a)
            for j in 1:p
                if typeof(a[j]).parameters[1] <: Number
                    for i = 1:n
                        res[i, j] = isna(a[i, j]) ? NA : ($scalarf)(a[i, j], v)
                    end
                else
                    for i = 1:n
                        res[i, j] = NA
                    end
                end
            end
            return res
        end
        function ($f){T <: String}(a::DataFrame, v::T)
            res = DataFrame(Array(Bool, size(a)))
            n, p = nrow(a), ncol(a)
            for j in 1:p
                if typeof(a[j]).parameters[1] <: String
                    for i in 1:n
                        res[i, j] = isna(a[i, j]) ? NA : ($scalarf)(a[i, j], v)
                    end
                else
                    for i in 1:n
                        res[i, j] = NA
                    end
                end
            end
            return res
        end
        function ($f){T <: String}(v::T, a::DataFrame)
            res = DataFrame(Array(Bool, size(a)))
            n, p = nrow(a), ncol(a)
            for j in 1:p
                if typeof(a[j]).parameters[1] <: String
                    for i = 1:n
                        res[i, j] = isna(a[i, j]) ? NA : ($scalarf)(a[i, j], v)
                    end
                else
                    for i = 1:n
                        res[i, j] = NA
                    end
                end
            end
            return res
        end
        function ($f)(a::DataFrame, v::NAtype)
            res = DataFrame(Array(Bool, size(a)))
            n, p = nrow(a), ncol(a)
            for j in 1:p
                for i in 1:n
                    res[i, j] = NA
                end
            end
            return res
        end
        function ($f)(v::NAtype, a::DataFrame)
            res = DataFrame(Array(Bool, size(a)))
            n, p = nrow(a), ncol(a)
            for j in 1:p
                for i = 1:n
                    res[i, j] = NA
                end
            end
            return res
        end
    end
end

for f in scalar_comparison_operators
    @eval begin
        function ($f){S, T}(a::AbstractDataVector{S}, b::AbstractDataVector{T})
            error(string(string($f), " not defined for DataVectors. Try .", string($f)))
        end
        function ($f){S, T}(a::DataMatrix{S}, b::DataMatrix{T})
            error(string(string($f), " not defined for DataMatrix's. Try .", string($f)))
        end
        function ($f)(a::AbstractDataFrame, b::AbstractDataFrame)
            error(string(string($f), " not defined for DataFrames. Try .", string($f)))
        end
    end
end

for (f, scalarf) in vectorized_comparison_operators
    @eval begin
        function ($f){S, T}(a::AbstractDataVector{S}, b::AbstractDataVector{T})
            res = DataArray(Array(Bool, length(a)), BitArray(length(a)))
            for i in 1:length(a)
                if isna(a[i]) || isna(b[i])
                    res[i] = NA
                else
                    res[i] = ($scalarf)(a[i], b[i])
                end
            end
            return res
        end
        function ($f){S, T}(a::AbstractDataVector{S}, b::Vector{T})
            res = DataArray(Array(Bool, length(a)), BitArray(length(a)))
            for i in 1:length(a)
                if isna(a[i])
                    res[i] = NA
                else
                    res[i] = ($scalarf)(a[i], b[i])
                end
            end
            return res
        end
        function ($f){S, T}(a::Vector{S}, b::AbstractDataVector{T})
            res = DataArray(Array(Bool, length(a)), BitArray(length(a)))
            for i in 1:length(a)
                if isna(b[i])
                    res[i] = NA
                else
                    res[i] = ($scalarf)(a[i], b[i])
                end
            end
            return res
        end
        function ($f){S, T}(a::DataMatrix{S}, b::DataMatrix{T})
            res = DataArray(Array(Bool, size(a)), BitArray(size(a)))
            for i in 1:length(a)
                if isna(a[i]) || isna(b[i])
                    res[i] = NA
                else
                    res[i] = ($scalarf)(a[i], b[i])
                end
            end
            return res
        end
        function ($f){S, T}(a::DataMatrix{S}, b::Matrix{T})
            res = DataArray(Array(Bool, size(a)), BitArray(size(a)))
            for i in 1:length(a)
                if isna(a[i])
                    res[i] = NA
                else
                    res[i] = ($scalarf)(a[i], b[i])
                end
            end
            return res
        end
        function ($f){S, T}(a::Matrix{S}, b::DataMatrix{T})
            res = DataArray(Array(Bool, size(a)), BitArray(size(a)))
            for i in 1:length(a)
                if isna(b[i])
                    res[i] = NA
                else
                    res[i] = ($scalarf)(a[i], b[i])
                end
            end
            return res
        end
        function ($f)(a::DataFrame, b::DataFrame)
            n, p = nrow(a), ncol(a)
            if n != nrow(b) || p != ncol(b)
                error("DataFrames must have matching sizes for comparisons")
            end
            # Tries to preserve types from a
            results = DataFrame(Array(Bool, size(a)))
            # TODO: Test that types match across a and b
            for j in 1:p
                for i in 1:n
                    results[i, j] = ($scalarf)(a[i, j], b[i, j])
                end
            end
            return results
        end
    end
end

for f in binary_operators
    @eval begin
        function ($f)(d::NAtype, e::NAtype)
            return NA
        end
        function ($f){T <: Union(Number, String)}(d::NAtype, x::T)
            return NA
        end
        function ($f){T <: Union(Number, String)}(x::T, d::NAtype)
            return NA
        end
        function ($f){T}(A::Number, B::DataVector{T})
            res = DataArray(Array(promote_type(typeof(A),T), length(B)), B.na)
            for i in 1:length(B)
                res.data[i] = ($f)(A, B.data[i])
            end
            return res
        end
        function ($f){T}(A::DataVector{T}, B::Number)
            res = DataArray(Array(promote_type(typeof(B),T), length(A)), A.na)
            for i in 1:length(A)
                res.data[i] = ($f)(A.data[i], B)
            end
            return res
        end
        function ($f){T}(A::NAtype, B::DataVector{T})
            res = DataArray(Array(typeof(B).parameters[1], length(B)), B.na)
            for i in 1:length(B)
                res.data[i] = B.data[i]
                res.na[i] = true
            end
            return res
        end
        function ($f){T}(A::DataVector{T}, B::NAtype)
            res = DataArray(Array(typeof(A).parameters[1], length(A)), A.na)
            for i in 1:length(A)
                res.data[i] = A.data[i]
                res.na[i] = true
            end
            return res
        end
        function ($f)(df::DataFrame, x::Union(Number, NAtype))
            n, p = nrow(df), ncol(df)
            # Tries to preserve types
            results = deepcopy(df)
            for j in 1:p
                if typeof(df[j]).parameters[1] <: Number
                    for i in 1:n
                        results[i, j] = ($f)(df[i, j], x)
                    end
                else
                    for i in 1:n
                        results[i, j] = NA
                    end
                end
            end
            return results
        end
        function ($f)(x::Union(Number, NAtype), df::DataFrame)
            n, p = nrow(df), ncol(df)
            # Tries to preserve types
            results = deepcopy(df)
            for j in 1:p
                if typeof(df[j]).parameters[1] <: Number
                    for i in 1:n
                        results[i, j] = ($f)(x, df[i, j])
                    end
                else
                    for i in 1:n
                        results[i, j] = NA
                    end
                end
            end
            return results
        end
        function ($f){T <: Union(String, Number)}(d::NAtype, x::T)
            return NA
        end
        function ($f){T <: Union(String, Number)}(x::T, d::NAtype)
            return NA
        end
        function ($f){T,N}(a::AbstractArray{T,N}, d::NAtype)
            return NA
        end
        function ($f){T,N}(d::NAtype, a::AbstractArray{T,N})
            return NA
        end
    end
end

for f in induced_binary_operators
    @eval begin
        function ($f)(d::NAtype, e::NAtype)
            return NA
        end
        function ($f)(d::Union(String, Number), e::NAtype)
            return NA
        end
        function ($f)(d::NAtype, e::FloatingPoint)
            return NA
        end
    end
end

# for arithmetic, NAs propagate
for f in array_arithmetic_operators
    @eval begin
        function ($f){S, T}(A::DataVector{S}, B::Vector{T})
            n_A, n_B = length(A), length(B)
            if n_A != n_B
                error("DataVector and Vector lengths must match")
            end
            res = DataArray(Array(promote_type(S, T), n_A), BitArray(n_A))
            for i in 1:n_A
                res.na[i] = A.na[i]
                res.data[i] = ($f)(A.data[i], B[i])
            end
            return res
        end
        function ($f){S, T}(A::Vector{S}, B::DataVector{T})
            n_A, n_B = length(A), length(B)
            if n_A != n_B
                error("Vector and DataVector lengths must match")
            end
            res = DataArray(Array(promote_type(S, T), n_A), BitArray(n_A))
            for i in 1:n_A
                res.na[i] = B.na[i]
                res.data[i] = ($f)(A[i], B.data[i])
            end
            return res
        end
        function ($f){S, T}(A::DataVector{S}, B::DataVector{T})
            if (length(A) != length(B))
                error("DataVector lengths must match")
            end
            res = DataArray(Array(promote_type(S, T),
                                  length(A)),
                            BitArray(length(A)))
            for i in 1:length(A)
                res.na[i] = (A.na[i] || B.na[i])
                res.data[i] = ($f)(A.data[i], B.data[i])
            end
            return res
        end
        function ($f){S, T}(A::DataMatrix{S}, B::DataMatrix{T})
            if size(A) != size(B)
                error("DataMatrix sizes must match")
            end
            res = DataArray(Array(promote_type(S, T), size(A)),
                            BitArray(size(A)))
            for i in 1:length(A)
                res.na[i] = (A.na[i] || B.na[i])
                res.data[i] = ($f)(A.data[i], B.data[i])
            end
            return res
        end
        function ($f)(a::DataFrame, b::DataFrame)
            n, p = nrow(a), ncol(a)
            if n != nrow(b) || p != ncol(b)
                error("DataFrames must have matching sizes for arithmetic")
            end
            # Tries to preserve types from a
            results = deepcopy(a)
            for j in 1:p
                if typeof(a[j]).parameters[1] <: Number
                    for i in 1:n
                        results[i, j] = ($f)(a[i, j], b[i, j])
                    end
                else
                    for i in 1:n
                        results[i, j] = NA
                    end
                end
            end
            return results
        end
    end
end

for f in biscalar_operators
    @eval begin
        function ($f)(d::NAtype, e::NAtype)
            return NA
        end
        function ($f){T <: Number}(d::NAtype, x::T)
            return NA
        end
        function ($f){T <: Number}(x::T, d::NAtype)
            return NA
        end
    end
end

for f in pairwise_vector_operators
    @eval begin
        function ($f)(dv::DataVector)
            n = length(dv)
            new_data = ($f)(dv.data)
            new_na = falses(n - 1)
            for i = 2:(n - 1)
                if isna(dv[i])
                    new_na[i - 1] = true
                    new_na[i] = true
                end
            end
            if isna(dv[n])
                new_na[n - 1] = true
            end
            return DataArray(new_data, new_na)
        end
    end
end

for f in cumulative_vector_operators
    @eval begin
        function ($f)(dv::DataVector)
            new_data = ($f)(dv.data)
            new_na = falses(length(dv))
            hitna = false
            for i = 1:length(dv)
                if isna(dv[i])
                    hitna = true
                end
                if hitna
                    new_na[i] = true
                end
            end
            return DataArray(new_data, new_na)
        end
    end
end

for f in unary_vector_operators
    @eval begin
        function ($f)(dv::DataVector)
            if any(isna(dv))
                return NA
            else
                return ($f)(dv.data)
            end
        end
    end
end

for f in ffts
    @eval begin
        function ($f)(dv::DataVector)
            if any(isna(dv))
                return NA
            else
                return ($f)(dv.data)
            end
        end
    end
end

for f in binary_vector_operators
    @eval begin
        function ($f)(dv1::DataVector, dv2::DataVector)
            if any(isna(dv1)) || any(isna(dv2))
                return NA
            else
                return ($f)(dv1.data, dv2.data)
            end
        end
    end
end

for (f, colf) in ((:min, :colmins),
                  (:max, :colmaxs),
                  (:prod, :colprods),
                  (:sum, :colsums),
                  (:mean, :colmeans),
                  (:median, :colmedians),
                  (:std, :colstds),
                  (:var, :colvars),
                  (:fft, :colffts), # TODO: Remove and/or fix
                  (:norm, :colnorms))
    @eval begin
        function ($colf)(df::DataFrame)
            p = ncol(df)
            res = DataFrame()
            for j in 1:p
                res[j] = DataArray(($f)(df[j]))
            end
            colnames!(res, colnames(df))
            return res
        end
        function ($colf)(dm::DataMatrix)
            n, p = nrow(dm), ncol(dm)
            res = datazeros(p)
            for j in 1:p
                res[j] = ($f)(dm[:, j])
            end
            return res
        end
    end
end

for (f, rowf) in ((:min, :rowmins),
                  (:max, :rowmaxs),
                  (:prod, :rowprods),
                  (:sum, :rowsums),
                  (:mean, :rowmeans),
                  (:median, :rowmedians),
                  (:std, :rowstds),
                  (:var, :rowvars),
                  (:fft, :rowffts), # TODO: Remove and/or fix
                  (:norm, :rownorms))
    @eval begin
        function ($rowf)(dm::DataMatrix)
            n, p = nrow(dm), ncol(dm)
            res = datazeros(n)
            for i in 1:n
                res[i] = ($f)(DataArray(reshape(dm.data[i, :], p), reshape(dm.na[i, :], p)))
            end
            return res
        end
    end
end

#
# Boolean operators
#

function all{T}(dv::DataArray{T})
    for i in 1:length(dv)
        if isna(dv[i])
            return NA
        end
        if !dv[i]
            return false
        end
    end
    return true
end

function any{T}(dv::DataArray{T})
    has_na = false
    for i in 1:length(dv)
        if !isna(dv[i])
            if dv[i]
                return true
            end
        else
            has_na = true
        end
    end
    if has_na
        return NA
    else
        return false
    end
end

function all(df::DataFrame)
    for i in 1:nrow(df)
        for j in 1:ncol(df)
            if isna(df[i, j])
                return NA
            end
            if !df[i, j]
                return false
            end
        end
    end
    return true
end

function any(df::DataFrame)
    has_na = false
    for i in 1:nrow(df)
        for j in 1:ncol(df)
            if !isna(df[i, j])
                if df[i, j]
                    return true
                end
            else
                has_na = true
            end
        end
    end
    if has_na
        return NA
    else
        return false
    end
end

# isequal() should for Data*
# * If missingness differs, underlying values are irrelevant
# * If both entries are NA, underlying values are irrelevant
function isequal{T}(a::DataVector{T}, b::DataVector{T})
    if length(a) != length(b)
        return false
    else
        for i = 1:length(a)
            if a.na[i] != b.na[i]
                return false
            elseif !a.na[i] && !b.na[i] && (a.data[i] != b.data[i])
                return false
            end
        end
    end
    return true
end
function isequal{T}(a::PooledDataVector{T}, b::PooledDataVector{T})
    if length(a) != length(b)
        return false
    else
        for i = 1:length(a)
            # Will we speed this up by looking under hood?
            if isna(a[i])
                if !isna(b[i])
                    return false
                end
            else
                if isna(b[i])
                    return false
                end
                if a[i] != b[i]
                    return false
                end
            end
        end
    end
    return true
end
function isequal{T}(a::AbstractDataVector{T}, b::AbstractDataVector{T})
    if length(a) != length(b)
        return false
    else
        for i = 1:length(a)
            if isna(a[i])
                if !isna(b[i])
                    return false
                end
            else
                if isna(b[i])
                    return false
                end
                if a[i] != b[i]
                    return false
                end
            end
        end
    end
    return true
end
function isequal(df1::AbstractDataFrame, df2::AbstractDataFrame)
    if ncol(df1) != ncol(df2)
        return false
    end
    for idx in 1:ncol(df1)
        if !isequal(df1[idx], df2[idx])
            return false
        end
    end
    return true
end

function range{T}(dv::DataVector{T})
    return DataArray([min(dv), max(dv)], falses(2))
end

function rle{T}(v::DataVector{T})
    n = length(v)
    current_value = v[1]
    current_length = 1
    values = DataArray(T, n)
    total_values = 1
    lengths = Array(Int16, n)
    total_lengths = 1
    for i in 2:n
        if isna(v[i]) || isna(current_value)
            if isna(v[i]) && isna(current_value)
                current_length += 1
            else
                values[total_values] = current_value
                total_values += 1
                lengths[total_lengths] = current_length
                total_lengths += 1
                current_value = v[i]
                current_length = 1
            end
        else
            if v[i] == current_value
                current_length += 1
            else
                values[total_values] = current_value
                total_values += 1
                lengths[total_lengths] = current_length
                total_lengths += 1
                current_value = v[i]
                current_length = 1
            end
        end
    end
    values[total_values] = current_value
    lengths[total_lengths] = current_length
    return (values[1:total_values], lengths[1:total_lengths])
end

## inverse run-length encoding
function inverse_rle{T}(values::DataVector{T}, lengths::Vector{Int16})
    total_n = sum(lengths)
    pos = 0
    res = DataArray(T, total_n)
    n = length(values)
    for i in 1:n
        v = values[i]
        l = lengths[i]
        for j in 1:l
            pos += 1
            res[pos] = v
        end
    end
    return res
end


text:
struct EnumerableGroupBySimple{T,TKey,TS,SO,ES<:Function} <: Enumerable
    source::SO
    elementSelector::ES
end

struct Grouping{TKey,T} <: AbstractArray{T,1}
    _key::TKey
    elements::Array{T,1}
end

key(g::Grouping) = getfield(g, :_key)

struct GroupColumnArrayView{T,G,INDEX} <: AbstractVector{T}
    grouping::G
end

function Base.size(gcav::GroupColumnArrayView)
    return size(gcav.grouping)
end

function Base.getindex(gcav::GroupColumnArrayView{T,G,INDEX}, i::Int) where {T,G,INDEX}
    return getproperty(gcav.grouping[i],INDEX)
end

Base.IndexStyle(::Type{GroupColumnArrayView}) = IndexLinear()

function Base.getproperty(g::Grouping{TKey,T}, name::Symbol) where {TKey,T}

    return GroupColumnArrayView{fieldtype(T,name),Grouping{TKey,T},name}(g)
end

Base.size(A::Grouping{TKey,T}) where {TKey,T} = size(getfield(A, :elements))
Base.getindex(A::Grouping{TKey,T},i) where {TKey,T} = getfield(A, :elements)[i]
Base.length(A::Grouping{TKey,T}) where {TKey,T} = length(getfield(A, :elements))

Base.eltype(::Type{EnumerableGroupBySimple{T,TKey,TS,SO,ES}}) where {T,TKey,TS,SO,ES} = T

function groupby(source::Enumerable, f_elementSelector::Function, elementSelector::Expr)
    TS = eltype(source)
    TKey = Base._return_type(f_elementSelector, Tuple{TS,})

    SO = typeof(source)

    T = Grouping{TKey,TS}

    ES = typeof(f_elementSelector)

    return EnumerableGroupBySimple{T,TKey,TS,SO,ES}(source,f_elementSelector)
end

function Base.iterate(iter::EnumerableGroupBySimple{T,TKey,TS,SO,ES}) where {T,TKey,TS,SO,ES}
    result = OrderedDict{TKey,T}()
    for i in iter.source
        key = iter.elementSelector(i)
        if !haskey(result, key)
            result[key] = T(key,Array{TS}(undef, 0))
        end
        push!(getfield(result[key], :elements),i)
    end

    groups = collect(values(result))
    if length(groups)==0
        return nothing
    else
        return groups[1], (groups, 2)
    end
end

function Base.iterate(iter::EnumerableGroupBySimple{T,TKey,TS,SO,ES}, state) where {T,TKey,TS,SO,ES}
    if state[2]>length(state[1])
        return nothing
    else
        return state[1][state[2]], (state[1], state[2]+1)
    end
end

struct EnumerableGroupBy{T,TKey,TR,SO,ES<:Function,RS<:Function} <: Enumerable
    source::SO
    elementSelector::ES
    resultSelector::RS
end

Base.eltype(::Type{EnumerableGroupBy{T,TKey,TR,SO,ES, RS}}) where {T,TKey,TR,SO,ES,RS} = T

function groupby(source::Enumerable, f_elementSelector::Function, elementSelector::Expr, f_resultSelector::Function, resultSelector::Expr)
    TS = eltype(source)
    TKey = Base._return_type(f_elementSelector, Tuple{TS,})

    SO = typeof(source)

    TR = Base._return_type(f_resultSelector, Tuple{TS,})

    T = Grouping{TKey,TR}

    ES = typeof(f_elementSelector)
    RS = typeof(f_resultSelector)

    return EnumerableGroupBy{T,TKey,TR,SO,ES,RS}(source,f_elementSelector,f_resultSelector)
end

function Base.iterate(iter::EnumerableGroupBy{T,TKey,TR,SO,ES}) where {T,TKey,TR,SO,ES}
    result = OrderedDict{TKey,T}()
    for i in iter.source
        key = iter.elementSelector(i)
        if !haskey(result, key)
            result[key] = T(key,Array{TR}(undef,0))
        end
        push!(getfield(result[key], :elements),iter.resultSelector(i))
    end

    groups = collect(values(result))
    if length(groups)==0
        return nothing
    else
        return groups[1], (groups, 2)
    end
end

function Base.iterate(iter::EnumerableGroupBy{T,TKey,TR,SO,ES}, state) where {T,TKey,TR,SO,ES}
    if state[2]>length(state[1])
        return nothing
    else
        return state[1][state[2]], (state[1], state[2]+1)
    end
end


text:
"""
    dataloader(name; header=true, dir="HorseMLdatasets")
Load a data for Machine Learning. `name` is either the name of the datasets or the full path of the data file to be loaded. The following three can be specified as the name of the dataset in `name`.
- `MNIST` : The MNIST Datasets
- `iris` : The iris Datasets
- `BostonHousing` : Boston Housing DataSets
And these datasets are downloaded and saved by creating a `dir` folder under the home directly(i.e. it is saved in the `/home_directory/HorseMLdatasets` by default).
When importing a data file, you can specify whether to read the header with `header`.

# Example
```@repl
julia> dataloader("MNIST");

julia> dataloader("/home/ubuntu/data/data.csv", header = false)
```
"""
function dataloader(name; header = true, dir = "HorseMLdatasets")
    hd = homedir()
    dir = joinpath(hd, dir)
    mkpath(dir)
    current = pwd()
    cd(dir)
    if name == "MNIST"
        if !(isfile("mnist_train.csv"))
            yn = Base.prompt("Can I download the mnist_train.csv(110MB)? (y/n)")
            yn == "y" && Downloads.download("https://pjreddie.com/media/files/mnist_train.csv", "mnist_train.csv")
        end
        df_r = CSV.read("mnist_train.csv", header = false, DataFrame)
        if !(isfile("mnist_test.csv"))
            yn = Base.prompt("Can I download the mnist_test.csv(18MB)? (y/n)")
            yn == "y" && Downloads.download("https://pjreddie.com/media/files/mnist_test.csv", "mnist_test.csv")
        end
        df_e = CSV.read("mnist_test.csv",header = false , DataFrame)
        cd(current)
        return df_r, df_e
    elseif name == "iris"
        url = "https://gist.githubusercontent.com/netj/8836201/raw/6f9306ad21398ea43cba4f7d537619d0e07d5ae3/iris.csv"
        if !(isfile("iris.csv"))
            yn = Base.prompt("Can I download the iris.csv(4KB)? (y/n)")
            yn == "y" && Downloads.download(url, "iris.csv")
        end
        df = CSV.read("iris.csv", DataFrame)
        cd(current)
        return df
    elseif name == "BostonHousing"
        url = "https://raw.githubusercontent.com/selva86/datasets/master/BostonHousing.csv"
        if !(isfile("BostonHousing.csv"))
            yn = Base.prompt("Can I download the BostonHousing.csv(36KB)? (y/n)")
            yn == "y" && Downloads.download(url, "BostonHousing.csv")
        end
        df = CSV.read("BostonHousing.csv", DataFrame)
        cd(current)
        return df
    end
    try
        df = CSV.read(name, header=header, DataFrame)
        return df
    finally
        cd(current)
    end
end

"""
    databuilder(x, y; batches=1)
`x` is the feature value, `y` is the teacher data, and `bacthes` is the batch size.
This function formats and returns the data used for the neural network(however, `x` and` y` must be Arrays or DataFrames, so please use this function  after encoding, normalization, etc.).

# Example
```jldoctest
julia> data = dataloader("iris");

julia> LE, OHE = LabelEncoder(), OneHotEncoder()
(LabelEncoder(Dict{Any, Any}()), OneHotEncoder())

julia> x, y = data[:, Not(:variety)], OHE(LE(data[:, :variety]));

julia> databuilder(x, y)
150-element Vector{Tuple{Matrix{Float64}, Matrix{Float64}}}:
 ([5.1; 3.5; 1.4; 0.2;;], [1.0; 0.0; 0.0;;])
 ([4.9; 3.0; 1.4; 0.2;;], [1.0; 0.0; 0.0;;])
 ([4.7; 3.2; 1.3; 0.2;;], [1.0; 0.0; 0.0;;])
 ([4.6; 3.1; 1.5; 0.2;;], [1.0; 0.0; 0.0;;])
 ([5.0; 3.6; 1.4; 0.2;;], [1.0; 0.0; 0.0;;])
 ([5.4; 3.9; 1.7; 0.4;;], [1.0; 0.0; 0.0;;])
 ([4.6; 3.4; 1.4; 0.3;;], [1.0; 0.0; 0.0;;])
 ([5.0; 3.4; 1.5; 0.2;;], [1.0; 0.0; 0.0;;])
 ([4.4; 2.9; 1.4; 0.2;;], [1.0; 0.0; 0.0;;])
 ([4.9; 3.1; 1.5; 0.1;;], [1.0; 0.0; 0.0;;])
 ([5.4; 3.7; 1.5; 0.2;;], [1.0; 0.0; 0.0;;])
 ([4.8; 3.4; 1.6; 0.2;;], [1.0; 0.0; 0.0;;])
 ([4.8; 3.0; 1.4; 0.1;;], [1.0; 0.0; 0.0;;])
 ⋮
 ([6.0; 3.0; 4.8; 1.8;;], [0.0; 0.0; 1.0;;])
 ([6.9; 3.1; 5.4; 2.1;;], [0.0; 0.0; 1.0;;])
 ([6.7; 3.1; 5.6; 2.4;;], [0.0; 0.0; 1.0;;])
 ([6.9; 3.1; 5.1; 2.3;;], [0.0; 0.0; 1.0;;])
 ([5.8; 2.7; 5.1; 1.9;;], [0.0; 0.0; 1.0;;])
 ([6.8; 3.2; 5.9; 2.3;;], [0.0; 0.0; 1.0;;])
 ([6.7; 3.3; 5.7; 2.5;;], [0.0; 0.0; 1.0;;])
 ([6.7; 3.0; 5.2; 2.3;;], [0.0; 0.0; 1.0;;])
 ([6.3; 2.5; 5.0; 1.9;;], [0.0; 0.0; 1.0;;])
 ([6.5; 3.0; 5.2; 2.0;;], [0.0; 0.0; 1.0;;])
 ([6.2; 3.4; 5.4; 2.3;;], [0.0; 0.0; 1.0;;])
 ([5.9; 3.0; 5.1; 1.8;;], [0.0; 0.0; 1.0;;])
```
"""
function databuilder(x::Matrix{TX}, t::Matrix{TY}; batches = 1) where {TX, TY}
    PT = promote_type(TX, TY)
    x, t = Matrix{PT}(x), Matrix{PT}(t)
    size(x, 1) % batches != 0 && @warn "the data size is not divisible by the batch size! some data will be omitted"
    N = round(Int64, size(x, 1) / batches)
    x, t = transpose(x), transpose(t)
    data = Vector{Tuple{Matrix{TX}, Matrix{TY}}}(undef, N)
    c = 1
    for n in 1 : N
        data[n] = (x[:, c:n*batches], t[:, c:n*batches])
        c += batches
    end
    return data
end

function databuilder(x::Array{TX, 4}, t::Matrix{TY}; batches = 1) where {TX, TY}
    PT = promote_type(TX, TY)
    x, t = Array{PT, 4}(x), Matrix{PT}(t)
    size(x, 1) % batches != 0 && @warn "the data size is not divisible by the batch size! some data will be omitted"
    N = round(Int64, size(x, 1) / batches)
    x, t = permutedims(x, (2, 3, 4, 1)), transpose(t)
    data = Vector{Tuple{Array{PT, 4}, Matrix{PT}}}(undef, N)
    c = 1
    for n in 1 : N
        data[n] = (x[:, :, :, c:n*batches], t[:, c:n*batches])
        c += batches
    end
    return data
end

databuilder(x::Array, t::Vector; batches = 1) = databuilder(x, t[:, :], batches = batches)
databuilder(x::DataFrame, t; batches=1) = databuilder(Matrix(x), t, batches = batches)
databuilder(x::DataFrame, t::DataFrame; batches=1) = databuilder(Matrix(x), Matrix(t), batches = batches)

"""
    DataSplitter(ndata; test_size=nothing, train_size=nothing)
Split the data into test data and training data. `ndata` is the number of the data, and you must specify either `test_size` or `train_size`. thease parameter can be proprtional or number of data.

!!! note
    If both `test_size` and `train_size` are specified, `test_size` takes precedence.

#Example
```jldoctest
julia> x = rand(20, 2);

julia> DS = DataSplitter(50, train_size = 0.3);

julia> train, test = DS(x, dims = 2);

julia> train |> size
(6, 2)

julia> test |> size
(14, 2)
```
"""
struct DataSplitter
    #TODO: make it possible to DataFrame split
    test_indices::Array{Int64, 1}
    train_indices::Array{Int64, 1}
    function DataSplitter(ndata; test_size=nothing, train_size=nothing)
        test_size == train_size == nothing && throw(ArgumentError("test_size and train_size wasn't specified. you must specify either one."))
        if test_size == nothing
            #number of train data is ceiled
            test_size = (0<=train_size<=1) ? ndata-ceil(train_size*ndata) : ndata-train_size
        else
            test_size = (0<=test_size<=1) ? ndata-ceil((1-test_size)*ndata) : test_size
        end
        new(sample(1:ndata, Int(test_size), nc=true)...)
    end
end

function (DS::DataSplitter)(xs; dims=1)
    test_index = fill!(Array{Union{Colon, Array{Int64, 1}}}(undef, ndims(xs)), :)
    train_index = fill!(Array{Union{Colon, Array{Int64, 1}}}(undef, ndims(xs)), :)
    train_index[dims] = DS.train_indices
    test_index[dims] = DS.test_indices
    return xs[train_index...], xs[test_index...]
end


text:
include("load_cheaplogic.jl")

magia=readclausefromfile("data/magia.cnf")
cmagi=createcore(magia)

data001=readclausefromfile("data/data001.cnf")
cd001=createcore(data001)

ctime=readcore("data/time.cnf")
printcore(ctime)

max_stars_repo_name: sazare/cheapViews

max_issues_repo_name: sazare/cheapViews

max_forks_repo_name: sazare/cheapViews


text:
using SEC, Test, GLMakie, LinearAlgebra
import SEC.Emulate: mean, cov, GaussianProcess, condition

μ(x) = 0.0
δ(x,y) = x==y ? 1 : 0
k(x,y) = 10exp(-(x-y)^2/2) + 0.1 * δ(x,y) # noise here
gp = GaussianProcess(mean = μ, covariance = k)
N = 10
X = [i for i in 1:N]
Y = [2*i for i in 1:N]
data = (; X, Y)
cgp = condition(gp, data, noise = 0.0) # noise here
# philosophical difference only
# just associates more noise at the given prediction point
##
fig = Figure(resolution = (700, 450))
ax = Axis(fig, xlabel = "x", ylabel = "y")
Xn = collect(range(0, N+1, 100))
Xn = setdiff(Xn, X)
μX = mean(cgp, Xn)
σXX =  [cgp.covariance(x,x) for x in Xn]
scatter!(ax, X, Y, color = :red)
lines!(ax, Xn, μX, color = :blue, linewidth =3)
# filled curve 1
band!(ax, Xn, μX - 3*σXX, μX + 3*σXX, color = (:blue, 0.2))
xlims!(ax, (0,11))
ylims!(ax, (-1,21))
ax.xlabelsize = 25
ax.xticklabelsize = 25
ax.ylabelsize = 25
ax.yticklabelsize = 25
update!(fig.scene)
fig[1,1] = ax
fig

max_stars_repo_name: bischtob/SEC

max_issues_repo_name: bischtob/SEC

max_forks_repo_name: bischtob/SEC


text:
using StructuralCausalModels
#using StatsPlots
#gr(size=(700,800))

ProjDir = @__DIR__
cd(ProjDir) #do

N = 100
b_AU = 0.5
b_AC = 3
b_UX = -1
b_UB = 2
b_CB = -1.5
b_CY = 1
b_XY = 5

df = DataFrame(
  :a => rand(Normal(), N)
);
df[!, :u] = rand(Normal(1, 2), N) + b_AU * df[:, :a]
df[!, :c] = rand(Normal(0, 1), N) + b_AC * df[:, :a]
df[!, :b] = rand(Normal(1, 1), N) + b_UB * df[:, :u]
df[!, :x] = rand(Normal(-2, 1), N) + b_UX * df[:, :u]
df[!, :y] = rand(Normal(1, 2), N) + b_XY * df[:, :x] + b_CY * df[:, :c]

#StatsPlots.cornerplot(Array(df), label=names(df))
#savefig("$(ProjDir)/SR6.4.2.png")

d = OrderedDict(
  :u => :a,
  :c => :a,
  :b => [:u, :c],
  :y => :c,
  :x => :u
);
u = [:u]

dag = DAG("sr6_4_2", d, df=df);
show(dag)

fname = joinpath(mktempdir(), "sr6.4.2.dot")
to_graphviz(dag, fname)
Sys.isapple() && run(`open -a GraphViz.app $(fname)`)
to_dagitty(dag) |> display
to_ggm(dag) |> display

display(dag.s); println()

bs = basis_set(dag)
display(bs); println()

t = shipley_test(dag)
display(t); println()

f = [:a]; s = [:b]; conditioning_set = [:u, :c]

e = d_separation(dag, f, s, c=conditioning_set)
println("d_separation($(dag.name), $f, $s, c=$conditioning_set) = $e")
println()

adjustmentsets = adjustment_sets(dag, :x, :y)
println("\nAdjustment sets:")
adjustmentsets |> display

#=
adjustmentsets = adjustment_sets(dag, :x, :y, u)
println("\nAdjustment sets:")
adjustmentsets |> display
=#

#end


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

"""
    georef(table, domain)

Georeference `table` on spatial `domain`.
"""
georef(table, domain) = SpatialData(domain, table)

"""
    georef(table, coords)

Georeference `table` on a `PointSet(coords)`.
"""
georef(table, coords::AbstractVecOrMat) = georef(table, PointSet(coords))

"""
    georef(table, coordnames)

Georeference `table` using columns `coordnames`.
"""
function georef(table, coordnames::NTuple)
  ctor = Tables.materializer(table)
  colnames = Tables.columnnames(table)
  @assert coordnames ⊆ colnames "invalid coordinates for table"
  @assert !(colnames ⊆ coordnames) "table must have at least one variable"
  varnames = setdiff(colnames, coordnames)
  vars = [v => Tables.getcolumn(table, v) for v in varnames]
  coords = reduce(hcat, [Tables.getcolumn(table, c) for c in coordnames])
  georef(ctor(vars), coords')
end

"""
    georef(tuple, domain)

Georeference named `tuple` on spatial `domain`.
"""
georef(tuple::NamedTuple, domain) =
  georef(DataFrame([var=>vec(val) for (var,val) in pairs(tuple)]), domain)

"""
    georef(tuple, coords)

Georefrence named `tuple` on `PointSet(coords)`.
"""
georef(tuple::NamedTuple, coords::AbstractVecOrMat) = georef(tuple, PointSet(coords))

"""
    georef(tuple; origin=(0.,0.,...), spacing=(1.,1.,...))

Georeference named `tuple` on `RegularGrid(size(tuple[1]), origin, spacing)`.
"""
georef(tuple;
       origin=ntuple(i->0., ndims(tuple[1])),
       spacing=ntuple(i->1., ndims(tuple[1]))) = georef(tuple, origin, spacing)

georef(tuple, origin, spacing) =
  georef(tuple, RegularGrid(size(tuple[1]), origin, spacing))


text:
"""
    pairclustering(𝚽::Matrix{Float64}, GP_star::GraphPart)

construct the GraphPart object of the primal graph via pair-clustering algorithm.

# Input Arguments
- `𝚽::Matrix{Float64}`: graph Laplacian eigenvectors 𝚽
- `GP_star::GraphPart`: GraphPart object of the dual graph

# Output Argument
- `GP::GraphPart`: GraphPart object of the primal graph via pair-clustering

"""
function pairclustering(𝚽::Matrix{Float64}, GP_star::GraphPart)
    rs_dual = GP_star.rs
    inds_dual = GP_star.inds
    (N, jmax) = Base.size(inds_dual)

    # TRACKING VARIABLES

    # `ind` records the way in which the nodes are indexed on each level
    ind = Vector{Int}(1:N)

    # `inds` records the way in which the nodes are indexed on all levels
    inds = zeros(Int, N, jmax)
    inds[:, 1] = ind

    # `rs` stands for regionstart, meaning that the index in `ind` of the first
    # point in region number `i` is `rs[i]`
    rs = deepcopy(rs_dual)

    j = 1
    regioncount = 0
    while regioncount < N
        regioncount = count(!iszero, rs[:, j]) - 1
        if j == jmax  # add a column to rs for level j+1, if necessary
            rs = hcat(rs, vcat(Int(1), zeros(Int, N)))
            inds = hcat(inds, zeros(Int, N))
            jmax = jmax + 1
        end
        # for tracking the child regions
        rr = 1
        for r = 1:regioncount
            rs1 = rs[r, j]
            rs2 = rs[r + 1, j]
            n = rs2 - rs1
            if n > 1            # regions with 2 or more nodes
                indrs = ind[rs1:(rs2 - 1)]
                n1 = rs[rr + 1, j + 1] - rs1
                indrs_dual1 = inds_dual[rs1:(rs1 + n1 - 1), j + 1]
                indrs_dual2 = inds_dual[(rs1 + n1):(rs2 - 1), j + 1]
                # partition the current region
                (indp1, indp2) = partition_pc(𝚽, indrs, indrs_dual1, indrs_dual2)
                # update the indexing
                ind[rs1:(rs1 + n1 - 1)] = indrs[indp1]
                ind[(rs1 + n1):(rs2 - 1)] = indrs[indp2]
                # update the child region tracking
                rr += 2
            elseif n == 1
                rr += 1
            end
        end
        j = j + 1
        inds[:, j] = ind
    end

    # get rid of excess columns in rs
    rs = rs[:, 1:(j - 1)]
    # get rid of excess columns in inds
    inds = inds[:, 1:(j - 1)]

    return GraphPart(ind, rs; inds = inds)
end

"""
    partition_pc(𝚽::Matrix{Float64}, indrs::Vector{Int}, indrs_dual1::Vector{Int}, indrs_dual2::Vector{Int})

based on the partition on dual graph, partition the primal graph region
(indexed by `indrs`) into two pieces via the pair-clustering algorithm.

# Input Arguments
- `𝚽::Matrix{Float64}`: graph Laplacian eigenvectors 𝚽
- `indrs::Vector{Int}`: indices of the primal graph region to be partitioned
- `indrs_dual1::Vector{Int}`: index1 of the dual graph region's partition result
- `indrs_dual2::Vector{Int}`: index2 of the dual graph region's partition result

# Output Argument
- `indp1::Vector{Int}`: index1 of `indrs` for the partition result
- `indp2::Vector{Int}`: index2 of `indrs` for the partition result

"""
function partition_pc(𝚽::Matrix{Float64}, indrs::Vector{Int}, indrs_dual1::Vector{Int}, indrs_dual2::Vector{Int})
    n = length(indrs)
    n1 = length(indrs_dual1)
    E = 𝚽[indrs, :].^2
    score = sum(E[:, indrs_dual1], dims = 2) - sum(E[:, indrs_dual2], dims = 2)
    indp1 = sortperm(score[:]; rev = true)[1:n1]
    indp2 = setdiff(1:n, indp1)
    return indp1, indp2
end

"""
    pc_ngwp(𝚽::Matrix{Float64}, GP_star::GraphPart, GP::GraphPart)

construct pair-clustering NGWP and GP.tag in place.

# Input Arguments
- `𝚽::Matrix{Float64}`: graph Laplacian eigenvectors 𝚽
- `GP_star::GraphPart`: GraphPart object of the dual graph
- `GP::GraphPart`: GraphPart object of the primal graph

# Output Argument
- `wavelet_packet::Array{Float64,3}`: the pair-clustering NGWP. The first
    index is for selecting wavelets at a fixed level; the second index is for
    selecting the level `j`; the third index is for selecting elements in the
    wavelet vector.

"""
function pc_ngwp(𝚽::Matrix{Float64}, GP_star::GraphPart, GP::GraphPart)
    rs = GP_star.rs
    inds_dual = GP_star.inds
    inds_primal = GP.inds
    (N, jmax) = Base.size(inds_dual)

    GP_star.tag = zeros(Int, N, jmax)
    GP_star.tag[:, 1] = Vector{Int}(0:(N - 1))

    wavelet_packet = zeros(N, jmax, N)
    wavelet_packet[:, 1, :] = Matrix{Float64}(I, N, N)

    for j = 2:jmax
        regioncount = count(!iszero, rs[:, j]) - 1
        for r = 1:regioncount
            indr = rs[r, j]:(rs[r + 1, j] - 1)
            GP_star.tag[indr, j] = Vector{Int}(0:(length(indr) - 1))
            wavelet_packet[indr, j, :] = const_proj_wavelets(𝚽, inds_primal[indr, j], inds_dual[indr, j])'
        end
    end

    return wavelet_packet
end


text:
#=
    MadNLP utils
=#

@kernel function _copy_diag!(dest, src)
    i = @index(Global)
    dest[i] = src[i, i]
end

function MadNLP.diag!(dest::CuVector{T}, src::CuMatrix{T}) where T
    @assert length(dest) == size(src, 1)
    ev = _copy_diag!(CUDADevice())(dest, src, ndrange=length(dest))
    wait(ev)
end

@kernel function _add_diagonal!(dest, src1, src2)
    i = @index(Global)
    dest[i, i] = src1[i] + src2[i]
end

function MadNLP.diag_add!(dest::CuMatrix, src1::CuVector, src2::CuVector)
    ev = _add_diagonal!(CUDADevice())(dest, src1, src2, ndrange=size(dest, 1))
    wait(ev)
end

#=
    MadNLP kernels
=#

# Overload is_valid to avoid fallback to default is_valid, slow on GPU
MadNLP.is_valid(src::CuArray) = true

# Constraint scaling
function MadNLP.set_con_scale!(con_scale::AbstractVector, jac::CuMatrix, nlp_scaling_max_gradient)
    # Compute reduction on the GPU with built-in CUDA.jl function
    d_con_scale = maximum(abs, jac, dims=2)
    copyto!(con_scale, d_con_scale)
    con_scale .= min.(1.0, nlp_scaling_max_gradient ./ con_scale)
end

@kernel function _treat_fixed_variable_kernell!(dest, ind_fixed)
    k, j = @index(Global, NTuple)
    i = ind_fixed[k]

    if i == j
        dest[i, i] = 1.0
    else
        dest[i, j] = 0.0
        dest[j, i] = 0.0
    end
end

function MadNLP.treat_fixed_variable!(kkt::MadNLP.AbstractKKTSystem{T, MT}) where {T, MT<:CuMatrix{T}}
    length(kkt.ind_fixed) == 0 && return
    aug = kkt.aug_com
    d_ind_fixed = kkt.ind_fixed |> CuVector # TODO: allocate ind_fixed directly on the GPU
    ndrange = (length(d_ind_fixed), size(aug, 1))
    ev = _treat_fixed_variable_kernell!(CUDADevice())(aug, d_ind_fixed, ndrange=ndrange)
    wait(ev)
end

#=
    DenseKKTSystem kernels
=#
function MadNLP.mul!(y::AbstractVector, kkt::MadNLP.DenseKKTSystem{T, VT, MT}, x::AbstractVector) where {T, VT<:CuVector{T}, MT<:CuMatrix{T}}
    # Load buffers
    haskey(kkt.etc, :hess_w1) || (kkt.etc[:hess_w1] = CuVector{T}(undef, size(kkt.aug_com, 1)))
    haskey(kkt.etc, :hess_w2) || (kkt.etc[:hess_w2] = CuVector{T}(undef, size(kkt.aug_com, 1)))

    d_x = kkt.etc[:hess_w1]::VT
    d_y = kkt.etc[:hess_w2]::VT

    # x and y can be host arrays. Copy them on the device to avoid side effect.
    copyto!(d_x, x)
    LinearAlgebra.mul!(d_y, kkt.aug_com, d_x)
    copyto!(y, d_y)
end

function MadNLP.jtprod!(y::AbstractVector, kkt::MadNLP.DenseKKTSystem{T, VT, MT}, x::AbstractVector) where {T, VT<:CuVector{T}, MT<:CuMatrix{T}}
    # Load buffers
    haskey(kkt.etc, :jac_w1) || (kkt.etc[:jac_w1] = CuVector{T}(undef, size(kkt.jac, 1)))
    haskey(kkt.etc, :jac_w2) || (kkt.etc[:jac_w2] = CuVector{T}(undef, size(kkt.jac, 2)))

    d_x = kkt.etc[:jac_w1]::VT
    d_y = kkt.etc[:jac_w2]::VT

    # x and y can be host arrays. Copy them on the device to avoid side effect.
    copyto!(d_x, x)
    LinearAlgebra.mul!(d_y, kkt.jac', d_x)
    copyto!(y, d_y)
end

function MadNLP.set_aug_diagonal!(kkt::MadNLP.DenseKKTSystem{T, VT, MT}, ips::MadNLP.InteriorPointSolver) where {T, VT<:CuVector{T}, MT<:CuMatrix{T}}
    haskey(kkt.etc, :pr_diag_host) || (kkt.etc[:pr_diag_host] = Vector{T}(undef, length(kkt.pr_diag)))
    pr_diag_h = kkt.etc[:pr_diag_host]::Vector{T}
    # Broadcast is not working as MadNLP array are allocated on the CPU,
    # whereas pr_diag is allocated on the GPU
    pr_diag_h .= ips.zl./(ips.x.-ips.xl) .+ ips.zu./(ips.xu.-ips.x)
    copyto!(kkt.pr_diag, pr_diag_h)
    fill!(kkt.du_diag, 0.0)
end

@kernel function _build_dense_kkt_system_kernel!(
    dest, hess, jac, pr_diag, du_diag, diag_hess, n, m, ns
)
    i, j = @index(Global, NTuple)
    if (i <= n)
        # Transfer Hessian
        if (i == j)
            dest[i, i] = pr_diag[i] + diag_hess[i]
        elseif j <= n
            dest[i, j] = hess[i, j]
            dest[j, i] = hess[j, i]
        end
    elseif i <= n + ns
        # Transfer slack diagonal
        dest[i, i] = pr_diag[i]
    elseif i <= n + ns + m
        # Transfer Jacobian
        i_ = i - n - ns
        dest[i, j] = jac[i_, j]
        dest[j, i] = jac[i_, j]
        # Transfer dual regularization
        dest[i, i] = du_diag[i_]
    end
end

function MadNLP._build_dense_kkt_system!(
    dest::CuMatrix, hess::CuMatrix, jac::CuMatrix,
    pr_diag::CuVector, du_diag::CuVector, diag_hess::CuVector, n, m, ns
)
    ndrange = (n+m+ns, n+ns)
    ev = _build_dense_kkt_system_kernel!(CUDADevice())(dest, hess, jac, pr_diag, du_diag, diag_hess, n, m, ns, ndrange=ndrange)
    wait(ev)
end

function MadNLP.compress_jacobian!(kkt::MadNLP.DenseKKTSystem{T, VT, MT}) where {T, VT<:CuVector{T}, MT<:CuMatrix{T}}
    m = size(kkt.jac, 1)
    n = size(kkt.hess, 1)
    # Extract diagonal terms corresponding to inequalities
    index = (LinearAlgebra.diagind(kkt.jac) .+ n * m)[kkt.ind_ineq]
    # Add slack indexes
    kkt.jac[index] .= -one(T)
    # Scale
    kkt.jac .*= kkt.jacobian_scaling
    return
end


text:
Base.@kwdef mutable struct SIPPSNode
    v::Int
    low::Int
    high::Int
    id::Int
    is_goal::Bool = false
    f::Float64
    c::Int
    parent::Union{Nothing,SIPPSNode} = nothing
end

const SafeIntervalTable = Vector{Vector{@NamedTuple{low::Int, high::Int}}}

function same_identity(n1::SIPPSNode, n2::SIPPSNode)
    return n1.v == n2.v && n1.id == n2.id && n1.is_goal == n2.is_goal
end

function build_safe_interval_table(g::AbstractGraph, soft_obstacles::Reservation)
    T = 1
    for (t, v) in soft_obstacles
        T = max(t, T)
    end
    safe_interval_table = [@NamedTuple{low::Int, high::Int}[] for v in 1:nv(g)]
    for v in 1:nv(g)
        t = 1
        while t <= T
            low = t
            status = (t, v) in soft_obstacles
            while t <= T && status == ((t + 1, v) in soft_obstacles)
                t += 1
            end
            push!(safe_interval_table[v], (low=low, high=t+1))
            t += 1
        end
        (; low, high) = safe_interval_table[v][end]
        safe_interval_table[v][end] = (low=low, high=typemax(Int) ÷ 100)
    end
    return safe_interval_table
end

function intervals_intersect((a1, b1)::Tuple{Int,Int}, (a2, b2)::Tuple{Int,Int})
    return (a2 <= a1 < b2) || (a2 <= b1 < b2)
end

function has_soft_obstacles(safe_interval_table::SafeIntervalTable, v::Integer, id::Integer)
    l = length(safe_interval_table[v])
    return (l - id) % 2 == 0  # the last interval is always obstacle-free
end

function extract_path(n::SIPPSNode)
    path = Path()
    while !isnothing(n)
        pushfirst!(path, (n.low, n.v))
        n = n.parent
    end
    return path
end

function compare_with_identical!(X::Vector{SIPPSNode}, n::SIPPSNode)
    indices_to_delete = Int[]
    n_v, n_low, n_high, n_id, n_is_goal, n_c = n.v, n.low, n.high, n.id, n.is_goal, n.c
    for i in eachindex(X)
        q = X[i]
        if n_v == q.v && n_id == q.id && n_is_goal == q.is_goal
            q_low, q_high, q_c = q.low, q.high, q.c
            if q_low <= n_low && q_c <= n_c
                return false, indices_to_delete
            elseif n_low <= q_low && n_c <= q_c
                push!(indices_to_delete, i)
            elseif n_low < q_high && q_low < n_high
                if n_low < q_low
                    n.high = q_low
                else
                    q.high = n_low
                end
            end
        end
    end
    return true, indices_to_delete
end

function insert_node!(Q::PriorityQueue{SIPPSNode}, P::Vector{SIPPSNode}, n::SIPPSNode)
    Q_keep_n, Q_indices_to_delete = compare_with_identical!(keys(Q), n)
    P_keep_n, P_indices_to_delete = compare_with_identical!(P, n)
    deleteat!(Q, Q_indices_to_delete)
    deleteat!(P, P_indices_to_delete)
    if P_keep_n && Q_keep_n
        enqueue!(Q, n, (n.c, n.f))
        return true
    end
end

function expand_node!(
    Q::PriorityQueue{SIPPSNode},
    P::Vector{SIPPSNode},
    n::SIPPSNode,
    g::AbstractGraph,
    heuristic,
    safe_interval_table::SafeIntervalTable,
)
    reachable = Tuple{Int,Int}[]
    for v in outneighbors(g, n.v)
        for id in 1:length(safe_interval_table[v])
            (; low, high) = safe_interval_table[v][id]
            if intervals_intersect((low, high), (n.low + 1, n.high + 1))
                push!(reachable, (v, id))
            end
        end
    end
    for (v, id) in reachable
        (; low, high) = safe_interval_table[v][id]
        n3 = SIPPSNode(;
            v=v,
            low=low,
            high=high,
            id=id,
            is_goal=false,
            f=low + heuristic(v),
            c=n.c + Int(has_soft_obstacles(safe_interval_table, v, id)),
            parent=n,
        )
        insert_node!(Q, P, n3)
    end
end

function SIPPS(
    g::AbstractGraph,
    s::Integer,
    d::Integer,
    t0::Integer;
    heuristic=v -> 0,
    soft_obstacles::Reservation,
)
    safe_interval_table = build_safe_interval_table(g, soft_obstacles)
    root_id = minimum(
        k for (k, (; low, high)) in enumerate(safe_interval_table[s]) if low <= t0 <= high
    )
    root = SIPPSNode(;
        v=s,
        low=t0,
        high=safe_interval_table[s][root_id].high,
        id=root_id,
        is_goal=false,
        f=safe_interval_table[s][root_id].low + heuristic(s),
        c=Int(has_soft_obstacles(safe_interval_table, s, root_id)),
        parent=nothing,
    )

    Q = PriorityQueue{SIPPSNode,Tuple{Int,Float64}}()
    P = SIPPSNode[]
    enqueue!(Q, root, (root.c, root.f))

    while !isempty(Q)
        n = dequeue!(Q)
        if n.is_goal
            return extract_path(n)
        elseif n.v == d
            c_future = count(t > n.low for (t, v) in soft_obstacles if v == d)
            if c_future == 0
                return extract_path(n)
            end
            n2 = SIPPSNode(;
                v=n.v,
                low=n.low,
                high=n.high,
                id=n.id,
                is_goal=true,
                f=n.f,
                c=n.c + c_future,
                parent=n.parent,
            )
            insert_node!(Q, P, n2)
        end
        expand_node!(Q, P, n, g, heuristic, safe_interval_table)
        if !(n in P)
            push!(P, n)
        end
    end

    return Path()
end

function cooperative_SIPPS!(solution::Solution, mapf::MAPF; agents=1:nb_agents(mapf))
    fixed_agents = [a for a = 1:nb_agents(mapf) if !(a in agents)]
    reservation = compute_reservation(solution, mapf; agents=fixed_agents)
    for a in agents
        s, d, t0 = mapf.sources[a], mapf.destinations[a], mapf.starting_times[a]
        dist = mapf.distances_to_destinations[d]
        heuristic(v) = dist[v]
        path = SIPPS(
            mapf.graph, s, d, t0; heuristic=heuristic, soft_obstacles=reservation
        )
        solution[a] = path
        update_reservation(reservation, path, mapf)
    end
end


text:
#-----------------------------------------------------------
#-----------------------------------------------------------
# Define useful windows types
# -- part of WinReg Module --
#
#File: wintypes.jl
#Author: Leonardo Cocco
#Creation Date: 11-04-2016
#Last update 22-04-2016
#ver: 0.10
#----------------------------------------------------------
#----------------------------------------------------------

typealias POINTER UInt
typealias LPVOID Ptr{UInt8}
typealias WINBOOL UInt
typealias BOOL UInt
typealias BYTE UInt8
typealias LONG Int32
typealias DWORD UInt32
typealias QWORD UInt64
typealias HANDLE Ptr{Void}
typealias LPDWORD Ptr{DWORD}
typealias LPBYTE Ptr{BYTE}
typealias LPCWSTR Ptr{UInt16}
typealias LPCSTR Ptr{UInt8}
#typealias LPCTSTR Ptr{UInt8}

type FILETIME
	dwLowDateTime::DWORD
	dwHighDateTime::DWORD
end
typealias PFILETIME Ptr{FILETIME}

type SECURITY_ATTRIBUTES
	nLength::DWORD
    lpSecurityDescriptor::LPVOID
    bInheritHandle::WINBOOL
end
typealias LPSECURITY_ATTRIBUTES Ptr{SECURITY_ATTRIBUTES}

const ERROR_SUCCESS = 0x00
const ERROR_MORE_DATA = 234
const ERROR_NO_MORE_ITEMS = 259

baremodule SAM
	const DELETE = 0x00010000
	const READ_CONTROL = 0x00020000
	const WRITE_DAC = 0x00040000
	const WRITE_OWNER = 0x00080000
	const SYNCHRONIZE = 0x00100000

	const STANDARD_RIGHTS_REQUIRED = 0x000F0000

	const STANDARD_RIGHTS_READ = (READ_CONTROL)
	const STANDARD_RIGHTS_WRITE = (READ_CONTROL)
	const STANDARD_RIGHTS_EXECUTE = (READ_CONTROL)

	const STANDARD_RIGHTS_ALL = 0x001F0000
	const SPECIFIC_RIGHTS_ALL = 0x0000FFFF
end