text:
#using Base.Threads
import StatsBase: BaseRadixSortSafeTypes
import Base: Forward, ForwardOrdering, Reverse, ReverseOrdering, sortperm, Ordering,
            setindex!, getindex, similar, sort!, Reverse

import Base.Threads: nthreads, threadid

# create bits types for easy loading of bytes of lengths up to 15
primitive type Bits24 24 end
primitive type Bits40 40 end
primitive type Bits48 48 end
primitive type Bits56 56 end
primitive type Bits72 72 end
primitive type Bits80 80 end
primitive type Bits88 88 end
primitive type Bits96 96 end
primitive type Bits104 104 end
primitive type Bits112 112 end
primitive type Bits120 120 end

# sort it by using sorttwo! on the pointer
function sort_spointer!(svec::Vector{String}, lo::Int, hi::Int, o::O) where O <: Union{ForwardOrdering, ReverseOrdering}
    if lo >= hi;  return svec;  end
    # the length subarray to sort
    l = hi - lo + 1

    # find the maximum string length
    lens = maximum(sizeof, svec)
    skipbytes = lens

    if lens > 0
        while lens > 4
            skipbytes = max(0, skipbytes - 8)
            bits64 = zeros(UInt64, l)
            if o == Reverse
                bits64[lo:hi] .= .~load_bits.(UInt64, @view(svec[lo:hi]), skipbytes)
                SortingLab.sorttwo!(bits64, svec, lo, hi)
            else
                bits64[lo:hi] .= load_bits.(UInt64, @view(svec[lo:hi]), skipbytes)
                SortingLab.sorttwo!(bits64, svec, lo, hi)
            end
            lens -= 8
        end
        if lens > 0
            skipbytes = max(0, skipbytes - 4)
            bits32 = zeros(UInt32, l)
            if o == Reverse
                bits32[lo:hi] .= .~load_bits.(UInt32, @view(svec[lo:hi]), skipbytes)
                SortingLab.sorttwo!(bits32, svec, lo, hi)
            else
                bits32[lo:hi] .= load_bits.(UInt32, @view(svec[lo:hi]), skipbytes)
                SortingLab.sorttwo!(bits32, svec, lo, hi)
            end
            lens -= 4
        end
    end
    # unsafe_pointer_to_objref.(ptrs - 8)
    svec
end

# loads `remaining_bytes_to_load` bytes from `ptrs` which is a C-style pointer to a string
# these functions assumes that remaining_bytes_to_load > 0
function load_bits_with_padding(::Type{UInt128}, ptrs::Ptr{UInt8}, remaining_bytes_to_load)::UInt128
    nbits_to_shift_away = 8(sizeof(UInt128) - remaining_bytes_to_load)

    # the below checks if the string is less than 16 bytes away from the page
    # boundary assuming the page size is 4kb
    # see https://discourse.julialang.org/t/is-there-a-way-to-check-how-far-away-a-pointer-is-from-a-page-boundary/8147/11?u=xiaodai
    if (UInt(ptrs) & 0xfff) > 0xff0
        if  remaining_bytes_to_load == 15
            return ntoh(Base.zext_int(UInt128, unsafe_load(Ptr{Bits120}(ptrs))))
        elseif  remaining_bytes_to_load == 14
            return ntoh(Base.zext_int(UInt128, unsafe_load(Ptr{Bits112}(ptrs))))
        elseif  remaining_bytes_to_load == 13
            return ntoh(Base.zext_int(UInt128, unsafe_load(Ptr{Bits104}(ptrs))))
        elseif  remaining_bytes_to_load == 12
            return ntoh(Base.zext_int(UInt128, unsafe_load(Ptr{Bits96}(ptrs))))
        elseif  remaining_bytes_to_load == 11
            return ntoh(Base.zext_int(UInt128, unsafe_load(Ptr{Bits88}(ptrs))))
        elseif  remaining_bytes_to_load == 10
            return ntoh(Base.zext_int(UInt128, unsafe_load(Ptr{Bits80}(ptrs))))
        elseif  remaining_bytes_to_load == 9
            return ntoh(Base.zext_int(UInt128, unsafe_load(Ptr{Bits72}(ptrs))))
        elseif  remaining_bytes_to_load == 8
            return ntoh(Base.zext_int(UInt128, unsafe_load(Ptr{UInt64}(ptrs))))
        elseif  remaining_bytes_to_load == 7
            return ntoh(Base.zext_int(UInt128, unsafe_load(Ptr{Bits56}(ptrs))))
        elseif  remaining_bytes_to_load == 6
            return ntoh(Base.zext_int(UInt128, unsafe_load(Ptr{Bits48}(ptrs))))
        elseif  remaining_bytes_to_load == 5
            return ntoh(Base.zext_int(UInt128, unsafe_load(Ptr{Bits40}(ptrs))))
        elseif  remaining_bytes_to_load == 4
            return ntoh(Base.zext_int(UInt128, unsafe_load(Ptr{UInt32}(ptrs))))
        elseif  remaining_bytes_to_load == 3
            return ntoh(Base.zext_int(UInt128, unsafe_load(Ptr{Bits24}(ptrs))))
        elseif  remaining_bytes_to_load == 2
            return ntoh(Base.zext_int(UInt128, unsafe_load(Ptr{UInt16}(ptrs))))
        else
            return ntoh(Base.zext_int(UInt128, unsafe_load(Ptr{UInt8}(ptrs))))
        end
    else
        return ntoh(unsafe_load(Ptr{UInt128}(ptrs))) >> nbits_to_shift_away << nbits_to_shift_away
    end
end

function load_bits_with_padding(::Type{UInt64}, ptrs::Ptr{UInt8}, remaining_bytes_to_load)::UInt64
    nbits_to_shift_away = 8(sizeof(UInt64) - remaining_bytes_to_load)

    # the below checks if the string is less than 8 bytes away from the page
    # boundary assuming the page size is 4kb
    # see https://discourse.julialang.org/t/is-there-a-way-to-check-how-far-away-a-pointer-is-from-a-page-boundary/8147/11?u=xiaodai
    if (UInt(ptrs) & 0xfff) > 0xff8
        if  remaining_bytes_to_load == 7
            return ntoh(Base.zext_int(UInt64, unsafe_load(Ptr{Bits56}(ptrs))))
        elseif  remaining_bytes_to_load == 6
            return ntoh(Base.zext_int(UInt64, unsafe_load(Ptr{Bits48}(ptrs))))
        elseif  remaining_bytes_to_load == 5
            return ntoh(Base.zext_int(UInt64, unsafe_load(Ptr{Bits40}(ptrs))))
        elseif  remaining_bytes_to_load == 4
            return ntoh(Base.zext_int(UInt64, unsafe_load(Ptr{UInt32}(ptrs))))
        elseif  remaining_bytes_to_load == 3
            return ntoh(Base.zext_int(UInt64, unsafe_load(Ptr{Bits24}(ptrs))))
        elseif  remaining_bytes_to_load == 2
            return ntoh(Base.zext_int(UInt64, unsafe_load(Ptr{UInt16}(ptrs))))
        else
            return ntoh(Base.zext_int(UInt64, unsafe_load(Ptr{UInt8}(ptrs))))
        end
    else
        return ntoh(unsafe_load(Ptr{UInt64}(ptrs))) >> nbits_to_shift_away << nbits_to_shift_away
    end
end

function load_bits_with_padding(::Type{UInt32}, ptrs::Ptr{UInt8}, remaining_bytes_to_load)::UInt32
    nbits_to_shift_away = 8(sizeof(UInt32) - remaining_bytes_to_load)

    # the below checks if the string is less than 4 bytes away from the page
    # boundary assuming the page size is 4kb
    # see https://discourse.julialang.org/t/is-there-a-way-to-check-how-far-away-a-pointer-is-from-a-page-boundary/8147/11?u=xiaodai
    if (UInt(ptrs) & 0xfff) > 0xffc
        if  remaining_bytes_to_load == 3
            return ntoh(Base.zext_int(UInt32, unsafe_load(Ptr{Bits24}(ptrs))))
        elseif  remaining_bytes_to_load == 2
            return ntoh(Base.zext_int(UInt32, unsafe_load(Ptr{UInt16}(ptrs))))
        else
            return ntoh(Base.zext_int(UInt32, unsafe_load(Ptr{UInt8}(ptrs))))
        end
    else
        return ntoh(unsafe_load(Ptr{UInt32}(ptrs))) >> nbits_to_shift_away << nbits_to_shift_away
    end
end

"""
    load_bits([type,] s, skipbytes)

Load the underlying bits of a string `s` into a `type` of the user's choosing.
The default is `UInt`, so on a 64 bit machine it loads 64 bits (8 bytes) at a time.
If the `String` is shorter than 8 bytes then it's padded with 0.

- `type`:       any bits type that has `>>`, `<<`, and `&` operations defined
- `s`:          a `String`
- `skipbytes`:  how many bytes to skip e.g. load_bits("abc", 1) will load "bc" as bits
"""
# Some part of the return result should be padded with 0s.
# To prevent any possibility of segfault we load the bits using
# successively smaller types
# it is assumed that the type you are trying to load into needs padding
# i.e. `remaining_bytes_to_load > 0`
function load_bits(::Type{T}, s::String, skipbytes = 0)::T where T <: Unsigned
    n = sizeof(s)
    if n < skipbytes
        res = zero(T)
    elseif n - skipbytes >= sizeof(T)
        res = ntoh(unsafe_load(Ptr{T}(pointer(s, skipbytes+1))))
    else
        ptrs  = pointer(s) + skipbytes
        remaining_bytes_to_load = n - skipbytes
        res = load_bits_with_padding(T, ptrs, remaining_bytes_to_load)
    end
    return res
end

# Radix sort for strings
"""
    radixsort!(svec)

Applies radix sort to the string vector, svec, and sort it in place.
"""
function radixsort!(svec::AbstractVector{String}, lo::Int, hi::Int, o::O; RADIX_SIZE = 16, RADIX_MASK = 0xffff) where O <: Union{ForwardOrdering, ReverseOrdering}
    if lo >= hi;  return svec;  end
    # the length subarray to sort
    l = hi - lo + 1

    # find the maximum string length
    lens = maximum(sizeof, svec)
    skipbytes = lens
    if lens > 0
        # while lens > 8
        #     skipbytes = max(0, skipbytes - 16)
        #     bits128 = zeros(UInt128, l)
        #     if o == Reverse
        #         bits128[lo:hi] .= .~load_bits.(UInt128, @view(svec[lo:hi]), skipbytes)
        #         sorttwo!(bits128, svec, lo, hi)
        #     else
        #         bits128[lo:hi] .= load_bits.(UInt128, @view(svec[lo:hi]), skipbytes)
        #         sorttwo!(bits128, svec, lo, hi)
        #     end
        #     lens -= 16
        # end
        while lens > 4
            skipbytes = max(0, skipbytes - 8)
            bits64 = zeros(UInt64, l)
            if o == Reverse
                bits64[lo:hi] .= .~load_bits.(UInt64, @view(svec[lo:hi]), skipbytes)
                sorttwo!(bits64, svec, lo, hi, RADIX_SIZE, RADIX_MASK)
            else
                bits64[lo:hi] .= load_bits.(UInt64, @view(svec[lo:hi]), skipbytes)
                sorttwo!(bits64, svec, lo, hi, RADIX_SIZE, RADIX_MASK)
            end
            lens -= 8
        end
        if lens > 0
            skipbytes = max(0, skipbytes - 4)
            bits32 = zeros(UInt32, l)
            if o == Reverse
                bits32[lo:hi] .= .~load_bits.(UInt32, @view(svec[lo:hi]), skipbytes)
                sorttwo!(bits32, svec, lo, hi, RADIX_SIZE, RADIX_MASK)
            else
                bits32[lo:hi] .= load_bits.(UInt32, @view(svec[lo:hi]), skipbytes)
                sorttwo!(bits32, svec, lo, hi, RADIX_SIZE, RADIX_MASK)
            end
            lens -= 4
        end
    end
    svec
end

"""
    radixsort(svec, rev = false)

Applies radix sort to the string vector, svec, and sort it in place.
"""
fsort(svec::Vector{String}, rev = false, radix_opts = (16, 0xffff)) = radixsort(svec, rev , (16, 0xffff))
fsort!(svec::Vector{String}, rev = false, radix_opts = (16, 0xffff)) = radixsort!(svec, rev , (16, 0xffff))
radixsort(svec::Vector{String}, rev = false, radix_opts = (16, 0xffff)) = radixsort!(copy(svec), rev, radix_opts)
radixsort!(svec::Vector{String}, rev = false, radix_opts = (16, 0xffff)) =  radixsort!(svec, 1, length(svec), rev ? Base.Reverse : Base.Forward, RADIX_SIZE = radix_opts[1], RADIX_MASK = radix_opts[2])


text:
#=
Conservative 7-color palette from Points of view: Color blindness, Bang Wong - Nature Methods
https://www.nature.com/articles/nmeth.1618?WT.ec_id=NMETH-201106
=#
const default_colors = [
    RGB(230/255, 159/255, 0/255),
    RGB(86/255, 180/255, 233/255),
    RGB(0/255, 158/255, 115/255),
    RGB(240/255, 228/255, 66/255),
    RGB(0/255, 114/255, 178/255),
    RGB(213/255, 94/255, 0/255),
    RGB(204/255, 121/255, 167/255),
]

function default_palettes()
    return arguments(
        color=default_colors,
        marker=[:circle, :xcross, :utriangle, :diamond, :dtriangle, :star8, :pentagon, :rect],
        linestyle=[:solid, :dash, :dot, :dashdot, :dashdotdot],
        side=[:left, :right],
        layout=wrap,
    )
end

apply_palette(p::AbstractVector, uv) = [cycle(p, idx) for idx in eachindex(uv)]
apply_palette(::Automatic, uv) = eachindex(uv)
apply_palette(p, uv) = map(p, eachindex(uv))

# TODO: add more customizations?
struct Wrap end

const wrap = Wrap()

function apply_palette(::Wrap, uv)
    ncols = ceil(Int, sqrt(length(uv)))
    return [fldmod1(idx, ncols) for idx in eachindex(uv)]
end

struct ContinuousScale{T, F}
    f::F
    extrema::Tuple{T, T}
end

rescale(values, c::ContinuousScale) = values # Is this ideal?

struct CategoricalScale{S, T}
    data::S
    plot::T
end

function rescale(values, c::CategoricalScale)
    idxs = indexin(values, c.data)
    return c.plot[idxs]
end

Base.length(c::CategoricalScale) = length(c.data)

function default_scale(summary, palette)
    iscont = summary isa Tuple
    return if iscont
        f = palette isa Function ? palette : identity
        ContinuousScale(f, summary)
    else
        data = sort!(collect(summary))
        plot = apply_palette(palette, data)
        return CategoricalScale(data, plot)
    end
end

# Logic to infer good scales
function default_scales(summaries, palettes)
    palettes = merge!(map(_ -> automatic, summaries), palettes)
    return map(default_scale, summaries, palettes)
end

ext: jl


text:
# EXCLUDE FROM TESTING

using BenchmarkTools
BenchmarkTools.DEFAULT_PARAMETERS.seconds = 10
BenchmarkTools.DEFAULT_PARAMETERS.gcsample = true

include("reduce.jl")

CUDAnative.initialize()
const dev = device()
const cap = capability(dev)
@assert(cap >= v"3.0", "this example requires a newer GPU")

len = 10^7
input = ones(Int32, len)

## CPU

benchmark_cpu = @benchmarkable begin
        reduce(+, input)
    end

@show run(benchmark_cpu)

## CUDAnative

# PTX generation
open(joinpath(@__DIR__, "reduce.jl.ptx"), "w") do f
    CUDAnative.code_ptx(f, reduce_grid, Tuple{typeof(+), CuDeviceVector{Int32,AS.Global},
                                              CuDeviceVector{Int32,AS.Global}, Int32};
                        cap=v"6.1.0")
end

benchmark_gpu = @benchmarkable begin
        gpu_reduce(+, gpu_input, gpu_output)
        val = Array(gpu_output)[1]
    end setup=(
        val = nothing;
        gpu_input = CuArray($input);
        gpu_output = similar(gpu_input)
    ) teardown=(
        gpu_input = nothing;
        gpu_output = nothing
    )

@show run(benchmark_gpu)

## CUDA

using CUDAapi
using Libdl

cd(@__DIR__) do
    toolkit = CUDAapi.find_toolkit()
    nvcc = CUDAapi.find_cuda_binary("nvcc", toolkit)
    toolchain = CUDAapi.find_toolchain(toolkit)
    flags = `-ccbin=$(toolchain.host_compiler) -arch=sm_$(cap.major)$(cap.minor)`
    run(`$nvcc $flags -ptx -o reduce.cu.ptx reduce.cu`)
    run(`$nvcc $flags -shared --compiler-options '-fPIC' -o reduce.so reduce.cu`)
end

# Entry-point wrappers
lib = Libdl.dlopen(joinpath(@__DIR__, "reduce.so"))
setup_cuda(input)    = ccall(Libdl.dlsym(lib, "setup"), Ptr{Cvoid},
                             (Ptr{Cint}, Csize_t), input, length(input))
run_cuda(state)      = ccall(Libdl.dlsym(lib, "run"), Cint,
                             (Ptr{Cvoid},), state)
teardown_cuda(state) = ccall(Libdl.dlsym(lib, "teardown"), Cvoid,
                             (Ptr{Cvoid},), state)

# Correctness check (not part of verify.jl which is meant to run during testing)
using Test
let
    cuda_state = setup_cuda(input)
    cuda_val = run_cuda(cuda_state)
    teardown_cuda(cuda_state)
    @assert cuda_val == reduce(+, input)
end

benchmark_cuda = @benchmarkable begin
        val = run_cuda(state)
    end setup=(
        val = nothing;
        state = setup_cuda($input);
    ) teardown=(
        teardown_cuda(state)
    )

@show run(benchmark_cuda)


text:
using PyPlot

# assumption: seller blue #00A2FF, buyer orange #FF8C00
# we have already run the simulation_at_exp script -
plot(data_array_at_exp[:,1], data_array_at_exp[:,2], lw=2, color="#FF8C00")
ax = gca()
fig = gcf()
ax.set_facecolor("black")
fig.set_facecolor("black")
grid(color="dimgray", linestyle="dashed")

ax.spines["bottom"].set_color("white")
ax.spines["top"].set_color("white")
ax.spines["right"].set_color("white")
ax.spines["left"].set_color("white")
ax.tick_params(axis="x", colors="white")
ax.tick_params(axis="y", colors="white")

xlabel("Stock Price (USD)", fontsize=16,fontname="Arial",color="white")
ylabel("Profit or Loss (USD)", fontsize=16,fontname="Arial",color="white")


text:
g = evolving_graph(Char, Int, is_directed = false)
add_edge!(g, 'a', 'b', 1)
add_edge!(g, 'b', 'g', 1)
add_edge!(g, 'b', 'f', 1)
add_edge!(g, 'c', 'e', 2)
add_edge!(g, 'e', 'g', 2)
add_edge!(g, 'a', 'b', 2)
add_edge!(g, 'b', 'f', 2)
add_edge!(g, 'f', 'd', 2)
add_edge!(g, 'a', 'b', 3)
add_edge!(g, 'c', 'f', 3)
add_edge!(g, 'e', 'g', 3)
add_edge!(g, 'g', 'f', 3)

p = shortest_temporal_path(g, 'a', 1, 'e', 3)
display(p)

@test shortest_temporal_distance(g, 'c', 2, 'a', 2) == Inf

g = int_evolving_graph(6, 3)
add_edge!(g, 1, 2, 1)
add_edge!(g, 1, 3, 2)
add_edge!(g, 4, 5, 2)
add_edge!(g, 2, 3, 3)
add_edge!(g, 5, 6, 3)

p = shortest_temporal_path(g, 1, 1, 3, 3)
@test length(p) - 1 == shortest_temporal_distance(g, 1, 1, 3, 3)
display(p)


text:
#####
##### `Symmetric`/`Hermitian`
#####

function frule((_, ΔA, _), T::Type{<:LinearAlgebra.HermOrSym}, A::AbstractMatrix, uplo)
    return T(A, uplo), T(ΔA, uplo)
end

function rrule(T::Type{<:LinearAlgebra.HermOrSym}, A::AbstractMatrix, uplo)
    Ω = T(A, uplo)
    @inline function HermOrSym_pullback(ΔΩ)
        return (NoTangent(), _symherm_back(typeof(Ω), ΔΩ, uplo), NoTangent())
    end
    return Ω, HermOrSym_pullback
end

function frule((_, ΔA), TM::Type{<:Matrix}, A::LinearAlgebra.HermOrSym)
    return TM(A), TM(_symherm_forward(A, ΔA))
end
function frule((_, ΔA), ::Type{Array}, A::LinearAlgebra.HermOrSym)
    return Array(A), Array(_symherm_forward(A, ΔA))
end

function rrule(TM::Type{<:Matrix}, A::LinearAlgebra.HermOrSym)
    function Matrix_pullback(Ω̄)
        ΔΩ = unthunk(Ω̄)
        TA = _symhermtype(A)
        T∂A = TA{eltype(ΔΩ),typeof(ΔΩ)}
        uplo = A.uplo
        ∂A = T∂A(_symherm_back(typeof(A), ΔΩ, Symbol(uplo)), uplo)
        return NoTangent(), ∂A
    end
    return TM(A), Matrix_pullback
end
rrule(::Type{Array}, A::LinearAlgebra.HermOrSym) = rrule(Matrix, A)

# for Ω = Matrix(A::HermOrSym), push forward ΔA to get ∂Ω
function _symherm_forward(A, Ȧ)
    ΔA = unthunk(Ȧ)
    TA = _symhermtype(A)
    return if typeof(ΔA) <: TA
        ΔA
    else
        TA{eltype(ΔA),typeof(ΔA)}(ΔA, A.uplo)
    end
end

# for Ω = HermOrSym(A, uplo), pull back ΔΩ to get ∂A
@inline function _symherm_back(::Type{T}, Ω̄, uplo::Symbol) where {T}
    ΔΩ = unthunk(Ω̄)
    if T <: Symmetric
        return _symmetric_back(ΔΩ, uplo)
    elseif T <: Hermitian
        if ΔΩ isa AbstractMatrix{<:Real}
            return _symmetric_back(ΔΩ, uplo)
        else
            return _hermitian_back(ΔΩ, uplo)
        end
    end
    error()
end

@inline function _symmetric_back(Ω̄, uplo::Symbol)
    ΔΩ = unthunk(Ω̄)
    if ΔΩ isa Diagonal
        return ΔΩ
    elseif ΔΩ isa LinearAlgebra.AbstractTriangular
        if istriu(ΔΩ)
            return Matrix(uplo === :U ? ΔΩ : transpose(ΔΩ))
        else
            return Matrix(uplo === :U ? transpose(ΔΩ) : ΔΩ)
        end
    end
    L, U, D = LowerTriangular(ΔΩ), UpperTriangular(ΔΩ), Diagonal(ΔΩ)
    return Matrix(uplo === :U ? U .+ transpose(L) - D : L .+ transpose(U) - D)
end

@inline function _hermitian_back(ΔΩ, uplo::Symbol)
    if ΔΩ isa Diagonal
        return real.(ΔΩ)
    elseif ΔΩ isa LinearAlgebra.AbstractTriangular
        ∂UL = ΔΩ .- Diagonal(_extract_imag.(diag(ΔΩ)))
        if istriu(ΔΩ)
            return Matrix(uplo === :U ? ∂UL : ∂UL')
        else
            return Matrix(uplo === :U ? ∂UL' : ∂UL)
        end
    end
    L, U, rD = LowerTriangular(ΔΩ), UpperTriangular(ΔΩ), real.(Diagonal(ΔΩ))
    return Matrix(uplo === :U ? U .+ L' - rD : L .+ U' - rD)
end

#####
##### `eigen!`/`eigen`
#####

# rule is old but the usual references are
# real rules:
# Giles M. B., An extended collection of matrix derivative results for forward and reverse
# mode algorithmic differentiation.
# https://people.maths.ox.ac.uk/gilesm/files/NA-08-01.pdf.
# complex rules:
# Boeddeker C., Hanebrink P., et al, On the Computation of Complex-valued Gradients with
# Application to Statistically Optimum Beamforming. arXiv:1701.00392v2 [cs.NA]
#
# accounting for normalization convention appears in Boeddeker && Hanebrink.
# account for phase convention is unpublished.
function frule(
    (_, Ȧ),
    ::typeof(eigen!),
    A::LinearAlgebra.RealHermSymComplexHerm{<:BLAS.BlasReal,<:StridedMatrix};
    kwargs...,
)
    ΔA = unthunk(Ȧ)
    F = eigen!(A; kwargs...)
    ΔA isa AbstractZero && return F, ΔA
    λ, U = F.values, F.vectors
    tmp = U' * ΔA
    ∂K = mul!(ΔA.data, tmp, U)
    ∂Kdiag = @view ∂K[diagind(∂K)]
    ∂λ = real.(∂Kdiag)
    ∂K ./= λ' .- λ
    fill!(∂Kdiag, 0)
    ∂U = mul!(tmp, U, ∂K)
    _eigen_norm_phase_fwd!(∂U, A, U)
    ∂F = Tangent{typeof(F)}(values = ∂λ, vectors = ∂U)
    return F, ∂F
end

function rrule(
    ::typeof(eigen),
    A::LinearAlgebra.RealHermSymComplexHerm{<:BLAS.BlasReal,<:StridedMatrix};
    kwargs...,
)
    F = eigen(A; kwargs...)
    function eigen_pullback(ΔF::Tangent)
        λ, U = F.values, F.vectors
        Δλ, ΔU = ΔF.values, ΔF.vectors
        ΔU = ΔU isa AbstractZero ? ΔU : copy(ΔU)
        ∂A = eigen_rev!(A, λ, U, Δλ, ΔU)
        return NoTangent(), ∂A
    end
    eigen_pullback(ΔF::AbstractZero) = (NoTangent(), ΔF)
    return F, eigen_pullback
end

# ∂U is overwritten if not an `AbstractZero`
function eigen_rev!(A::LinearAlgebra.RealHermSymComplexHerm, λ, U, ∂λ, ∂U)
    ∂λ isa AbstractZero && ∂U isa AbstractZero && return ∂λ + ∂U
    Ā = similar(parent(A), eltype(U))
    tmp = ∂U
    if ∂U isa AbstractZero
        mul!(Ā, U, real.(∂λ) .* U')
    else
        _eigen_norm_phase_rev!(∂U, A, U)
        ∂K = mul!(Ā, U', ∂U)
        ∂K ./= λ' .- λ
        ∂K[diagind(∂K)] .= real.(∂λ)
        mul!(tmp, ∂K, U')
        mul!(Ā, U, tmp)
    end
    ∂A = _hermitrizelike!(Ā, A)
    return ∂A
end

# NOTE: for small vₖ, the derivative of sign(vₖ) explodes, causing the tangents to become
# unstable even for phase-invariant programs. So for small vₖ we don't account for the phase
# in the gradient. Then derivatives are accurate for phase-invariant programs but inaccurate
# for phase-dependent programs that have low vₖ.

_eigen_norm_phase_fwd!(∂V, ::Union{Symmetric{T,S},Hermitian{T,S}}, V) where {T<:Real,S} = ∂V
function _eigen_norm_phase_fwd!(∂V, A::Hermitian{<:Complex}, V)
    k = A.uplo === 'U' ? size(A, 1) : 1
    ϵ = sqrt(eps(real(eltype(V))))
    @inbounds for i in axes(V, 2)
        v = @view V[:, i]
        vₖ = real(v[k])
        if abs(vₖ) > ϵ
            ∂v = @view ∂V[:, i]
            ∂v .-= v .* (im * (imag(∂v[k]) / vₖ))
        end
    end
    return ∂V
end

_eigen_norm_phase_rev!(∂V, ::Union{Symmetric{T,S},Hermitian{T,S}}, V) where {T<:Real,S} = ∂V
function _eigen_norm_phase_rev!(∂V, A::Hermitian{<:Complex}, V)
    k = A.uplo === 'U' ? size(A, 1) : 1
    ϵ = sqrt(eps(real(eltype(V))))
    @inbounds for i in axes(V, 2)
        v = @view V[:, i]
        vₖ = real(v[k])
        if abs(vₖ) > ϵ
            ∂v = @view ∂V[:, i]
            ∂c = dot(v, ∂v)
            ∂v[k] -= im * (imag(∂c) / vₖ)
        end
    end
    return ∂V
end

#####
##### `eigvals!`/`eigvals`
#####

function frule(
    (_, Ȧ),
    ::typeof(eigvals!),
    A::LinearAlgebra.RealHermSymComplexHerm{<:BLAS.BlasReal,<:StridedMatrix};
    kwargs...,
)
    ΔA = unthunk(Ȧ)
    ΔA isa AbstractZero && return eigvals!(A; kwargs...), ΔA
    F = eigen!(A; kwargs...)
    λ, U = F.values, F.vectors
    tmp = ΔA * U
    # diag(U' * tmp) without computing matrix product
    ∂λ = similar(λ)
    @inbounds for i in eachindex(λ)
        ∂λ[i] = @views real(dot(U[:, i], tmp[:, i]))
    end
    return λ, ∂λ
end

function rrule(
    ::typeof(eigvals),
    A::LinearAlgebra.RealHermSymComplexHerm{<:BLAS.BlasReal,<:StridedMatrix};
    kwargs...,
)
    F, eigen_back = rrule(eigen, A; kwargs...)
    λ = F.values
    function eigvals_pullback(Δλ)
        ∂F = Tangent{typeof(F)}(values = Δλ)
        _, ∂A = eigen_back(∂F)
        return NoTangent(), ∂A
    end
    return λ, eigvals_pullback
end

#####
##### `svd`
#####

# NOTE: rrule defined because the `svd` primal mutates after calling `eigen`.
# otherwise, this rule just applies the chain rule and can be removed when mutation
# is supported by reverse-mode AD packages
function rrule(::typeof(svd), A::LinearAlgebra.RealHermSymComplexHerm{<:BLAS.BlasReal,<:StridedMatrix})
    F = svd(A)
    function svd_pullback(ΔF::Tangent)
        U, V = F.U, F.V
        c = _svd_eigvals_sign!(similar(F.S), U, V)
        λ = F.S .* c
        ∂λ = ΔF.S isa AbstractZero ? ΔF.S : ΔF.S .* c
        if all(x -> x isa AbstractZero, (ΔF.U, ΔF.V, ΔF.Vt))
            ∂U = ΔF.U + ΔF.V + ΔF.Vt
        else
            ∂U = ΔF.U .+ (ΔF.V .+ ΔF.Vt') .* c'
        end
        ∂A = eigen_rev!(A, λ, U, ∂λ, ∂U)
        return NoTangent(), ∂A
    end
    svd_pullback(ΔF::AbstractZero) = (NoTangent(), ΔF)
    return F, svd_pullback
end

# given singular vectors, compute sign of eigenvalues corresponding to singular values
function _svd_eigvals_sign!(c, U, V)
    n = size(U, 1)
    @inbounds broadcast!(c, eachindex(c)) do i
        u = @views U[:, i]
        # find element not close to zero
        # at least one element satisfies abs2(x) ≥ 1/n > 1/(n + 1)
        k = findfirst(x -> (n + 1) * abs2(x) ≥ 1, u)
        return sign(real(u[k]) * real(V[k, i]))
    end
    return c
end

#####
##### `svdvals`
#####

# NOTE: rrule defined because `svdvals` calls mutating `svdvals!` internally.
# can be removed when mutation is supported by reverse-mode AD packages
function rrule(::typeof(svdvals), A::LinearAlgebra.RealHermSymComplexHerm{<:BLAS.BlasReal,<:StridedMatrix})
    λ, back = rrule(eigvals, A)
    S = abs.(λ)
    p = sortperm(S; rev=true)
    permute!(S, p)
    function svdvals_pullback(ΔS)
        ∂λ = real.(ΔS)
        invpermute!(∂λ, p)
        ∂λ .*= sign.(λ)
        _, ∂A = back(∂λ)
        return NoTangent(), unthunk(∂A)
    end
    svdvals_pullback(ΔS::AbstractZero) = (NoTangent(), ΔS)
    return S, svdvals_pullback
end

#####
##### matrix functions
#####

# Formula for frule (Fréchet derivative) from Daleckiĭ-Kreĭn theorem given in Theorem 3.11 of
# Higham N.J. Functions of Matrices: Theory and Computation. 2008. ISBN: 978-0-898716-46-7.
# rrule is derived from frule. These rules are more stable for degenerate matrices than
# applying the chain rule to the rules for `eigen`.

for func in (:exp, :log, :sqrt, :cos, :sin, :tan, :cosh, :sinh, :tanh, :acos, :asin, :atan, :acosh, :asinh, :atanh)
    @eval begin
        function frule((_, Ȧ), ::typeof($func), A::LinearAlgebra.RealHermSymComplexHerm)
            ΔA = unthunk(Ȧ)
            Y, intermediates = _matfun($func, A)
            Ȳ = _matfun_frechet($func, ΔA, A, Y, intermediates)
            # If ΔA was hermitian, then ∂Y has the same structure as Y
            ∂Y = if ishermitian(ΔA) && (isa(Y, Symmetric) || isa(Y, Hermitian))
                _symhermlike!(Ȳ, Y)
            else
                Ȳ
            end
            return Y, ∂Y
        end

        function rrule(::typeof($func), A::LinearAlgebra.RealHermSymComplexHerm)
            Y, intermediates = _matfun($func, A)
            function $(Symbol(func, :_pullback))(ΔY)
                # for Hermitian Y, we don't need to realify the diagonal of ΔY, since the
                # effect is the same as applying _hermitrizelike! at the end
                ∂Y = eltype(Y) <: Real ? real(ΔY) : ΔY
                Ā = _matfun_frechet_adjoint($func, ∂Y, A, Y, intermediates)
                # the cotangent of Hermitian A should be Hermitian
                ∂A = _hermitrizelike!(Ā, A)
                return NoTangent(), ∂A
            end
            return Y, $(Symbol(func, :_pullback))
        end
    end
end

function frule((_, Ȧ), ::typeof(sincos), A::LinearAlgebra.RealHermSymComplexHerm)
    ΔA = unthunk(Ȧ)
    sinA, (λ, U, sinλ, cosλ) = _matfun(sin, A)
    cosA = _symhermtype(sinA)((U * Diagonal(cosλ)) * U')
    # We will overwrite tmp matrix several times to hold different values
    tmp = mul!(similar(U, Base.promote_eltype(ΔA, U)), ΔA, U)
    ∂Λ = mul!(similar(U), U', tmp)
    ∂sinΛ = _muldiffquotmat!!(similar(∂Λ), sin, λ, sinλ, cosλ, ∂Λ)
    ∂cosΛ = _muldiffquotmat!!(∂Λ, cos, λ, cosλ, -sinλ, ∂Λ)
    ∂sinA = _symhermlike!(mul!(∂sinΛ, U, mul!(tmp, ∂sinΛ, U')), sinA)
    ∂cosA = _symhermlike!(mul!(∂cosΛ, U, mul!(tmp, ∂cosΛ, U')), cosA)
    Y = (sinA, cosA)
    ∂Y = Tangent{typeof(Y)}(∂sinA, ∂cosA)
    return Y, ∂Y
end

function rrule(::typeof(sincos), A::LinearAlgebra.RealHermSymComplexHerm)
    sinA, (λ, U, sinλ, cosλ) = _matfun(sin, A)
    cosA = typeof(sinA)((U * Diagonal(cosλ)) * U', sinA.uplo)
    Y = (sinA, cosA)
    function sincos_pullback((ΔsinA, ΔcosA)::Tangent)
        ΔsinA isa AbstractZero && ΔcosA isa AbstractZero && return NoTangent(), ΔsinA + ΔcosA
        if eltype(A) <: Real
            ΔsinA, ΔcosA = real(ΔsinA), real(ΔcosA)
        end
        if ΔcosA isa AbstractZero
            Ā = _matfun_frechet_adjoint(sin, ΔsinA, A, sinA, (λ, U, sinλ, cosλ))
        elseif ΔsinA isa AbstractZero
            Ā = _matfun_frechet_adjoint(cos, ΔcosA, A, cosA, (λ, U, cosλ, -sinλ))
        else
            # we will overwrite tmp with various temporary values during this computation
            tmp = mul!(similar(U, Base.promote_eltype(U, ΔsinA, ΔcosA)), ΔsinA, U)
            ∂sinΛ = mul!(similar(tmp), U', tmp)
            ∂cosΛ = U' * mul!(tmp, ΔcosA, U)
            ∂Λ = _muldiffquotmat!!(∂sinΛ, sin, λ, sinλ, cosλ, ∂sinΛ)
            ∂Λ = _muldiffquotmat!!(∂Λ, cos, λ, cosλ, -sinλ, ∂cosΛ, true)
            Ā = mul!(∂Λ, U, mul!(tmp, ∂Λ, U'))
        end
        ∂A = _hermitrizelike!(Ā, A)
        return NoTangent(), ∂A
    end
    return Y, sincos_pullback
end

"""
    _matfun(f, A::LinearAlgebra.RealHermSymComplexHerm)

Compute the matrix function `f(A)` for real or complex hermitian `A`.
The function returns a tuple containing the result and a tuple of intermediates to be
reused by [`_matfun_frechet`](@ref) to compute the Fréchet derivative.

Note any function `f` used with this **must** have a `frule` defined on it.
"""
function _matfun(f, A::LinearAlgebra.RealHermSymComplexHerm)
    λ, U = eigen(A)
    if all(λi -> _isindomain(f, λi), λ)
        fλ_df_dλ = map(λi -> frule((ZeroTangent(), true), f, λi), λ)
    else  # promote to complex if necessary
        fλ_df_dλ = map(λi -> frule((ZeroTangent(), true), f, complex(λi)), λ)
    end
    fλ = first.(fλ_df_dλ)
    df_dλ = last.(unthunk.(fλ_df_dλ))
    fA = (U * Diagonal(fλ)) * U'
    Y = if eltype(A) <: Real
        Symmetric(fA)
    elseif eltype(fλ) <: Complex
        fA
    else
        Hermitian(fA)
    end
    intermediates = (λ, U, fλ, df_dλ)
    return Y, intermediates
end

# Computes ∂Y = U * (P .* (U' * ΔA * U)) * U' with fewer allocations
function _matfun_frechet(f, ΔA, A::LinearAlgebra.RealHermSymComplexHerm, Y, (λ, U, fλ, df_dλ))
    # We will overwrite tmp matrix several times to hold different values
    tmp = mul!(similar(U, Base.promote_eltype(U, ΔA)), ΔA, U)
    ∂Λ = mul!(similar(tmp), U', tmp)
    ∂fΛ = _muldiffquotmat!!(∂Λ, f, λ, fλ, df_dλ, ∂Λ)
    # reuse intermediate if possible
    if eltype(tmp) <: Real && eltype(∂fΛ) <: Complex
        tmp2 = ∂fΛ * U'
    else
        tmp2 = mul!(tmp, ∂fΛ, U')
    end
    ∂Y = mul!(∂fΛ, U, tmp2)
    return ∂Y
end

# difference quotient, i.e. Pᵢⱼ = (f(λⱼ) - f(λᵢ)) / (λⱼ - λᵢ), with f'(λᵢ) when λᵢ=λⱼ
function _diffquot(f, λi, λj, fλi, fλj, ∂fλi, ∂fλj)
    T = Base.promote_typeof(λi, λj, fλi, fλj, ∂fλi, ∂fλj)
    Δλ = λj - λi
    iszero(Δλ) && return T(∂fλi)
    # handle round-off error using Maclaurin series of (f(λᵢ + Δλ) - f(λᵢ)) / Δλ wrt Δλ
    # and approximating f''(λᵢ) with forward difference (f'(λᵢ + Δλ) - f'(λᵢ)) / Δλ
    # so (f(λᵢ + Δλ) - f(λᵢ)) / Δλ = (f'(λᵢ + Δλ) + f'(λᵢ)) / 2 + O(Δλ^2)
    # total error on the order of f(λᵢ) * eps()^(2/3)
    abs(Δλ) < cbrt(eps(real(T))) && return T((∂fλj + ∂fλi) / 2)
    Δfλ = fλj - fλi
    return T(Δfλ / Δλ)
end

# broadcast multiply Δ by the matrix of difference quotients P, storing the result in PΔ.
# If β is is nonzero, then @. PΔ = β*PΔ + P*Δ
# if type of PΔ is incompatible with result, new matrix is allocated
function _muldiffquotmat!!(PΔ, f, λ, fλ, ∂fλ, Δ, β = false)
    if eltype(PΔ) <: Real && eltype(fλ) <: Complex
        PΔ2 = similar(PΔ, complex(eltype(PΔ)))
        return _muldiffquotmat!!(PΔ2, f, λ, fλ, ∂fλ, Δ, β)
    else
        PΔ .= β .* PΔ .+ _diffquot.(f, λ, λ', fλ, transpose(fλ), ∂fλ, transpose(∂fλ)) .* Δ
        return PΔ
    end
end

_isindomain(f, x) = true
_isindomain(::Union{typeof(acos),typeof(asin)}, x::Real) = -1 ≤ x ≤ 1
_isindomain(::typeof(acosh), x::Real) = x ≥ 1
_isindomain(::Union{typeof(log),typeof(sqrt)}, x::Real) = x ≥ 0

#####
##### utilities
#####

# Get type (Symmetric or Hermitian) from type or matrix
_symhermtype(::Type{<:Symmetric}) = Symmetric
_symhermtype(::Type{<:Hermitian}) = Hermitian
_symhermtype(A) = _symhermtype(typeof(A))
_symhermtype(a::AbstractThunk) = _symhermtype(unthunk(a))

function _realifydiag!(A)
    for i in diagind(A)
        @inbounds A[i] = real(A[i])
    end
    return A
end

function _symhermlike!(A, S::Union{Symmetric,Hermitian})
    A isa Hermitian{<:Complex} && _realifydiag!(A)
    return typeof(S)(A, S.uplo)
end

# in-place hermitrize matrix
function _hermitrizelike!(A_, S::LinearAlgebra.RealHermSymComplexHerm)
    A = eltype(S) <: Real ? real(A_) : A_
    n = size(A, 1)
    for i in 1:n
        for j in (i + 1):n
            A[i, j] = (A[i, j] + conj(A[j, i])) / 2
            A[j, i] = conj(A[i, j])
        end
        A[i, i] = real(A[i, i])
    end
    return _symhermtype(S)(A, Symbol(S.uplo))
end


text:
# Synchronization Functions

export barrier

const cl_mem_fence_flags = Cuint
const CLK_LOCAL_MEM_FENCE = cl_mem_fence_flags(1)
const CLK_GLOBAL_MEM_FENCE = cl_mem_fence_flags(2)

barrier(flags=0) = @builtin_ccall("barrier", Cvoid, (Cuint,), flags)


text:
# Deps
using Plots
using Random
using StatsBase
using StatsPlots
using Distributions
using KernelDensity
using LinearAlgebra
using DelimitedFiles
default(titlefont=(14,"arial"),guidefont=(16,"arial"),tickfont=(12,"arial"))

######################################################################################
# Functions
######################################################################################
function get_percs(d::UnivariateDistribution,centers::Vector{Float64})

    # Get intervals
    win_size = centers[2]-centers[1]
    intv_vec = [(0.5+win_size*i):(0.5+win_size*(i+1)) for i=0:(length(centers)-1)]

    # Loop over centers
    percs = zeros(length(intv_vec))
    for i=1:length(percs)
        if i==1
            percs[1] = cdf(d,maximum(intv_vec[1]))
        else
            percs[i] = cdf(d,maximum(intv_vec[i]))-cdf(d,maximum(intv_vec[i-1]))
        end
    end

    # Return percentages
    return percs

end
function plot_dist(lens,reso)

    # Histograms
    len_range = 0.0:reso:50000.0
    sts = collect(len_range)
    centers = sts.+step(len_range)/2.0
    histo = fit(Histogram,lens,len_range)
    histo = normalize(histo, mode=:probability)

    # Plot histogram
    len_plot = plot(sts,histo.weights.*100,seriestype=:bar,bins=100,xlim=(0,maximum(len_range)),
        alpha=0.5,xlabel="Read length (bp)",ylabel="Percentage(%)",label="Observed",size=(700,400));

    # Fit distributions
    # gam_fit = fit_mle(Gamma,lens)
    # gam_logpdf = sum(logpdf.(gam_fit,lens))
    # println("Gamma: $(gam_fit.α), $(gam_fit.θ). AIC: $(4-2*gam_logpdf)")
    exp_fit = fit_mle(Exponential,lens)
    exp_logpdf = sum(logpdf.(exp_fit,lens))
    println("Exp: $(exp_fit.θ). AIC: $(2-2*exp_logpdf)")
    # nb_fit = fit_mle(NegativeBinomial,lens)
    # nb_logpdf = sum(logpdf.(nb_fit,lens))
    # println("NB: $(nb_fit.r),$(nb_fit.p). AIC: $(4-2*nb_logpdf)")
    # gam_perc = get_percs(gam_fit,centers)
    exp_perc = get_percs(exp_fit,centers)
    # nb_perc = get_percs(nb_fit,centers)

    # Plot density
    plot!(len_plot,centers,exp_perc.*100,seriestype=:scatter,lw=2,label="Exp");
    # plot!(len_plot,centers,gam_perc.*100,seriestype=:scatter,lw=2,label="Gamma");
    # plot!(len_plot,centers,nb_perc.*100,seriestype=:scatter,lw=2,label="NB");

    # Return plot
    return len_plot

end
######################################################################################
# Raw lengths
######################################################################################

# IO
data_dir = "/Users/jordiabante/OneDrive - Johns Hopkins/CpelNano/Data/Read-Length/"

# Read in histogram data
gm_lens = readdlm("$(data_dir)/gm12878_rel_6_read_lens.txt")[:,1]

# Histograms
p1 = plot_dist(gm_lens,500)
savefig(p1,"$(data_dir)/gm12878_read_lens_plot.png")
savefig(p1,"$(data_dir)/gm12878_read_lens_plot.pdf")

######################################################################################
# Simulated lengths
######################################################################################

# IO
data_dir = "/Users/jordiabante/OneDrive - Johns Hopkins/CpelNano/Data/Read-Length/"

# Read in histogram data
read_len = readdlm("$(data_dir)/simulated_read_lengths.txt")[:,1]

# Histograms
p1 = plot_dist(read_len,500)
savefig(p1,"$(data_dir)/simulated_read_lengths_plot.png")
savefig(p1,"$(data_dir)/simulated_read_lengths_plot.pdf")


text:
"""
    lms(setting; iters = nothing, crit = 2.5)

Perform Least Median of Squares regression estimator with random sampling.

# Arguments
- `setting::RegressionSetting`: A regression setting object.
- `iters::Int`: Number of random samples.
- `crit::Float64`: Critical value for standardized residuals.

# Description
LMS (Least Median of Squares) estimator is highly robust with 50% breakdown property. The algorithm
searches for regression coefficients which minimize (h)th ordered squared residual where h is Int(floor((n + 1.0) / 2.0))

# Output
- `["stdres"]`: Array of standardized residuals
- `["S"]`: Standard error of regression
- `["outliers"]`: Array of indices of outliers
- `["objective"]`: LMS objective value
- `["coef"]`: Estimated regression coefficients
- `["crit"]`: Threshold value.

# Examples
```julia-repl
julia> reg = createRegressionSetting(@formula(calls ~ year), phones);

julia> lms(reg)
Dict{Any,Any} with 6 entries:
  "stdres"    => [2.28328, 1.55551, 0.573308, 0.608843, 0.220321, -0.168202, -0.471913, -0.860435, -0.31603, -0.110871  …  85.7265, 88.9849, 103.269, 116.705, 135.229, 159.69,…
  "S"         => 1.17908
  "outliers"  => [14, 15, 16, 17, 18, 19, 20, 21]
  "objective" => 0.515348
  "coef"      => [-56.1972, 1.1581]
  "crit"      => 2.5
```

# References
Rousseeuw, Peter J. "Least median of squares regression." Journal of the American
statistical association 79.388 (1984): 871-880.
"""
function lms(setting::RegressionSetting; iters=nothing, crit=2.5)
    X, y = @extractRegressionSetting setting
    return lms(X, y, iters=iters, crit=crit)
end

function lms(X::Array{Float64,2}, y::Array{Float64,1}; iters=nothing, crit=2.5)
    n, p = size(X)
    h = Int(floor((n + 1.0) / 2.0))
    if iters === nothing
        iters = minimum([500 * p, 3000])
    end
    bestobjective = Inf
    bestparamaters = []
    bestres = []
    indices = collect(1:n)
    kindices = collect(p:n)
    for iter in 1:iters
        try
            k = rand(kindices, 1)[1]
            sampledindices = sample(indices, k, replace=false)
            olsreg = ols(X[sampledindices,:], y[sampledindices])
            betas = coef(olsreg)
            origres = y .- X * betas
            res = sort(origres.^2.0)
            m2 = res[h]
            if m2 < bestobjective
                bestparamaters = betas
                bestobjective = m2
                bestres = origres
            end
        catch e
            @warn e
        end
    end
    s = 1.4826 * sqrt((1.0 + (5.0 / (n - p))) * bestobjective)
    standardizedres = bestres / s
    d = Dict()
    d["coef"] = bestparamaters
    d["objective"] = bestobjective
    d["S"] = s
    d["stdres"] = standardizedres
    d["crit"] = crit
    d["outliers"] = filter(i -> abs(standardizedres[i]) > crit, 1:n)
    return d
end

max_stars_repo_name: jbytecode/LinRegOutliers

max_issues_repo_name: jbytecode/LinRegOutliers

max_forks_repo_name: jbytecode/LinRegOutliers


text:
using Test
using LogicCircuits
using ProbabilisticCircuits
using DataFrames
using CUDA

@testset "Random Region Graph tests" begin

    num_vars = 10
    depth = 2
    replicas = 2

    region_graph = random_region_graph([Var(i) for i=1:num_vars]; depth, replicas);

    @test isnothing(region_graph.parent)
    @test typeof(region_graph) == ProbabilisticCircuits.RegionGraphInnerNode
    @test length(region_graph.partitions) == 2 # 2 replices means 2 partitions in root
    @test variables(region_graph) == BitSet([Var(i) for i=1:num_vars])

    @test length(region_graph.partitions[1]) == 2 # Parition into two subsets
    @test length(variables(region_graph.partitions[1][1])) == 5
    @test length(variables(region_graph.partitions[1][2])) == 5

    # Each parition should include the same set of variables
    prev_scope = nothing
    for cur_partition in region_graph.partitions
        cur_scope = variables(cur_partition[1])
        for i = 2:length(cur_partition)
            cur_scope = cur_scope ∪ variables(cur_partition[i])
        end

        if !isnothing(prev_scope)
            @test prev_scope == cur_scope
        else
            prev_scope = cur_scope
        end
    end

end

@testset "RAT SPN tests" begin

    data = DataFrame(BitArray([1 0 1 0 1 0 1 0 1 0;
        1 1 1 1 1 1 1 1 1 1;
        0 0 0 0 0 0 0 0 0 0;
        0 1 1 0 1 0 0 1 0 1]), :auto)

    pseudocount = 0.1
    num_vars = num_features(data)

    num_vars = 10
    depth = 2
    replicas = 2

    num_nodes_region = 3
    num_nodes_leaf   = 4

    region_graph = random_region_graph([Var(i) for i=1:num_vars]; depth, replicas);
    circuit = region_graph_2_pc(region_graph; num_nodes_root = 1, num_nodes_region, num_nodes_leaf)[1];

    @test typeof(circuit) <: ProbCircuit

    estimate_parameters_em(circuit, data; pseudocount, use_gpu=false, update_per_batch = false)
    estimate_parameters_em(circuit, data; pseudocount, use_gpu=false, update_per_batch = true)

    if CUDA.functional()
        estimate_parameters_em(circuit, data; pseudocount, use_gpu=true, update_per_batch = false)
        estimate_parameters_em(circuit, data; pseudocount, use_gpu=true, update_per_batch = true)
    end

end


text:
# -----------------------------------------------------------------------

type Layout <: Widget
  widget::PyObject  # QFrame with the given layout
  children::Vector{Widget}
end

function Layout(isvert::Bool, children::Widget...)

  layout = Layout(QT.QFrame(), Widget[])
  qlayout = isvert ? QT.QVBoxLayout() : QT.QHBoxLayout()
  layout.widget[:setLayout](qlayout)
  for child in children
    push!(layout.children, child)
    qlayout[:addWidget](child.widget)
  end
  layout
end

vbox(args...) = Layout(true, args...)
hbox(args...) = Layout(false, args...)

type Splitter <: Widget
  widget::PyObject  # QSplitter
  children::Vector{Widget}
end

function Splitter(isvert::Bool, children::Widget...)
  splitter = Splitter(QT.QSplitter(isvert ? 2 : 1), Widget[])
  for child in children
    push!(splitter.children, child)
    splitter.widget[:addWidget](child.widget)
  end
  splitter
end

vsplitter(args...) = Splitter(true, args...)
hsplitter(args...) = Splitter(false, args...)


text:
export
  # CloudGraphs helper functions
  insertnodefromcv!,
  checkandinsertedges!,
  getbinarraymongo,
  gettopoint,
  getdotwothree,
  bin2arr,
  fetchsubgraph!,
  getVertNeoIDs!,
  insertrobotdatafirstpose!,
  tryunpackalltypes!,
  fetchrobotdatafirstpose,
  getExVertexNeoIDs,
  db2jld,

  # solver service SLAMinDB
  getcredentials,
  startSlamInDb,
  runSlamInDbOnSession,
  slamindb,
  convertdb,
  resetconvertdb,
  getmaxfactorid,
  # would be CloudGraphs calls
  hasBigDataElement,
  getBigDataElement,
  removeNeo4jID,

  # webserver
  SolverStatus,
  CaesarConfig,
  IterationStatistics,
  VisualizationConfig

    # # multisession utils
    # multisessionquery,
    # parsemultisessionqueryresult!,
    # getLandmOtherSessNeoIDs,
    # getAllLandmarkNeoIDs,
    # getLocalSubGraphMultisession,
    # findExistingMSConstraints,
    # getprpt2kde,
    # rmInstMultisessionPriors!,
    # removeMultisessions!,


text:
# https://github.com/morhetz/gruvbox
const _gruvbox_dark_palette =  (:bright_green, :bright_yellow, :bright_blue, :bright_aqua, :bright_purple, :bright_red, :bright_orange)

const _gruvbox_light_palette = (:faded_green, :faded_yellow, :faded_blue, :faded_aqua, :faded_purple, :faded_red, :faded_orange)

const _gruvbox_colors = Dict(
	:dark0_hard     => RGB255(29,32,33),
	:dark0          => RGB255(40,40,40),
	:dark0_soft     => RGB255(50,48,47),
	:dark1          => RGB255(60,56,54),
	:dark2          => RGB255(80,73,69),
	:dark3          => RGB255(102,92,84),
	:dark4          => RGB255(124,111,100),
	:dark4_256      => RGB255(124,111,100),
	:gray_245       => RGB255(146,131,116),
	:gray_244       => RGB255(146,131,116),
	:light0_hard    => RGB255(249,245,215),
	:light0         => RGB255(253,244,193),
	:light0_soft    => RGB255(242,229,188),
	:light1         => RGB255(235,219,178),
	:light2         => RGB255(213,196,161),
	:light3         => RGB255(189,174,147),
	:light4         => RGB255(168,153,132),
	:light4_256     => RGB255(168,153,132),
	:bright_red     => RGB255(251,73,52),
	:bright_green   => RGB255(184,187,38),
	:bright_yellow  => RGB255(250,189,47),
	:bright_blue    => RGB255(131,165,152),
	:bright_purple  => RGB255(211,134,155),
	:bright_aqua    => RGB255(142,192,124),
	:bright_orange  => RGB255(254,128,25),
	:neutral_red    => RGB255(204,36,29),
	:neutral_green  => RGB255(152,151,26),
	:neutral_yellow => RGB255(215,153,33),
	:neutral_blue   => RGB255(69,133,136),
	:neutral_purple => RGB255(177,98,134),
	:neutral_aqua   => RGB255(104,157,106),
	:neutral_orange => RGB255(214,93,14),
	:faded_red      => RGB255(157,0,6),
	:faded_green    => RGB255(121,116,14),
	:faded_yellow   => RGB255(181,118,20),
	:faded_blue     => RGB255(7,102,120),
	:faded_purple   => RGB255(143,63,113),
	:faded_aqua     => RGB255(66,123,88),
	:faded_orange   => RGB255(175,58,3),
)

_themes[:gruvbox_dark] = PlotTheme(
	bg = _gruvbox_colors[:dark2],
	bginside = _gruvbox_colors[:dark0],
	fg = _gruvbox_colors[:light3],
	fgtext = _gruvbox_colors[:light3],
	fgguide = _gruvbox_colors[:dark1],
	fglegend = _gruvbox_colors[:light3],
	palette = expand_palette(_gruvbox_colors[:dark3], [_gruvbox_colors[c] for c in _gruvbox_dark_palette]),
    gradient = cgrad(:YlOrRd_r).colors
)

_themes[:gruvbox_light] = PlotTheme(
	bg = _gruvbox_colors[:light1],
	bginside = _gruvbox_colors[:light0],
	fg = _gruvbox_colors[:dark1],
	fgtext = _gruvbox_colors[:dark1],
	fgguide = _gruvbox_colors[:dark1],
	fglegend = _gruvbox_colors[:dark1],
	palette = expand_palette(_gruvbox_colors[:light3], [_gruvbox_colors[c] for c in _gruvbox_light_palette]),
    gradient = cgrad(:YlOrRd_r).colors
)


text:
function readcddb(fnm::AbstractString)
    nms = ASCIIString[]
    vv = ByteString[]
    open(fnm,"r") do io
        for ln in eachline(io)
            ismatch(r"^\s*#",ln) && next
            vals = map(strip,split(chomp(ln),"="))
            if length(vals) == 2 && length(vals[2]) > 0
                push!(nms,vals[1])
                push!(vv,vals[2])
            end
        end
    end
    nms,vv
end

function readinf(nm::ByteString)
    res = Dict{ASCIIString,ByteString}()
    if isfile(nm)
        open(nm,"r") do io
            for ln in eachline(io)
                vals = map(strip,split(split(chomp(ln),"#")[1],"="))
                if length(vals) == 2 && length(vals[2]) > 0
                    res[vals[1]] = strip(vals[2],'\'')
                end
            end
        end
    end
    res
end

function process_dir(dnm::ByteString)
    cd(dnm) do
        for nm in readdir()
            bb,ee = splitext(nm)
            if ee == ".wav"
                meta = readinf(string(bb,".inf"))
                hdr,


text:
module Storages

using ...PCG.Types: Turn

export Storage, StorageIndex, StorageNode, get_node, set_node!, apply_changes! #, record_state!, recorders

abstract type Storage end
abstract type StorageIndex end

# TODO: mark parent classes for template parameters
struct StorageNode{P}
    current::P
    new::P
    changed_at::Turn

    StorageNode(current::P) where P = new{P}(current, current, 0)
    StorageNode(current::P, turn::Turn) where P = new{P}(current, current, turn)
    StorageNode(current::P, new::P, turn::Turn) where P = new{P}(current, new, turn)
end

# TODO: does that enough to get plain data layout in array?
const StorageNodes{P} = Vector{StorageNode{P}}

# TODO: specify types?

function get_node end

function set_node! end

function apply_changes! end

module LinearStorages

using ...PCG.Storages
using ...PCG.Storages: StorageIndex, Storage, StorageNodes, StorageNode, Turn

export LinearStorageIndex, LinearStorage

struct LinearStorageIndex <: StorageIndex
    i::Int64
end

const LinearStorageIndexes = Vector{LinearStorageIndex}

struct LinearStorage{P} <: Storage
    nodes::StorageNodes{P}
    new_nodes::LinearStorageIndexes

    function LinearStorage(base_properties::P, size::Int64) where P
        base_node = StorageNode(base_properties)

        nodes = fill(base_node, size)

        new_nodes = LinearStorageIndexes()
        sizehint!(new_nodes, size)

        return new{P}(nodes, new_nodes)
    end
end

function Storages.get_node(storage::LinearStorage, i::LinearStorageIndex)
    return storage.nodes[i.i]
end

function Storages.set_node!(storage::LinearStorage, i::LinearStorageIndex, node::StorageNode)
    storage.nodes[i.i] = node
    push!(storage.new_nodes, i)
    return
end

function Storages.apply_changes!(storage::LinearStorage, turn::Turn)
    for i in storage.new_nodes
        if storage.nodes[i.i].changed_at == turn
            storage.nodes[i.i] = StorageNode(storage.nodes[i.i].new, turn)
        end
    end

    resize!(storage.new_nodes, 0)
end

end

end

max_stars_repo_name: Tiendil/pcg-julia

max_issues_repo_name: Tiendil/pcg-julia

max_forks_repo_name: Tiendil/pcg-julia


text:
module DropboxCLI

using ArgParse
using Base.Iterators

using DropboxSDK

const arg_settings = ArgParseSettings()

add_arg_table(
    arg_settings,
    "account", Dict(:help => "show account information",
                    :action => :command),
    "cmp", Dict(:help => "compare (content hash of) files or directories",
                :action => :command),
    "du", Dict(:help => "show disk usage",
               :action => :command),
    "get", Dict(:help => "get files or directories",
                :action => :command),
    "ls", Dict(:help => "list directory content",
               :action => :command),
    "mkdir", Dict(:help => "create new directory",
                  :action => :command),
    "put", Dict(:help => "put files or directories",
                :action => :command),
    "rm", Dict(:help => "delete file or directory",
               :action => :command),
    "version", Dict(:help => "delete file or directory",
                    :action => :command),
    ["--verbose", "-v"], Dict(:help => "output progress information",
                              :nargs => 0),
    ["--debug", "-g"], Dict(:help => "output detailed progress information",
                            :nargs => 0),
)

add_arg_table(
    arg_settings["cmp"],
    "filename", Dict(:help => "name of local file or directory to compare",
                     :nargs => '*'),
    "destination", Dict(:help => "destination Dropbox file or directory",
                        :nargs => 'A'),
)

add_arg_table(
    arg_settings["get"],
    "filename", Dict(:help => "name of Dropbox file or directory to get",
                     :nargs => '*'),
    "destination", Dict(:help => "local destination file or directory",
                        :nargs => 'A'),
)

add_arg_table(
    arg_settings["ls"],
    ["--long", "-l"], Dict(:help => "use a long (detailed) format",
                           :nargs => 0),
    ["--recursive", "-R"], Dict(:help => "recursively list subdirectories",
                                :nargs => 0),
    "filename", Dict(:help => "file or directory name",
                     :nargs => '*'),
)

add_arg_table(
    arg_settings["mkdir"],
    "directoryname", Dict(:help => "name of directory to create",
                       :nargs => '*'),
)

add_arg_table(
    arg_settings["put"],
    "filename", Dict(:help => "name of local file or directory to put",
                     :nargs => '*'),
    "destination", Dict(:help => "destination Dropbox file or directory",
                        :nargs => 'A'),
)

add_arg_table(
    arg_settings["rm"],
    "filename", Dict(:help => "name of file or directory to remove",
                     :nargs => '*'),
    ["--recursive", "-r"], Dict(:help => "recursively delete subdirectories",
                                :nargs => 0),
)

# 0: no output
# 1: some output
# 2: detailed output
const verbose = Ref{Int}(0)
function parse_verbose(opts)
    if opts["verbose"]
        verbose[] = 1
    end
    if opts["debug"]
        verbose[] = 2
    end
end

function quote_string(str::AbstractString)::AbstractString
    repr(str)[2:end-1]
end

function cmd_account(args)
    auth = get_authorization()
    account = users_get_current_account(auth)
    first = account.name.given_name
    last = account.name.surname
    display = account.name.display_name
    println("Account: Name: $first $last ($display)")
    return 0
end

function cmd_cmp(args)
    filenames = args["filename"]
    if length(filenames) == 0
        println("Error: No file names given")
        return 1
    elseif length(filenames) == 1
        println("Error: Destination missing")
        return 1
    end
    # The last file name is the destination
    @assert !haskey(args, "target")
    destination = filenames[end]
    sources = filenames[1:end-1]

    exit_code = 0

    auth = get_authorization()

    # Add leading and remove trailing slashes
    if !startswith(destination, "/")
        destination = "/$destination"
    end
    while endswith(destination, "/")
        destination = destination[1:end-1]
    end

    if destination == ""
        # The root directory does not support files_get_metadata
        ispath_destination = true
        isdir_destination = true
    else
        try
            metadata = files_get_metadata(auth, destination)
            ispath_destination = true
            isdir_destination = metadata isa FolderMetadata
        catch ex
            if ex isa DropboxError
                if ex.dict["error"][".tag"] == "path" &&
                    ex.dict["error"]["path"][".tag"] == "not_found"
                    ispath_destination = false
                    isdir_destination = false
                else
                    println("$(quote_string(filename)):",
                            " $(ex.dict["error_summary"])")
                    return 1
                end
            else
                rethrow(ex)
            end
        end
    end

    if !isdir_destination && length(sources) != 1
        # Multiple sources: Destination is not a directory
        println("Destination directory \"$destination\" does not exist")
        return 1
    end

    for source in sources

        # Distinguish between files and directories
        if !ispath(source)
            println("$source: File not found")
            continue
        end
        isdir_source = isdir(source)

        if isdir_destination
            # If the destination is a directory, the files will be
            # downloaded into that directory.
            filename = joinpath(destination, basename(source))
        elseif length(sources) == 1 && !isdir_source
            # If there is only a single source
            if ispath_destination
                # If the destination exists and is not a directory,
                # then it is overwritten.
                @assert !isdirpath(destination)
                filename = destination
            else
                # If the destination does not exist, then a file with
                # that name will be created.
                filename = destination
                pathname = dirname(filename)
            end
        else
            # Multiple sources: Destination is not a directory
            @assert false
        end

        # Add leading and remove trailing slashes
        if !startswith(filename, "/")
            filename = "/$filename"
        end
        while endswith(filename, "/")
            filename = filename[1:end-1]
        end

        # TODO: Handle sources that are a directory
        @assert !isdir_source

        # Compare content hash
        metadata = try
            files_get_metadata(auth, filename)
        catch ex
            if ex isa DropboxError
                println("$(quote_string(filename)):",
                        " $(ex.dict["error_summary"])")
                exit_code = 1
                continue
            else
                rethrow(ex)
            end
        end
        if !(metadata isa FileMetadata)
            println("$(quote_string(source)): Not a file")
            exit_code = 2
        else
            size = filesize(source)
            if metadata.size != size
                println("$(quote_string(source)): File size differs")
                exit_code = 2
            else
                content = read(source)
                content_hash = calc_content_hash(content)
                if metadata.content_hash != content_hash
                    println("$(quote_string(source)): Content hash differs")
                    exit_code = 2
                else
                    # no difference
                end
            end
        end

    end

    return exit_code
end

const metric_prefixes = ["", "k", "M", "G", "T", "P", "E", "Z", "Y"]
function find_prefix(x::Real)::Tuple{Int64, String}
    for (exp3, prefix) in enumerate(metric_prefixes)
        scale = 1000.0^(exp3-1)
        if abs(x) < 1000.0*scale
            return scale, prefix
        end
    end
    return 1, ""
end

function cmd_du(args)
    auth = get_authorization()
    usage = users_get_space_usage(auth)
    used = usage.used
    allocated = usage.allocation.allocated

    max_bytes = max(used, allocated)
    bytes_digits = length(string(max_bytes))

    scale, prefix = find_prefix(allocated)
    println("allocated:",
            " $(lpad(allocated, bytes_digits)) bytes",
            " ($(round(allocated / scale, sigdigits=3)) $(prefix)Byte)")

    scale, prefix = find_prefix(used)
    pct = round(100 * used / allocated, digits=1)
    println("used:     ",
            " $(lpad(used, bytes_digits)) bytes",
            " ($(round(used / scale, sigdigits=3)) $(prefix)Byte, $pct%)")

    return 0
end

function cmd_get(args)
    filenames = args["filename"]
    if length(filenames) == 0
        println("Error: No file names given")
        return 1
    elseif length(filenames) == 1
        println("Error: Destination missing")
        return 1
    end
    # The last file name is the destination
    @assert !haskey(args, "target")
    destination = filenames[end]
    sources = filenames[1:end-1]

    exit_code = 0

    auth = get_authorization()

    for source in sources

        # Add leading and remove trailing slashes
        if !startswith(source, "/")
            source = "/$source"
        end
        @assert !endswith(source, "/")

        # Distinguish between files and directories
        if source == ""
            # The root directory does not support files_get_metadata
            isdir_source = true
        else
            metadata = try
                files_get_metadata(auth, source)
            catch ex
                if ex isa DropboxError
                    println("$(quote_string(source)):",
                            " $(ex.dict["error_summary"])")
                    exit_code = 1
                    continue
                else
                    rethrow(ex)
                end
            end
            isdir_source = metadata isa FolderMetadata
        end

        if isdir(destination)
            # If the destination is a directory, the files will be
            # downloaded into that directory.
            filename = joinpath(destination, basename(source))
        elseif length(sources) == 1 && !isdir_source
            # If there is only a single source
            if ispath(destination)
                # If the destination exists and is not a directory,
                # then it is overwritten.
                @assert !isdirpath(destination)
                filename = destination
            else
                # If the destination does not exist, then a file with
                # that name will be created.
                filename = destination
                pathname = dirname(filename)
                @assert isdir(pathname)
            end
        else
            # Multiple sources: Destination is not a directory
            @assert false
        end

        # TODO: Handle sources that are a directory
        @assert !isdir_source

        # Compare content hash before downloading
        need_download = true
        if metadata.size == 0
            open(filename, "w") do io
                truncate(io, 0)
            end
            need_download = false
        elseif isfile(filename)
            size = filesize(filename)
            if size == metadata.size
                # Don't download if content hash matches
                content = read(filename)
                content_hash = calc_content_hash(content)
                if content_hash == metadata.content_hash
                    if verbose[] >= 2
                        println("Info: $(quote_string(filename)):",
                                " content hash matches; skipping download")
                    end
                    need_download = false
                end
            elseif size < metadata.size
                # TODO: Download only missing fraction
            else
                # Truncate if necessary
                content = read(filename, metadata.size)
                content_hash = calc_content_hash(content)
                if content_hash == metadata.content_hash
                    if verbose[] >= 2
                        println("Info: $(quote_string(filename)):",
                                " content hash matches;",
                                " truncating local file and skipping download")
                    end
                    open(filename, "w") do io
                        truncate(io, metadata.size)
                    end
                    need_download = false
                end
            end
        end

        # TODO: touch file if download is skipped?
        if need_download
            # TODO: download all files simultaneously
            metadata, content = files_download(auth, source)
            open(filename, "w") do io
                write(io, content)
            end
        end

    end

    return exit_code
end

const mode_strings = Dict{Type, String}(
    FileMetadata => "-",
    FolderMetadata => "d",
    DeletedMetadata => "?",
)

metadata_size(metadata) = Int64(0)
metadata_size(metadata::FileMetadata) = metadata.size

# e.g. "2019-03-15T23:05:18Z"
metadata_modified(metadata) = "                    " # 20 spaces
metadata_modified(metadata::FileMetadata) = metadata.server_modified

function metadata_path(metadata, prefix)
    path = metadata.path_display
    if startswith(path, prefix)
        path = path[length(prefix)+1 : end]
    end
    if startswith(path, "/")
        path = path[2:end]
    end
    if path == ""
        path = "."
    end
    path
end

function cmd_ls(args)
    long = args["long"]
    recursive = args["recursive"]
    filenames = args["filename"]
    if isempty(filenames)
        filenames = [""]        # show root
    end

    exit_code = 0

    auth = get_authorization()

    # TODO: Sort filenames by type. First show all Files, then all
    # Directories prefixed with their names. Also sort everything
    # alphabetically.
    for filename in sort(filenames)

        # Add leading and remove trailing slashes
        if !startswith(filename, "/")
            filename = "/$filename"
        end
        while endswith(filename, "/")
            filename = filename[1:end-1]
        end

        # Distinguish between files and directories
        if filename == ""
            # The root directory does not support files_get_metadata
            isdirectory = true
        else
            metadata = try
                files_get_metadata(auth, filename)
            catch ex
                if ex isa DropboxError
                    println("$(quote_string(filename)):",
                            " $(ex.dict["error_summary"])")
                    exit_code = 1
                    continue
                else
                    rethrow(ex)
                end
            end
            isdirectory = metadata isa FolderMetadata
        end

        if isdirectory
            metadatas = files_list_folder(auth, filename, recursive=recursive)
            prefix_to_hide = filename
        else
            metadatas = [metadata]
            prefix_to_hide = ""
        end

        # Output directory name if there are multiple directories
        if length(filenames) > 1
            println()
            println("$(quote_string(filename)):")
        end

        # Sort
        sort!(metadatas, by = m -> metadata_path(m, prefix_to_hide))

        # Output
        if long
            max_size =
                maximum(metadata_size(metadata) for metadata in metadatas)
            size_digits = length(string(max_size))
            for metadata in metadatas
                isdir = metadata isa FolderMetadata
                mode = mode_strings[typeof(metadata)]
                size = isdir ? "" : metadata_size(metadata)
                modified = metadata_modified(metadata)
                path = metadata_path(metadata, prefix_to_hide)
                println("$mode $(lpad(size, size_digits)) $modified",
                        " $(quote_string(path))")
            end
        else
            for metadata in metadatas
                path = metadata_path(metadata, prefix_to_hide)
                println("$(quote_string(path))")
            end
        end

    end

    return exit_code
end

function cmd_mkdir(args)
    directorynames = args["directoryname"]

    exit_code = 0

    auth = get_authorization()
    for directoryname in directorynames

        # Add leading and remove trailing slashes
        if !startswith(directoryname, "/")
            directoryname = "/" * directoryname
        end
        while endswith(directoryname, "/")
            directoryname = directoryname[1:end-1]
        end

        try
            files_create_folder(auth, directoryname)
        catch ex
            if ex isa DropboxError
                println("$(quote_string(directoryname)):",
                        " $(ex.dict["error_summary"])")
                exit_code = 1
                continue
            else
                rethrow(ex)
            end
        end

    end

    return exit_code
end

function cmd_put(args)
    filenames = args["filename"]
    if length(filenames) == 0
        println("Error: No file names given")
        return 1
    elseif length(filenames) == 1
        println("Error: Destination missing")
        return 1
    end
    # The last file name is the destination
    @assert !haskey(args, "target")
    destination = filenames[end]
    sources = filenames[1:end-1]

    exit_code = 0

    auth = get_authorization()

    # Add leading and remove trailing slashes
    if !startswith(destination, "/")
        destination = "/$destination"
    end
    while endswith(destination, "/")
        destination = destination[1:end-1]
    end

    if destination == ""
        # The root directory does not support files_get_metadata
        ispath_destination = true
        isdir_destination = true
    else
        try
            metadata = files_get_metadata(auth, destination)
            ispath_destination = true
            isdir_destination = metadata isa FolderMetadata
        catch ex
            if ex isa DropboxError
                if ex.dict["error"][".tag"] == "path" &&
                    ex.dict["error"]["path"][".tag"] == "not_found"
                    ispath_destination = false
                    isdir_destination = false
                else
                    println("$(quote_string(destination)):",
                            " $(ex.dict["error_summary"])")
                    return 1
                end
            else
                rethrow(ex)
            end
        end
    end

    # If the destination is not a directory, i.e. if it is a regular
    # file, or if it does not exist, then there can be only a single
    # source.
    if !isdir_destination && length(sources) > 1
        println("Destination \"$destination\" exists, but is not a directory")
        return 1
    end

    metadata_dict = Dict{String, Metadata}()

    upload_channel = Channel{Tuple{String, String}}(1000)
    upload_task =
        start_task(() -> upload_many_files(auth, metadata_dict, upload_channel),
                   (:cmd_put, :upload_many_files, sources, destination))
    for source in sources
        # Distinguish between files and directories
        if !ispath(source)
            println("$source: File not found")
            continue
        end
        isdir_source = isdir(source)

        if isdir_destination
            # If the destination is a directory, the files will be
            # uploaded into that directory.
            filename = joinpath(destination, basename(source))
        else
            # If the destination is not a directory, then it specifies
            # the destination file name (which may or may not exist)
            @assert length(sources) == 1
            filename = destination
        end

        # Add leading and remove trailing slashes
        if !startswith(filename, "/")
            filename = "/$filename"
        end
        while endswith(filename, "/")
            filename = filename[1:end-1]
        end

        if !haskey(metadata_dict, lowercase(filename))
            metadata =
                if filename != ""
                    try
                        files_get_metadata(auth, filename)
                    catch ex
                        if ex isa DropboxError
                            nothing
                        else
                            rethrow(ex)
                        end
                    end
                else
                    nothing
                end
            if metadata !== nothing
                metadata_dict[lowercase(filename)] = metadata
            end
            if filename == "" || metadata isa FolderMetadata
                try
                    metadatas =
                        files_list_folder(auth, filename, recursive=true)
                    for m in metadatas
                        # @assert m.path_display !== nothing
                        # metadata_dict[lowercase(m.path_display)] = m
                        @assert m.path_lower !== nothing
                        metadata_dict[m.path_lower] = m
                    end
                catch ex
                    if ex isa DropboxError
                        nothing
                    else
                        rethrow(ex)
                    end
                end
            end
        end

        if !isdir_source
            # The source is a regular file
            put!(upload_channel, (source, filename))
        else
            # Both source and destination are directories
            function walk_upload_entry(source_prefix::String,
                                       filename_prefix::String)
                for entry in readdir(source)
                    source = joinpath(source_prefix, entry)
                    filename = "$filename_prefix/$entry"
                    if isdir(source)
                        walk_upload_entry(source, filename)
                    else
                        put!(upload_channel, (source, filename))
                    end
                end
            end
            walk_upload_entry(source, filename)
        end
    end
    close(upload_channel)
    wait(upload_task)

    return exit_code
end

function upload_many_files(auth::Authorization,
                           metadata_dict::Dict{String, Metadata},
                           upload_channel::Channel{Tuple{String, String}}
                           )::Nothing
    # Use batches of at most N files, and which upload in at most S
    # seconds
    max_files = 500             # Dropbox limit is 1000
    max_seconds = 300.0         # 5 minutes
    max_tasks = 10

    finalize_task = nothing
    old_tasks = Task[]

    i = 0
    n = Ref(0)
    while isready(upload_channel) || isopen(upload_channel)
        # Where n is still changing
        np = isopen(upload_channel)

        nuploads = 0
        upload_spec_channel = Channel{UploadSpec}(0)
        upload_spec_task =
            start_task(() -> files_upload(auth, upload_spec_channel),
                       (:upload_many_files, :files_upload))

        tasks = Task[]

        start_time = time()
        for (source, destination) in upload_channel
            i += 1
            n[] = i + Base.n_avail(upload_channel)

            j,m = i,n[]
            push!(tasks,
                  start_task(() -> upload_one_file(auth, metadata_dict,
                                                   j, n, np,
                                                   source, destination,
                                                   upload_spec_channel),
                             (:upload_many_files, :upload_one_file,
                              j, m, np, source, destination)))
            while length(tasks) + length(old_tasks) >= max_tasks
                yield()
                function check_task(t)
                    if istaskdone(t)
                        did_upload = fetch(t)
                        nuploads += did_upload
                        return false
                    end
                    return true
                end
                filter!(check_task, tasks)
            end

            upload_cond = nuploads + length(tasks) >= max_files
            time_cond = time() - start_time >= max_seconds
            (upload_cond || time_cond) && break
        end

        function finalize(i::Int, n::Ref{Int}, np::Bool)
            while !isempty(old_tasks)
                yield()
                filter!(!istaskdone, old_tasks)
            end
            if verbose[] >= 1
                println("Info: Flushing upload ($i/$(n[])$(np ? "+" : ""))")
                flush(stdout)
            end
            close(upload_spec_channel)
            wait(upload_spec_task)
        end
        finalize_task !== nothing && wait(finalize_task)
        @assert isempty(old_tasks)
        # copy!(old_tasks, tasks)
        resize!(old_tasks, length(tasks))
        old_tasks[:] = tasks[:]
        empty!(tasks)
        j,m = i,n[]
        finalize_task = start_task(() -> finalize(j, n, np),
                                   (:upload_many_files, :finalize, j, m, np))

    end

    finalize_task !== nothing && wait(finalize_task)
    @assert isempty(old_tasks)

    return nothing
end

function upload_one_file(auth::Authorization,
                         metadata_dict::Dict{String, Metadata},
                         i::Int, n::Ref{Int}, np::Bool,
                         source::String, destination::String,
                         upload_spec_channel::Channel{UploadSpec})::Bool
    # Compare content hash before uploading
    # TODO: remember missing parent directories, and don't try to get
    # a content hash from Dropox in that case
    need_upload = true
    content = nothing

    metadata = get(metadata_dict, lowercase(destination), nothing)

    if metadata isa FileMetadata
        size = filesize(source)
        if metadata.size == size
            # Don't upload if content hash matches
            if verbose[] >= 2
                println("Info: Comparing content hash ($i/$(n[])):",
                        " $(quote_string(source))")
            end
            # content = read(source)
            # content_hash = calc_content_hash(content)
            # Read in chunks
            data_channel = Channel{Vector{UInt8}}(0)
            content_hash_task =
                start_task(() -> calc_content_hash(data_channel),
                           (:upload_one_file, :calc_content_hash,
                            i, n, source, destination))
            open(source, "r") do io
                while !eof(io)
                    chunksize = 4 * 1024 * 1024
                    chunk = read(io, chunksize)
                    put!(data_channel, chunk)
                end
            end
            close(data_channel)
            content_hash = fetch(content_hash_task)
            if metadata.content_hash == content_hash
                if verbose[] >= 2
                    println(
                        "Info: Content hash matches;",
                        " skipping upload ($i/$(n[])$(np ? "+" : "")):",
                        " $(quote_string(source))")
                end
                need_upload = false
            end
        elseif metadata.size < size
            # TODO: Upload only missing fraction
        end
        if need_upload
            # The file exists, but is different. We need to delete it
            # before it can be re-uploaded.
            # TODO: Upload to a temporary name, then rename.
            # TODO: Use WriteMode=overwrite instead of deleting.
            if verbose[] >= 2
                println(
                    "Info: Content hash differs;",
                    " deleting file first ($i/$(n[])$(np ? "+" : "")):",
                    " $(quote_string(source))")
            end
            files_delete(auth, destination)
        end
    end

    # TODO: touch Dropbox file if upload is skipped?
    if need_upload
        if verbose[] >= 1
            println("Info: Uploading ($i/$(n[])$(np ? "+" : "")):",
                    " $(quote_string(source))")
        end
        try
            open(source, "r") do io
                content_channel = Channel{Vector{UInt8}}(0)
                put!(upload_spec_channel,
                     UploadSpec(destination, content_channel))
                # Read in chunks of 150 MByte
                chunksize = 150 * 1024 * 1024
                bytes_total = filesize(source)
                bytes_read = Int64(0)
                start_time = time()
                while !eof(io)
                    # println("Info: Debug: reading chunk ($i/$(n[]))...")
                    chunk = read(io, chunksize)
                    bytes_read += length(chunk)
                    # println("Info: Debug: sending chunk ($i/$(n[]))...")
                    put!(content_channel, chunk)
                    if bytes_total == 0 || bytes_read == 0
                        ratio = 1.0
                        rate = 0
                    else
                        ratio = bytes_read / bytes_total
                        rate = bytes_read / max(1.0, time() - start_time)
                    end
                    bscale, bprefix = find_prefix(bytes_read)
                    rscale, rprefix = find_prefix(rate)
                    if verbose[] >= 2
                        println(
                            "Info: Uploading ($i/$(n[])$(np ? "+" : "")):",
                            " $(quote_string(source))",
                            " ($(round(bytes_read / bscale, sigdigits=3))",
                            " $(bprefix)B,",
                            " $(round(Int, 100 * ratio))%,",
                            " $(round(rate / rscale, sigdigits=3))",
                            " $(rprefix)B/s)")
                    end
                end
                close(content_channel)
            end
        catch ex
            println("Warning: Could not read file $(quote_string(source))")
            need_upload = false
        end
    end

    return need_upload
end

function cmd_rm(args)
    filenames = args["filename"]

    exit_code = 0

    auth = get_authorization()
    for filename in filenames

        # Add leading and remove trailing slashes
        if !startswith(filename, "/")
            filename = "/" * filename
        end
        while endswith(filename, "/")
            filename = filename[1:end-1]
        end

        try
            files_delete(auth, filename)
        catch ex
            if ex isa DropboxError
                println("$(quote_string(filename)):",
                        " $(ex.dict["error_summary"])")
                exit_code = 1
                continue
            else
                rethrow(ex)
            end
        end

    end

    return exit_code
end

const re_name_to_string    = r"^\s*name\s*=\s*\"(.*)\"\s*(?:#|$)"
const re_version_to_string = r"^\s*version\s*=\s*\"(.*)\"\s*(?:#|$)"

function cmd_version(args)
    project_filename = joinpath(dirname(dirname(pathof(DropboxSDK))),
                                "Project.toml")
    name = nothing
    version = nothing
    open(project_filename) do io
        for line in eachline(io)
            if (m = match(re_name_to_string, line)) != nothing
                name = String(m.captures[1])
            elseif (m = match(re_version_to_string, line)) != nothing
                version = VersionNumber(m.captures[1])
            end
        end
    end
    @assert name == "DropboxSDK"

    println("Version $version")

    return 0
end

const cmds = Dict(
    "account" => cmd_account,
    "cmp" => cmd_cmp,
    "du" => cmd_du,
    "get" => cmd_get,
    "ls" => cmd_ls,
    "mkdir" => cmd_mkdir,
    "put" => cmd_put,
    "rm" => cmd_rm,
    "version" => cmd_version,
)

function main(args)
    println("Julia Dropbox client   ",
            "<https://github.com/eschnett/DropboxSDK.jl>")
    println("(modified)")
    println()
    opts = parse_args(args, arg_settings)
    parse_verbose(opts)
    cmd = opts["%COMMAND%"]
    fun = cmds[cmd]
    code = fun(opts[cmd])
    return code
end

end