text:
export jacobian, derivative, @posdef, @odd, @skewsym, square
export makeposdef, makeodd, makeskewsym

function ltril(v::AbstractVector)
    N = length(v)
    M = isqrt(N)
    LowerTriangular(reshape(v,M,M))
end
function ltril(V::AbstractMatrix)
    cat(map(ltril, eachcol(V))..., dims=(1,2))
end
function posdef(v::AbstractVecOrMat,::Any=nothing)
    L = ltril(v)
    L*L'
end

function skewsym(v::AbstractVector,::Any=nothing)
    N = floor(Int, sqrt(length(v)))
    A = reshape(v,N,N)
    S = A - A'
end
function skewsym(V::AbstractMatrix,::Any=nothing)
    cat(map(skewsym, eachcol(V))..., dims=(1,2))
end

const SubMatrix = SubArray{T,2} where {T}

function hadamard(A::AbstractMatrix,B::M) where {M<:Union{AbstractMatrix,SubMatrix}}
    nx = size(A,2)
    m,n = size(B)
    res = zeros(eltype(A), (m*nx, n*nx))
    for i = 1:nx
        res[m*(i-1)+1:m*i,n*(i-1)+1:n*i] = A[:,i] .* B
    end
    return res
end
hadamard(A::AbstractVector, B::M) where {M<:Union{AbstractMatrix,SubMatrix}} = A .* B

delu(x::Real, α=one(x)) = ifelse(x > 0.0, one(x), α*exp(x) )
dtanh(x) = one(x) - tanh(x)*tanh(x)
square(x,::Any=nothing) = x.^2

derivative(::typeof(elu)) = delu
derivative(::typeof(tanh)) = dtanh
derivative(::typeof(identity)) = (x)->one(eltype(x))
derivative(::typeof(square)) = (x,::Any=nothing)->2x

jacobian(::typeof(square), x::AbstractVector, ::Any=nothing) = diagm(derivative(square)(x))
jacobian(l::Dense, x) = hadamard(derivative(l.σ).(x), l.W)
function jacobian(l::FastDense, x, p)
    W, _ = param2Wb(l, p)
    σbar = derivative(l.σ).(x)
    return hadamard(σbar, W)
end
function jacobian(::F, x, ::Any=nothing) where {F<:Function}
    error("Analytical jacobian of $F required, but not implemented.")
end

struct ActivationMode{T} end
const NoActivation = ActivationMode{:NoActivation}()
const DefaultActivation = ActivationMode{:Default}()

# Flux ####################################################

const DenseLayers = NTuple{N,Dense} where {N}
const InputDense = Tuple{<:Function,Vararg{Dense,N}} where {N}

function (a::Dense)(x::AbstractVecOrMat, s::ActivationMode)
    s === DefaultActivation && return a(x)
    W, b = a.weight, a.bias
    if b === Flux.Zeros()
        return W*x
    else
        return W*x .+ b
    end
end

function _applychain(fs::DenseLayers, x)
    y = Flux.applychain(front(fs),x)
    last(fs)(y, NoActivation)
end
_applychain(fs::InputDense, x) = _applychain(tail(fs), first(fs)(x))
_applychain(::Tuple{<:Function}, x) = x

chainGrad(::Tuple{}, x) = LinearAlgebra.I
function chainGrad(fs::T, x)  where {T<:Union{DenseLayers,InputDense}}
    y = _applychain(fs,x)
    jacobian(last(fs), y) * chainGrad(front(fs), x)
end

function jacobian(m::Chain{T}, x) where {T<:Union{DenseLayers,InputDense}}
    chainGrad(m.layers,x)
end
function jacobian(c::Chain, x)
    fs = front(c.layers)
    J = chainGrad(fs, x)
    fout = typeof(last(c.layers))
    if fout === typeof(posdef)
        y = Flux.applychain(fs, x)
        L = ltril(y)
        ∂L∂x = map(eachcol(J)) do col
            dL = ltril(col)
            (L * dL') + (dL * L')
        end
    elseif fout === typeof(skewsym)
        ∂L∂x = map(eachcol(J)) do col
            dL = skewsym(col)
        end
    else
        error("Analytical Jacobian not supported for this Chain type. Use AD.")
    end
end

linear(m::Chain, x) = extraChain(Tuple(m.layers), x)
function extraChain(fs::Tuple, x)
    res = first(fs)(x, NoActivation)
    σ = first(fs).σ
    return (res, extraChain(Base.tail(fs), σ.(res))...)
end
extraChain(::Tuple{}, x) = ()

makeodd(m::Chain) = Chain(square, m.layers...)
makeposdef(m::Chain) = Chain(m.layers..., posdef)
makeskewsym(m::Chain) = Chain(m.layers..., skewsym)

# # DiffEqFlux ##############################################

const FastDenseLayers = NTuple{N,FastDense} where {N}
const InputFastDense = Tuple{<:Function,Vararg{FastDense,N}} where {N}

function (a::FastDense)(x::AbstractVecOrMat, p, s::ActivationMode)
    s === DefaultActivation && return a(x, p)
    W, b = param2Wb(a, p)
    if a.bias
        return W*x .+ b
    else
        return W*x
    end
end
function _applychain(fs::FastDenseLayers, x, p)
    ps = paramslice(fs,p)
    psvec = reduce(vcat,front(ps))
    y = DiffEqFlux.applychain(front(fs),x,psvec)
    last(fs)(y, last(ps), NoActivation)
end
_applychain(fs::Tuple{FastDense}, x, p) = last(fs)(x, p, NoActivation)
_applychain(fs::InputFastDense, x, p) = _applychain(tail(fs), first(fs)(x), p)
_applychain(::Tuple{<:Function}, x, ::Any) = x

function chainGrad(fs::InputFastDense, x, p)
    y = _applychain(fs,x,p)
    ps = last(paramslice(fs,p))
    return jacobian(last(fs), y, ps) * chainGrad(front(fs), x, p[1:end-length(ps)])
end
chainGrad(::Tuple{}, ::Any, ::Any) = LinearAlgebra.I

function jacobian(m::FastChain{T}, x, p) where {T<:Union{FastDenseLayers,InputFastDense}}
    chainGrad(m.layers,x,p)
end
function jacobian(c::FastChain, x, p)
    fs = front(c.layers)
    J = chainGrad(fs, x, p)
    fout = typeof(last(c.layers))
    if fout === typeof(posdef)
        y = DiffEqFlux.applychain(fs, x, p)
        L = ltril(y)
        ∂L∂x = map(eachcol(J)) do col
            dL = ltril(col)
            (L * dL') + (dL * L')
        end
    elseif fout === typeof(skewsym)
        ∂L∂x = map(eachcol(J)) do col
            dL = skewsym(col)
        end
    else
        error("Analytical Jacobian not supported for this Chain type. Use AD.")
    end
end

paramslice(c::FastChain{<:FastDenseLayers}, p::Vector) = paramslice(c.layers,p)
function paramslice(fs::Tuple, p)   # layers fs must be in increasing order
    f = first(fs)
    N = DiffEqFlux.paramlength(f)
    res = @view(p[1:N])
    return (res, paramslice(tail(fs), p[N+1:end])...)
end
paramslice(::Tuple{}, ::Any) = ()

function param2Wb(f::FastDense, p)
    np = DiffEqFlux.paramlength(f)
    W = @view p[reshape(1:(f.out*f.in),f.out,f.in)]
    b = @view p[(f.out*f.in + 1):np]
    W, b
end

makeodd(m::FastChain) = FastChain(square, m.layers...)
makeposdef(m::FastChain) = FastChain(m.layers..., posdef)
makeskewsym(m::FastChain) = FastChain(m.layers..., skewsym)

ext: jl


text:
"""
See [`SimpleFileServer.serve`](@ref) to run the file server.
"""
module SimpleFileServer

using Deno_jll

export serve

"""
    serve(dir="."; <keyword arguments>)

Run a simple HTTP file server, serving files from the given local directory.

The server is running on [deno](https://deno.land/), using the [`file_server` standard library](https://deno.land/manual/examples/file_server).

# Arguments
- `cors::Bool=false`: Enable CORS via the "Access-Control-Allow-Origin" header.
- `port::Integer=4507`: Set port.
- `host="0.0.0.0"`: Hostname (default is 0.0.0.0, another common choice is 127.0.0.1).
- `dir_listing::Bool=true`: Show a directory listing for folders with no `index.html` file.
- `dotfiles::Bool=true`: Also show dotfiles in directory listings.
- `tls_certificate_file::Union{Nothing,String}=nothing`: TLS certificate file (enables TLS).
- `tls_key_file::Union{Nothing,String}=nothing`: TLS key file (enables TLS).

All TLS options are required when one is provided.
"""
function serve(dir::String=".";
    port::Integer=4507,
    host="0.0.0.0",
    cors::Bool=false,
    dir_listing::Bool=true,
    dotfiles::Bool=true,
    tls_certificate_file::Union{Nothing,String}=nothing,
    tls_key_file::Union{Nothing,String}=nothing,
)

    args = []

    push!(args, "--port", string(port))
    push!(args, "--host", string(host))

    if cors
        push!(args, "--cors")
    end
    if !dir_listing
        push!(args, "--no-dir-listing")
    end
    if !dotfiles
        push!(args, "--no-dotfiles")
    end

    if tls_certificate_file !== nothing
        push!(args, "--cert", tls_certificate_file)
    end
    if tls_key_file !== nothing
        push!(args, "--key", tls_key_file)
    end

    run(`$(deno()) run --allow-net --allow-read https://deno.land/std@0.97.0/http/file_server.ts $(dir) $(args)`)
end
end


text:
struct Woodbury{T,AType,UType,VType,CType,CpType} <: AbstractWoodbury{T}
    A::AType
    U::UType
    C::CType
    Cp::CpType
    V::VType
    tmpN1::Vector{T}
    tmpN2::Vector{T}
    tmpk1::Vector{T}
    tmpk2::Vector{T}
end

"""
    W = Woodbury(A, U, C, V)

Represent a matrix `W = A + UCV`.
Equations `Wx = b` will be solved using the
[Woodbury matrix identity](https://en.wikipedia.org/wiki/Woodbury_matrix_identity).

If your main goal is to solve equations, it's often advantageous to supply
`A` as a factorization (e.g., `Woodbury(lu(A), U, C, V)`).

See also [SymWoodbury](@ref).
"""
function Woodbury(A, U::AbstractMatrix, C, V::AbstractMatrix)
    @noinline throwdmm1(U, V, A) = throw(DimensionMismatch("Sizes of U ($(size(U))) and/or V ($(size(V))) are inconsistent with A ($(size(A)))"))
    @noinline throwdmm2(k) = throw(DimensionMismatch("C should be $(k)x$(k)"))

    N = size(A, 1)
    k = size(U, 2)
    if size(A, 2) != N || size(U, 1) != N || size(V, 1) != k || size(V, 2) != N
        throwdmm1(U, V, A)
    end
    if k > 1
        if size(C, 1) != k || size(C, 2) != k
            throwdmm2(k)
        end
    end
    Cp = safeinv(safeinv(C) .+ V*(A\U))
    # temporary space for allocation-free solver (vector RHS only)
    T = typeof(float(zero(eltype(A)) * zero(eltype(U)) * zero(eltype(C)) * zero(eltype(V))))
    tmpN1 = Vector{T}(undef, N)
    tmpN2 = Vector{T}(undef, N)
    tmpk1 = Vector{T}(undef, k)
    tmpk2 = Vector{T}(undef, k)

    Woodbury(A, U, C, Cp, V, tmpN1, tmpN2, tmpk1, tmpk2)
end

Woodbury(A, U::AbstractVector{T}, C, V::AbstractMatrix{T}) where {T} = Woodbury(A, reshape(U, length(U), 1), C, V)

# Woodbury(A, U::AbstractVector, C, V::Adjoint) = Woodbury(A, U, C, Matrix(V))

function show(io::IO, W::Woodbury)
    println(io, "Woodbury factorization:\nA:")
    show(io, MIME("text/plain"), W.A)
    print(io, "\nU:\n")
    Base.print_matrix(IOContext(io, :compact=>true), W.U)
    if isa(W.C, Matrix)
        print(io, "\nC:\n")
        Base.print_matrix(IOContext(io, :compact=>true), W.C)
    else
        print(io, "\nC: ", W.C)
    end
    print(io, "\nV:\n")
    Base.print_matrix(IOContext(io, :compact=>true), W.V)
end

copy(W::Woodbury) = Woodbury(W.A, W.U, W.C, W.V)

function inv(W::AbstractWoodbury)
    A′ = inv(W.A)
    U′ = W.A \ W.U
    C′ = -W.Cp
    V′ = W.V*A′
    return Woodbury(A′, U′, C′, V′)
end

# multiplying everything out would be O(m^2 k + m k^2) operations
# best I could think of was O(m^2 k) operations
diag(W::AbstractWoodbury) = diag(W.A) + vec(sum(W.U .* (W.C * W.V)', dims=2))

+(W::Woodbury, X::Woodbury)       = Woodbury(W.A + X.A, [W.U X.U],
                                             cat(W.C, X.C; dims=(1,2)), [W.V; X.V])
*(α::Real, W::Woodbury)           = Woodbury(α*W.A, W.U, α*W.C, W.V)
+(W::Woodbury, M::AbstractMatrix) = Woodbury(W.A + M, W.U, W.C, W.V)

adjoint(W::Woodbury) = Woodbury(adjoint(W.A), adjoint(W.V), adjoint(W.C), adjoint(W.U))
transpose(W::Woodbury) = Woodbury(transpose(W.A), transpose(W.V), transpose(W.C), transpose(W.U))


text:
function BicubicInterpolator()
    return BicubicInterpolator(())
end

function interpolate(obj::BicubicInterpolator, arg0::Vector{jdouble}, arg1::Vector{jdouble}, arg2::Vector{Vector{jdouble}})
    return jcall(obj, "interpolate", BicubicInterpolatingFunction, (Vector{jdouble}, Vector{jdouble}, Vector{Vector{jdouble}}), arg0, arg1, arg2)
end


text:
"""
    function generate_bracket(
        teams::Vector{Team},
        drawsize::Int64
    )::Vector{Any}

Generates a tournament bracket.

Arguments
=========
- `teams`:
    a vector of teams participating in the tournament.
- `drawsize`
    tournament draw size.
"""

function generate_bracket(
    teams::Vector{Team},
    drawsize::Int64
)::Vector{Any}
    @assert (drawsize != 0) && (drawsize-1 != 0)
    matches = Int(drawsize / 2)
    bracket = teams
    while matches > 1
        newbracket = Vector{Any}([[] for _ in 1:matches])
        for m in 1:matches
            push!(newbracket[m], bracket[m])
            push!(newbracket[m], bracket[end-m+1])
        end
        matches = Int(matches/2)
        bracket = newbracket
    end
    return bracket
end

"""
    function execute_bracket(
        current_bracket::Union{Vector{Any}, Vector{Team}, Team},
        tournament::Tournament,
        round::Int64,
        strength_method::String
    )::Team

Simulates a tournament bracket.

Arguments
=========
- `current_bracket`:
    a vector of vectors representing the tournament bracket.
- `tournament`:
    the current tournament being played.
- `round`:
    the current round being played.
- `strength_method`:
    method used for determing probability team 1 beats team 2.
"""

function execute_bracket(
    current_bracket::Union{Vector{Any}, Vector{Team}, Team},
    tournament::Tournament,
    round::Int64,
    strength_method::String
)::Team
    if isa(current_bracket, Vector) == false
        return current_bracket
    end
    team1 = execute_bracket(
        current_bracket[1],
        tournament,
        round * 2,
        strength_method
    )
    team2 = execute_bracket(
        current_bracket[2],
        tournament,
        round * 2,
        strength_method
    )
    winning_team = play_doubles_match(
        team1,
        team2,
        tournament,
        round,
        strength_method
    )
    return winning_team
end

max_stars_repo_name: Danial-Hussain/Tennis-Ranking

max_issues_repo_name: Danial-Hussain/Tennis-Ranking

max_forks_repo_name: Danial-Hussain/Tennis-Ranking


text:
const XFER_BUF_SIZE = 32767
const MB = 1048576

struct SFTPPath
    host::SSHHost
    path::String
end

function Base.show(io::IO, p::SFTPPath)
    print(io, "SFTPPath(\"$(p.host.address):$(p.host.port):$(p.path)\")")
end

SFTPPath(hostname::String, path) = SFTPPath(SSHHost(hostname), path)
SFTPPath(host::SSHHost, path) = SFTPPath(host, path)
Base.joinpath(p1::SFTPPath, p2::SFTPPath) = @assert p1.host === p2.host && return SFTPPath(p1.host, joinpath(p1.path, p2.path))
Base.joinpath(p1::SFTPPath, p2::AbstractString) = SFTPPath(p1.host, joinpath(p1.path, p2))
Base.basename(p::SFTPPath) = basename(p.path)
Base.endswith(p::SFTPPath, ending) = endswith(p.path, ending)
Base.isless(a::SFTPPath, b::SFTPPath) = isless(a.host, b.host) && isless(a.path, b.path)
Base.convert(T::Type{<:String}, p::SFTPPath) = p.path

struct SSHDirDescriptor
    handle::Ptr{Cvoid}
    session::SFTPSession
end

Base.convert(::Type{Ptr{LIBSSH2_SFTP_HANDLE}}, sess::SSHDirDescriptor) = sess.handle
Base.unsafe_convert(::Type{Ptr{LIBSSH2_SFTP_HANDLE}}, sess::SSHDirDescriptor) = convert(Ptr{LIBSSH2_SFTP_HANDLE}, sess.handle)

function sftp_open_dir(session::SFTPSession, path::AbstractString)
    dir_descriptor = libssh2_sftp_opendir(session, path)
    return SSHDirDescriptor(dir_descriptor, session)
end

sftp_open_dir(session::SFTPSession, path::SFTPPath) = sftp_open_dir(session, path.path)

Base.close(dir::SSHDirDescriptor) = libssh2_sftp_close_handle(dir)

Base.readdir(path::SFTPPath; kwargs...) =
    readdir(current_sftp_session(), path; kwargs...)

function Base.readdir(session::SFTPSession, path::SFTPPath; join::Bool=false, sort::Bool=true)
    dir = sftp_open_dir(session, path)

    entries = SFTPPath[]
    is_eof = false
    while true
        is_eof, name, info, attr = libssh2_sftp_readdir(dir, 2^8)
        is_eof && break
        (name == "." || name == "..") && continue
        fname = SFTPPath(path.host, name)
        push!(entries, join ? joinpath(path, fname) : fname)
    end

    sort && sort!(entries)

    close(dir)

    return entries
end

##
mutable struct SSHFileDescriptor <: IO
    handle::Ptr{Cvoid}
    session::SFTPSession
    eof::Bool
end

Base.convert(::Type{Ptr{LIBSSH2_SFTP_HANDLE}}, file::SSHFileDescriptor) = file.handle
Base.unsafe_convert(::Type{Ptr{LIBSSH2_SFTP_HANDLE}}, file::SSHFileDescriptor) = convert(Ptr{LIBSSH2_SFTP_HANDLE}, file.handle)

Base.close(file::SSHFileDescriptor) = libssh2_sftp_close_handle(file)
Base.eof(file::SSHFileDescriptor) = file.eof

Base.open(path::SFTPPath, args...; kwargs...) =
    Base.open(current_sftp_session(), path, args...; kwargs...)

function Base.open(session::SFTPSession, path::SFTPPath, mode::AbstractString; lock = true)
    mode == "r"  ? Base.open(session, path, lock = lock, read = true)                  :
    mode == "r+" ? Base.open(session, path, lock = lock, read = true, write = true)    :
    mode == "w"  ? Base.open(session, path, lock = lock, truncate = true)              :
    mode == "w+" ? Base.open(session, path, lock = lock, truncate = true, read = true) :
    mode == "a"  ? Base.open(session, path, lock = lock, append = true)                :
    mode == "a+" ? Base.open(session, path, lock = lock, append = true, read = true)   :
    throw(ArgumentError("invalid open mode: $mode"))
end

function Base.open(session::SFTPSession, path::SFTPPath; lock = true,
    read     :: Union{Bool,Nothing} = nothing,
    write    :: Union{Bool,Nothing} = nothing,
    create   :: Union{Bool,Nothing} = nothing,
    truncate :: Union{Bool,Nothing} = nothing,
    append   :: Union{Bool,Nothing} = nothing,
)
    flags = UInt32(0)

    flags = flags | ( read     == true ? LIBSSH2_FXF_READ   : 0x0 )
    flags = flags | ( write    == true ? LIBSSH2_FXF_WRITE  : 0x0 )
    flags = flags | ( create   == true ? LIBSSH2_FXF_CREAT  : 0x0 )
    flags = flags | ( truncate == true ? LIBSSH2_FXF_TRUNC  : 0x0 )
    flags = flags | ( append   == true ? LIBSSH2_FXF_APPEND : 0x0 )

    permissions = LIBSSH2_SFTP_S_IRUSR | LIBSSH2_SFTP_S_IWUSR | LIBSSH2_SFTP_S_IRGRP | LIBSSH2_SFTP_S_IROTH

    f_handle = libssh2_sftp_open(session, path.path, flags, permissions)

    return SSHFileDescriptor(f_handle, session, false)
end

Base.read(file::SFTPPath, args...) = read(current_sftp_session(), file, args...)
Base.read(session::SFTPSession, file::SFTPPath, args...) = open(io->read(io, args...), session, file)

Base.read(file::SFTPPath, ::Type{T}) where {T} = read(current_sftp_session(), file, T)
Base.read(session::SFTPSession, file::SFTPPath, ::Type{T}) where {T} = open(io->read(io, T), session, file)

Base.filesize(file::SFTPPath) = filesize(current_sftp_session(), file)
Base.filesize(session::SFTPSession, file::SFTPPath) = begin
    f = open(session, file)
    fs = filesize(f)
    close(f)
    return fs
end
Base.filesize(file::SSHFileDescriptor) = filesize(libssh2_sftp_fstat(file))

function unsafe_read_atmost(file::SSHFileDescriptor, p::Ptr{UInt8}, n::Integer)
    nr = libssh2_sftp_read(file, p, n)
    if nr == 0
    	file.eof = true
    end
    @assert nr >= 0
    return nr
end

Base.readbytes_all!(s::SSHFileDescriptor, b::Array{UInt8}, nb) = begin
    GC.@preserve b begin
        ret = unsafe_read_atmost(s, pointer(b), length(b))
    end
    return ret
end

function Base.unsafe_read(file::SSHFileDescriptor, p::Ptr{UInt8}, n::UInt)
    nr = libssh2_sftp_read(file, p, n)
    if nr != n
    	file.eof = true
        throw(EOFError())
    end

    return nr
end

function Base.read(file::SSHFileDescriptor, ::Type{UInt8})
	mem = Ref(UInt8(0))
	b = unsafe_read(file, pointer_from_objref(mem), 1)
    return mem[]
end

function Base.readbytes!(s::SSHFileDescriptor, b::AbstractArray{UInt8}, nb=length(b))
    Base.require_one_based_indexing(b)
    olb = lb = length(b)
    nr = 0

    while !eof(s)
		if nr == lb
			lb = nr * 2
			resize!(b, lb)
		end

	    b_read = unsafe_read_atmost(s, pointer(b, nr+1), lb-nr)
	    nr += b_read
	end
    if nr != nb
    	resize!(b, nr)
    end
    return nr
end

function Base.unsafe_write(file::SSHFileDescriptor, p::Ptr{UInt8}, n::UInt)
    written::Int = 0
    written = libssh2_sftp_write(file, p, n)
    return written
end

function Base.write(file::SSHFileDescriptor, b::UInt8)
    #iswritable(s) || throw(ArgumentError("write failed, IOStream is not writeable"))
    mem = Ref(b)
    GC.@preserve mem begin
        retcode = unsafe_write(file, pointer_from_objref(mem), 1)
    end
    @assert retcode == 1
    return retcode
end

Base.seek(file::SSHFileDescriptor, pos::Integer) = libssh2_sftp_seek64(file, pos)
Base.position(file::SSHFileDescriptor) = convert(Int, libssh2_sftp_tell64(file))
Base.skip(file::SSHFileDescriptor, delta::Integer) = seek(file, position(file) + delta)


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

"""
The ArrowTypes module provides the [`ArrowTypes.Arrowtype`](@ref) interface trait that objects can define
in order to signal how they should be serialized in the arrow format.
"""
module ArrowTypes

using UUIDs

export ArrowKind, NullKind, PrimitiveKind, BoolKind, ListKind, FixedSizeListKind, MapKind, StructKind, UnionKind, DictEncodedKind, toarrow, arrowname, fromarrow, ToArrow

"""
    ArrowTypes.ArrowKind(T)

For a give type `T`, define it's "arrow type kind", or the general category of arrow types it should be treated as. Must be one of:
  * [`ArrowTypes.NullKind`](@ref): `Missing` is the only type defined as `NullKind`
  * [`ArrowTypes.PrimitiveKind`](@ref): `<:Integer`, `<:AbstractFloat`, along with `Arrow.Decimal`, and the various `Arrow.ArrowTimeType` subtypes
  * [`ArrowTypes.BoolKind`](@ref): only `Bool`
  * [`ArrowTypes.ListKind`](@ref): any `AbstractString` or `AbstractArray`
  * [`ArrowTypes.FixedSizeList`](@ref): `NTuple{N, T}`
  * [`ArrowTypes.MapKind`](@ref): any `AbstractDict`
  * [`ArrowTypes.StructKind`](@ref): any `NamedTuple` or plain struct (mutable or otherwise)
  * [`ArrowTypes.UnionKind`](@ref): any `Union`
  * [`ArrowTypes.DictEncodedKind`](@ref): array types that implement the `DataAPI.refarray` interface

The list of `ArrowKind`s listed above translate to different ways to physically store data as supported by the arrow data format.
See the docs for each for an idea of whether they might be an appropriate fit for a custom type.
Note that custom types need to satisfy any additional "interface methods" as required by the various `ArrowKind`
types. By default, if a type in julia is declared like `primitive type ...` it is considered a `PrimitiveKind`
and if `struct` or `mutable struct` it's considered a `StructKind`. Also note that types will rarely need to define `ArrowKind`;
much more common is to define `ArrowType(T)` and `toarrow(x::T)` to transform `T` to a natively supported arrow type, which will
already have its `ArrowKind` defined.
"""
abstract type ArrowKind end

ArrowKind(x::T) where {T} = ArrowKind(T)
ArrowKind(::Type{T}) where {T} = isprimitivetype(T) ? PrimitiveKind() : StructKind()

"""
    ArrowTypes.ArrowType(T) = S

Interface method to define the natively supported arrow type `S` that a given type `T` should be converted to before serializing.
Useful when a custom type wants a "serialization hook" or otherwise needs to be transformed/converted into a natively
supported arrow type for serialization. If a type defines `ArrowType`, it must also define a corresponding
[`ArrowTypes.toarrow(x::T)`](@ref) method which does the actual conversion from `T` to `S`.
Note that custom structs defined like `struct T` or `mutable struct T` are natively supported in serialization, so unless
_additional_ transformation/customization is desired, a custom type `T` can serialize with no `ArrowType` definition (by default,
each field of a struct is serialized, using the results of `fieldnames(T)` and `getfield(x, i)`).
Note that defining these methods only deal with custom _serialization_ to the arrow format; to be able to _deserialize_ custom
types at all, see the docs for [`ArrowTypes.arrowname`](@ref), [`ArrowTypes.arrowmetadata`](@ref), [`ArrowTypes.JuliaType`](@ref),
and [`ArrowTypes.fromarrow`](@ref).
"""
function ArrowType end
ArrowType(::Type{T}) where {T} = T
ArrowType(::Type{Any}) = Any
ArrowType(::Type{Union{Missing, T}}) where {T} = Union{Missing, ArrowType(T)}
ArrowType(::Type{Missing}) = Missing

"""
    ArrowTypes.toarrow(x::T) => S

Interface method to perform the actual conversion from an object `x` of type `T` to the type `S`. `T` and `S` must match the
types used when defining `ArrowTypes.ArrowType(::Type{T}) = S`. Hence, `S` is the natively supported arrow type that `T`
desires to convert to to enable serialization. See [`ArrowTypes.ArrowType`](@ref) docs for more details.
This enables custom objects to be serialized as a natively supported arrow data type.
"""
function toarrow end
toarrow(x) = x

"""
    ArrowTypes.arrowname(T) = Symbol(name)

Interface method to define the logical type "label" for a custom Julia type `T`. Names will be global for an entire arrow dataset,
and conventionally, custom types will just use their type name along with a Julia-specific prefix; for example, for a custom type
`Foo`, I would define `ArrowTypes.arrowname(::Type{Foo}) = Symbol("JuliaLang.Foo")`. This ensures other language implementations
won't get confused and are safe to ignore the logical type label.
When arrow stores non-native data, it must still be _stored_ as a native data type, but can have type metadata tied to the data that
labels the original _logical_ type it originated from. This enables the conversion of native data back to the logical type when
deserializing, as long as the deserializer has the same definitions when the data was serialized. Namely, the current Julia
session will need the appropriate [`ArrowTypes.JuliaType`](@ref) and [`ArrowTypes.fromarrow`](@ref) definitions in order to know
how to convert the native data to the original logical type. See the docs for those interface methods in order to ensure a complete
implementation. Also see the accompanying [`ArrowTypes.arrowmetadata`](@ref) docs around providing additional metadata about a custom
logical type that may be necessary to fully re-create a Julia type (e.g. non-field-based type parameters).
"""
function arrowname end
const EMPTY_SYMBOL = Symbol()
arrowname(T) = EMPTY_SYMBOL
hasarrowname(T) = arrowname(T) !== EMPTY_SYMBOL

"""
    ArrowTypes.arrowmetadata(T) => String

Interface method to provide additional logical type metadata when serializing extension types. [`ArrowTypes.arrowname`](@ref)
provides the logical type _name_, which may be all that's needed to return a proper Julia type from [`ArrowTypes.JuliaType`](@ref),
but some custom types may, for example have type parameters that aren't inferred/based on fields. In order to fully recreate these
kinds of types when deserializing, these type parameters can be stored by defining `ArrowTypes.arrowmetadata(::Type{T}) = "type_param"`.
This will then be available to access by overloading `ArrowTypes.JuliaType(::Val{Symbol(name)}, S, arrowmetadata::String)`.
"""
function arrowmetadata end
arrowmetadata(T) = ""

"""
    ArrowTypes.JuliaType(::Val{Symbol(name)}, ::Type{S}, arrowmetadata::String) = T

Interface method to define the custom Julia logical type `T` that a serialized metadata label should be converted to when
deserializing. When reading arrow data, and a logical type label is encountered for a column, it will call
`ArrowTypes.JuliaType(Val(Symbol(name)), S, arrowmetadata)` to see if a Julia type has been "registered" for deserialization. The `name`
used when defining the method *must* correspond to the same `name` when defining `ArrowTypes.arrowname(::Type{T}) = Symbol(name)`.
The use of `Val(Symbol(...))` is to allow overloading a method on a specific logical type label. The `S` 2nd argument passed to
`JuliaType` is the native arrow serialized type. This can be useful for parametric Julia types that wish to correctly parameterize
their custom type based on what was serialized. The 3rd argument `arrowmetadata` is any metadata that was stored when the logical
type was serialized as the result of calling `ArrowTypes.arrowmetadata(T)`. Note the 2nd and 3rd arguments are optional when
overloading if unneeded.
When defining [`ArrowTypes.arrowname`](@ref) and `ArrowTypes.JuliaType`, you may also want to implement [`ArrowTypes.fromarrow`]
in order to customize how a custom type `T` should be constructed from the native arrow data type. See its docs for more details.
"""
function JuliaType end
JuliaType(val) = nothing
JuliaType(val, S) = JuliaType(val)
JuliaType(val, S, meta) = JuliaType(val, S)

"""
    ArrowTypes.fromarrow(::Type{T}, x::S) => T

Interface method that provides a "deserialization hook" for a custom type `T` to be constructed from the native arrow type `S`.
The `T` and `S` types must correspond to the definitions used in `ArrowTypes.ArrowType(::Type{T}) = S`. This is a paired method
with [`ArrowTypes.toarrow`](@ref).

The default definition is `ArrowTypes.fromarrow(::Type{T}, x) = T(x)`, so if that works for a custom type already, no additional
overload is necessary.
A few `ArrowKind`s have/allow slightly more custom overloads for their `fromarrow` methods:
  * `ListKind{true}`: for `String` types, they may overload `fromarrow(::Type{T}, ptr::Ptr{UInt8}, len::Int) = ...` to avoid
     materializing a `String`
  * `StructKind`: must overload `fromarrow(::Type{T}, x...)` where individual fields are passed as separate
     positional arguments; so if my custom type `Interval` has two fields `first` and `last`, then I'd overload like
     `ArrowTypes.fromarrow(::Type{Interval}, first, last) = ...`. Note the default implementation is
     `ArrowTypes.fromarrow(::Type{T}, x...) = T(x...)`, so if your type already accepts all arguments in a constructor
     no additional `fromarrow` method should be necessary (default struct constructors have this behavior).
"""
function fromarrow end
fromarrow(::Type{T}, x::T) where {T} = x
fromarrow(::Type{T}, x...) where {T} = T(x...)
fromarrow(::Type{Union{Missing, T}}, ::Missing) where {T} = missing
fromarrow(::Type{Union{Missing, T}}, x::T) where {T} = x
fromarrow(::Type{Union{Missing, T}}, x) where {T} = fromarrow(T, x)

"NullKind data is actually not physically stored since the data is constant; just the length is needed"
struct NullKind <: ArrowKind end

ArrowKind(::Type{Missing}) = NullKind()
ArrowKind(::Type{Nothing}) = NullKind()
ArrowType(::Type{Nothing}) = Missing
toarrow(::Nothing) = missing
const NOTHING = Symbol("JuliaLang.Nothing")
arrowname(::Type{Nothing}) = NOTHING
JuliaType(::Val{NOTHING}) = Nothing
fromarrow(::Type{Nothing}, ::Missing) = nothing

"PrimitiveKind data is stored as plain bits in a single contiguous buffer"
struct PrimitiveKind <: ArrowKind end

ArrowKind(::Type{<:Integer}) = PrimitiveKind()
ArrowKind(::Type{<:AbstractFloat}) = PrimitiveKind()

ArrowType(::Type{Char}) = UInt32
toarrow(x::Char) = convert(UInt32, x)
const CHAR = Symbol("JuliaLang.Char")
arrowname(::Type{Char}) = CHAR
JuliaType(::Val{CHAR}) = Char
fromarrow(::Type{Char}, x::UInt32) = Char(x)

"BoolKind data is stored with values packed down to individual bits; so instead of a traditional Bool being 1 byte/8 bits, 8 Bool values would be packed into a single byte"
struct BoolKind <: ArrowKind end
ArrowKind(::Type{Bool}) = BoolKind()

"ListKind data are stored in two separate buffers; one buffer contains all the original data elements flattened into one long buffer; the 2nd buffer contains an offset into the 1st buffer for how many elements make up the original array element"
struct ListKind{stringtype} <: ArrowKind end

ListKind() = ListKind{false}()
isstringtype(::ListKind{stringtype}) where {stringtype} = stringtype
isstringtype(::Type{ListKind{stringtype}}) where {stringtype} = stringtype

ArrowKind(::Type{<:AbstractString}) = ListKind{true}()

fromarrow(::Type{T}, ptr::Ptr{UInt8}, len::Int) where {T} = fromarrow(T, unsafe_string(ptr, len))

ArrowType(::Type{Symbol}) = String
toarrow(x::Symbol) = String(x)
const SYMBOL = Symbol("JuliaLang.Symbol")
arrowname(::Type{Symbol}) = SYMBOL
JuliaType(::Val{SYMBOL}) = Symbol
_symbol(ptr, len) = ccall(:jl_symbol_n, Ref{Symbol}, (Ptr{UInt8}, Int), ptr, len)
fromarrow(::Type{Symbol}, ptr::Ptr{UInt8}, len::Int) = _symbol(ptr, len)

ArrowKind(::Type{<:AbstractArray}) = ListKind()
fromarrow(::Type{A}, x::AbstractVector{T}) where {A <: AbstractVector{T}} where {T} = x
ArrowKind(::Type{<:AbstractSet}) = ListKind()
ArrowType(::Type{T}) where {T <: AbstractSet{S}} where {S} = Vector{S}
toarrow(x::AbstractSet) = collect(x)
const SET = Symbol("JuliaLang.Set")
arrowname(::Type{<:AbstractSet}) = SET
JuliaType(::Val{SET}, ::Type{T}) where {T <: AbstractVector{S}} where {S} = Set{S}
fromarrow(::Type{T}, x) where {T <: AbstractSet} = T(x)

"FixedSizeListKind data are stored in a single contiguous buffer; individual elements can be computed based on the fixed size of the lists"
struct FixedSizeListKind{N, T} <: ArrowKind end
gettype(::FixedSizeListKind{N, T}) where {N, T} = T
getsize(::FixedSizeListKind{N, T}) where {N, T} = N

ArrowKind(::Type{NTuple{N, T}}) where {N, T} = FixedSizeListKind{N, T}()

ArrowKind(::Type{UUID}) = FixedSizeListKind{16, UInt8}()
ArrowType(::Type{UUID}) = NTuple{16, UInt8}
toarrow(x::UUID) = _cast(NTuple{16, UInt8}, x.value)
const UUIDSYMBOL = Symbol("JuliaLang.UUID")
arrowname(::Type{UUID}) = UUIDSYMBOL
JuliaType(::Val{UUIDSYMBOL}) = UUID
fromarrow(::Type{UUID}, x::NTuple{16, UInt8}) = UUID(_cast(UInt128, x))
# for backwards compatibility
arrowconvert(::Type{UUID}, u::NamedTuple{(:value,),Tuple{UInt128}}) = UUID(u.value)
arrowconvert(::Type{UUID}, u::UInt128) = UUID(u)

function _cast(::Type{Y}, x)::Y where {Y}
    y = Ref{Y}()
    _unsafe_cast!(y, Ref(x), 1)
    return y[]
end

function _unsafe_cast!(y::Ref{Y}, x::Ref, n::Integer) where {Y}
    X = eltype(x)
    GC.@preserve x y begin
        ptr_x = Base.unsafe_convert(Ptr{X}, x)
        ptr_y = Base.unsafe_convert(Ptr{Y}, y)
        unsafe_copyto!(Ptr{X}(ptr_y), ptr_x, n)
    end
    return y
end

"StructKind data are stored in separate buffers for each field of the struct"
struct StructKind <: ArrowKind end

ArrowKind(::Type{<:NamedTuple}) = StructKind()

fromarrow(::Type{NamedTuple{names, types}}, x::NamedTuple{names, types}) where {names, types <: Tuple} = x
fromarrow(::Type{T}, x::NamedTuple) where {T} = fromarrow(T, Tuple(x)...)

ArrowKind(::Type{<:Tuple}) = StructKind()
const TUPLE = Symbol("JuliaLang.Tuple")
# needed to disambiguate the FixedSizeList case for NTuple
arrowname(::Type{NTuple{N, T}}) where {N, T} = EMPTY_SYMBOL
arrowname(::Type{T}) where {T <: Tuple} = TUPLE
JuliaType(::Val{TUPLE}, ::Type{NamedTuple{names, types}}) where {names, types <: Tuple} = types
fromarrow(::Type{T}, x::NamedTuple) where {T <: Tuple} = Tuple(x)

"MapKind data are stored similarly to ListKind, where elements are flattened, and a 2nd offsets buffer contains the individual list element length data"
struct MapKind <: ArrowKind end

ArrowKind(::Type{<:AbstractDict}) = MapKind()

"UnionKind data are stored either in a separate, compacted buffer for each union type (dense), or in full-length buffers for each union type (sparse)"
struct UnionKind <: ArrowKind end

ArrowKind(::Union) = UnionKind()

"DictEncodedKind store a small pool of unique values in one buffer, with a full-length buffer of integer offsets into the small value pool"
struct DictEncodedKind <: ArrowKind end

"""
There are a couple places when writing arrow buffers where
we need to write a "dummy" value; it doesn't really matter
what we write, but we need to write something of a specific
type. So each supported writing type needs to define `default`.
"""
function default end

default(T) = zero(T)
default(::Type{Symbol}) = Symbol()
default(::Type{Char}) = '\0'
default(::Type{<:AbstractString}) = ""
default(::Type{Any}) = nothing
default(::Type{Union{T, Missing}}) where {T} = default(T)

function default(::Type{A}) where {A <: AbstractVector{T}} where {T}
    a = similar(A, 1)
    a[1] = default(T)
    return a
end

default(::Type{NTuple{N, T}}) where {N, T} = ntuple(i -> default(T), N)
default(::Type{T}) where {T <: Tuple} = Tuple(default(fieldtype(T, i)) for i = 1:fieldcount(T))
default(::Type{T}) where {T <: AbstractDict} = T()
default(::Type{NamedTuple{names, types}}) where {names, types} = NamedTuple{names}(Tuple(default(fieldtype(types, i)) for i = 1:length(names)))

function promoteunion(T, S)
    new = promote_type(T, S)
    return isabstracttype(new) ? Union{T, S} : new
end

# lazily call toarrow(x) on getindex for each x in data
struct ToArrow{T, A} <: AbstractVector{T}
    data::A
end

concrete_or_concreteunion(T) = isconcretetype(T) ||
    (T isa Union && concrete_or_concreteunion(T.a) && concrete_or_concreteunion(T.b))

function ToArrow(x::A) where {A}
    S = eltype(A)
    T = ArrowType(S)
    if S === T && concrete_or_concreteunion(S)
        return x
    elseif !concrete_or_concreteunion(T)
        # arrow needs concrete types, so try to find a concrete common type, preferring unions
        if isempty(x)
            return Missing[]
        end
        T = typeof(toarrow(x[1]))
        for i = 2:length(x)
            @inbounds T = promoteunion(T, typeof(toarrow(x[i])))
        end
    end
    return ToArrow{T, A}(x)
end

Base.IndexStyle(::Type{<:ToArrow}) = Base.IndexLinear()
Base.size(x::ToArrow) = (length(x.data),)
Base.eltype(x::ToArrow{T, A}) where {T, A} = T
Base.getindex(x::ToArrow{T}, i::Int) where {T} = toarrow(getindex(x.data, i))

#### DEPRECATED
function arrowconvert end
arrowconvert(T, x) = convert(T, x)
arrowconvert(::Type{Union{T, Missing}}, x) where {T} = arrowconvert(T, x)
arrowconvert(::Type{Union{T, Missing}}, ::Missing) where {T} = missing

const JULIA_TO_ARROW_TYPE_MAPPING = Dict{Type, Tuple{String, Type}}()

istyperegistered(::Type{T}) where {T} = haskey(JULIA_TO_ARROW_TYPE_MAPPING, T)

function getarrowtype!(meta, ::Type{T}) where {T}
    arrowname, arrowtype = JULIA_TO_ARROW_TYPE_MAPPING[T]
    meta["ARROW:extension:name"] = arrowname
    meta["ARROW:extension:metadata"] = ""
    return arrowtype
end

const ARROW_TO_JULIA_TYPE_MAPPING = Dict{String, Tuple{Type, Type}}()

function extensiontype(f, meta)
    if haskey(meta, "ARROW:extension:name")
        typename = meta["ARROW:extension:name"]
        if haskey(ARROW_TO_JULIA_TYPE_MAPPING, typename)
            T = ARROW_TO_JULIA_TYPE_MAPPING[typename][1]
            return f.nullable ? Union{T, Missing} : T
        end
    end
    return nothing
end

function registertype!(juliatype::Type, arrowtype::Type, arrowname::String=string("JuliaLang.", string(juliatype)))
    @warn "Arrow.ArrowTypes.registertype! is deprecated in favor of defining `ArrowTypes.ArrowType`, `ArrowTypes.arrowname`, and `ArrowTypes.JuliaType`; see their docs for more information on how to migrate"
    # TODO: validate that juliatype isn't already default arrow type
    JULIA_TO_ARROW_TYPE_MAPPING[juliatype] = (arrowname, arrowtype)
    ARROW_TO_JULIA_TYPE_MAPPING[arrowname] = (juliatype, arrowtype)
    return
end

end # module ArrowTypes


text:
export FPNLearner, FullyParametrizedNet

using Flux
using CUDA
using Random
using Zygote
using Statistics: mean
using LinearAlgebra: dot

"""
    ImplicitQuantileNet(;ψ, ϕ, header)

```
        quantiles (n_action, n_quantiles, batch_size)
           ↑
         header
           ↑
feature ↱  ⨀   ↰ transformed embedding
       ψ       ϕ
       ↑       ↑
       s        τ
```
"""
Base.@kwdef struct FullyParametrizedNet{A,B,C} <: AbstractApproximator
    ψ::A
    ϕ::B
    header::C
end

Flux.@functor FullyParametrizedNet

function (net::FullyParametrizedNet)(s, emb)
    #fpn = net.fpn_app(s)
    features = net.ψ(s)  # (n_feature, batch_size)
    emb_aligned = net.ϕ(emb)  # (n_feature, N * batch_size)
    merged =
        Flux.unsqueeze(features, 2) .*
        reshape(emb_aligned, size(features, 1), :, size(features, 2))  # (n_feature, N, batch_size)
    quantiles = net.header(flatten_batch(merged))
    reshape(quantiles, :, size(merged, 2), size(merged, 3))  # (n_action, N, batch_size)
end

"""
    FPNLearner(;kwargs)

See [paper](https://arxiv.org/abs/1806.06923)

# Keyworkd arugments
- `approximator`, a [`FullyParametrizedNet`](@ref)
- `target_approximator`, a [`FullyParametrizedNet`](@ref), must have the same structure as `approximator`
- `κ = 1.0f0`,
- `N = 32`,
- `N′ = 32`,
- `Nₑₘ = 64`,
- `K = 32`,
- `γ = 0.99f0`,
- `stack_size = 4`,
- `batch_size = 32`,
- `update_horizon = 1`,
- `min_replay_history = 20000`,
- `update_freq = 4`,
- `target_update_freq = 8000`,
- `update_step = 0`,
- `default_priority = 1.0f2`,
- `β_priority = 0.5f0`,
- `rng = Random.GLOBAL_RNG`,
- `device_seed = nothing`,
"""
Base.@kwdef mutable struct FPNLearner{A,T,B<:NeuralNetworkApproximator,R,D} <: AbstractLearner
    approximator::A
    target_approximator::T
    fpn_app::B
    sampler::NStepBatchSampler
    κ::Float32
    N::Int
    N′::Int
    Nₑₘ::Int
    K::Int
    γ::Float32
    stack_size::Union{Nothing,Int}
    batch_size::Int
    update_horizon::Int
    min_replay_history::Int
    update_freq::Int
    target_update_freq::Int
    update_step::Int
    default_priority::Float32
    β_priority::Float32
    rng::R
    device_rng::D
    loss::Float32
end

Flux.functor(x::FPNLearner) =
    (Z = x.approximator, Zₜ = x.target_approximator,f = x.fpn_app, device_rng = x.device_rng),
    y -> begin
        x = @set x.approximator = y.Z
        x = @set x.target_approximator = y.Zₜ
        x = @set x.fpn_app = y.f
        x = @set x.device_rng = y.device_rng
        x
    end

Flux.gpu(rng::MersenneTwister) = CUDA.CURAND.RNG()
Flux.cpu(rng::CUDA.CURAND.RNG) = MersenneTwister()

function FPNLearner(;
    approximator,
    target_approximator,
    fpn_app,
    κ = 1.0f0,
    N = 32,
    N′ = 32,
    Nₑₘ = 64,
    K = 32,
    γ = 0.99f0,
    stack_size = 4,
    batch_size = 32,
    update_horizon = 1,
    min_replay_history = 20000,
    update_freq = 4,
    target_update_freq = 8000,
    update_step = 0,
    default_priority = 1.0f2,
    β_priority = 0.5f0,
    rng = Random.GLOBAL_RNG,
    device_rng = CUDA.CURAND.RNG(),
    traces = SARTS,
    loss = 0.0f0,
)
    copyto!(approximator, target_approximator )  # force sync
    if device(approximator) !== device(device_rng)
        throw(ArgumentError("device of `approximator` doesn't match with the device of `device_rng`: $(device(approximator)) !== $(device_rng)"))
    end
    sampler = NStepBatchSampler{traces}(;γ=γ, n=update_horizon,stack_size=stack_size,batch_size=batch_size)
    FPNLearner(
        approximator,
        target_approximator,
        fpn_app,
        sampler,
        κ,
        N,
        N′,
        Nₑₘ,
        K,
        γ,
        stack_size,
        batch_size,
        update_horizon,
        min_replay_history,
        update_freq,
        target_update_freq,
        update_step,
        default_priority,
        β_priority,
        rng,
        device_rng,
        loss,
    )
end

function (learner::FPNLearner)(env)
    state = send_to_device(device(learner.approximator), get_state(env))
    state = Flux.unsqueeze(state, ndims(state) + 1)
    tau_i = cumsum(softmax(learner.fpn_app(state); dims=1),dims=1)
    ze = zeros(Float32,1,size(tau_i,2))
    tau = cat(ze,tau_i,dims=1)
    N = learner.N
    tau_hat = ones(Float32,N,size(tau,2))
    for i=1:N
        tau_hat[i,:]=(tau[i,:].+tau[i+1,:])./2
    end
    τ = tau #rand(learner.device_rng, Float32, learner.K, 1)
    τₑₘ = embed(τ, learner.Nₑₘ)
    quantiles = learner.approximator(state, τₑₘ)
    vec(mean(quantiles; dims = 2)) |> send_to_host
end

embed(x, Nₑₘ) = cos.(Float32(π) .* (1:Nₑₘ) .* reshape(x, 1, :))

function RLBase.update!(learner::FPNLearner, t::AbstractTrajectory)
    length(t[:terminal]) < learner.min_replay_history && return

    learner.update_step += 1

    if learner.update_step % learner.target_update_freq == 0
        copyto!(learner.target_approximator, learner.approximator)
    end

    learner.update_step % learner.update_freq == 0 || return

    inds, batch = sample(learner.rng, t, learner.sampler)

    if t isa PrioritizedTrajectory
        priorities = update!(learner, batch)
        t[:priority][inds] .= priorities
    else
        update!(learner,batch)
    end
end

function RLBase.update!(learner::FPNLearner, batch::NamedTuple)
    Z = learner.approximator
    Zₜ = learner.target_approximator
    f = learner.fpn_app
    N = learner.N
    N′ = learner.N′
    Nₑₘ = learner.Nₑₘ
    κ = learner.κ
    β = learner.β_priority
    batch_size = learner.batch_size

    #D = device(Z)
    D = device(Z)
    s, r, t, s′ = (send_to_device(D, batch[x]) for x in (:state, :reward, :terminal, :next_state))

    tau_i_p = cumsum(softmax(f(s′); dims=1),dims=1)
    ze_p = zeros(Float32,1,batch_size)
    tau_p = cat(ze_p,tau_i_p,dims=1)
    tau_hat_p = ones(Float32,N,size(tau_p,2))
    for i=1:N
        tau_hat_p[i,:]=(tau_p[i,:].+tau_p[i+1,:])./2
    end
    τ′ = tau_hat_p[1:N,:] #rand(learner.device_rng, Float32, N′, batch_size)  # TODO: support β distribution
    τₑₘ′ = embed(τ′, Nₑₘ)
    zₜ = Zₜ(s′, τₑₘ′)
    #z = reshape(zₜ,size(zₜ,1),N+1,batch_size)
    #avg_zₜ = zeros(Float32,size(zₜ,1),1,batch_size)
    p = zeros(Float32,1,batch_size)
    #tau_pf = tau_p[1:N,:]
    #for j=1:size(zₜ,1)
    #for k=1:size(avg_zₜ,1)
        #for j=1:batch_size
            #c = 0
            #for i=1:N
                #c = c .+ (tau_p[i+1,j]-tau_p[i,j]).*zₜ[k,i,j]
            #end
            #avg_zₜ[k,1,j] = c
        #end
    #end
    for i=1:N
        zₜ[:,i,:] = reshape((tau_p[i+1,:].-tau_p[i,:]),1,1,:) .*zₜ[:,i,:]
    end
    avg_zₜ = sum(zₜ,dims=2)

    aₜ = argmax(avg_zₜ, dims = 1)
    aₜ = aₜ .+ typeof(aₜ)(CartesianIndices((0, 0:N-1, 0)))
    qₜ = reshape(zₜ[aₜ], :, batch_size)
    target = reshape(r, 1, batch_size) .+ learner.sampler.γ * reshape(1 .- t, 1, batch_size) .* qₜ

        #mean(zₜ, dims = 2)

    if haskey(batch, :next_legal_actions_mask)
        masked_value = fill(typemin(Float32), size(batch.next_legal_actions_mask))
        masked_value[batch.next_legal_actions_mask] .= 0
        avg_zₜ .+= send_to_device(D, masked_value)
    end
     # reshape to allow broadcast

    tau_i = cumsum(softmax(f(s); dims=1),dims=1)
    ze = zeros(Float32,1,batch_size)
    tau = cat(ze,tau_i,dims=1)
    tau_hat = ones(Float32,N,size(tau,2))
    for i=1:N
        tau_hat[i,:]=(tau[i,:].+tau[i+1,:])./2
    end
    τ = tau_hat[1:N,:]#rand(learner.device_rng, Float32, N, batch_size)
    τₑₘ = embed(τ, Nₑₘ)
    a = CartesianIndex.(repeat(batch.actions, inner = N), 1:(N*batch_size))

    tau_em = embed(tau[1:N,:],Nₑₘ)
    #q_t = Z(s, tau_em)
    #q_t = reshape(q_t_i,:,batch_size)
    #z = reshape(zₜ,size(zₜ,1),N+1,batch_size)

    is_use_PER = haskey(batch.priorities)  # is use Prioritized Experience Replay
    if is_use_PER
        updated_priorities = Vector{Float32}(undef, batch_size)
        weights = 1.0f0 ./ ((batch.priorities .+ 1f-10) .^ β)
        weights ./= maximum(weights)
        weights = send_to_device(D, weights)
    end
    gs1 = Zygote.gradient(Flux.params(f)) do
        z1 = flatten_batch(Z(s, tau_em))
        z2 = flatten_batch(Z(s, τₑₘ))
        dW1 = zeros(Float32,N-1,batch_size)
        for i=2:N
            dW1[i,:] = 2 .* reshape(z1[a],N,batch_size)[i,:] .- reshape(z2[a],N,batch_size)[i,:] .- reshape(z2[a],N,batch_size)[i-1,:]
        end
        loss = Zygote.dropgrad(mean(sum(dW1,dims=1),dims=2)[1,1]) .* tau_em
    end
    update!(f, gs1)

    gs = Zygote.gradient(Flux.params(Z)) do
        z = flatten_batch(Z(s, τₑₘ))
        q = z[a]

        TD_error = reshape(target, N, 1, batch_size) .- reshape(q, 1, N, batch_size)
        # can't apply huber_loss in RLCore directly here
        abs_error = abs.(TD_error)
        quadratic = min.(abs_error, κ)
        linear = abs_error .- quadratic
        huber_loss = 0.5f0 .* quadratic .* quadratic .+ κ .* linear

        # dropgrad
        raw_loss =
            abs.(reshape(τ, 1, N, batch_size) .- Zygote.dropgrad(TD_error .< 0)) .*
            huber_loss ./ κ
        loss_per_quantile = reshape(sum(raw_loss; dims = 1), N, batch_size)
        loss_per_element = mean(loss_per_quantile; dims = 1)  # use as priorities
        loss =
            is_use_PER ? dot(vec(weights), vec(loss_per_element)) * 1 // batch_size :
            mean(loss_per_element)
        ignore() do
            # @assert all(loss_per_element .>= 0)
            is_use_PER && (
                updated_priorities .=
                    send_to_host(vec((loss_per_element .+ 1f-10) .^ β))
            )
            learner.loss = loss
        end
        loss
    end

    update!(Z, gs)

    is_use_PER ? updated_priorities : nothing
end


text:
if !@isdefined(CGS)
  include("../../../src/CGS.jl")
  using Main.CGS
end

using JLD2; using PyPlot
using DelimitedFiles
using Statistics

#rhostar = 51.14; sig_rho=0.67 # In units of solar density
#rhostar = 53.7; sig_rho=0.85 # In units of solar density
#rhostar = 53.308; sig_rho=0.0208 # In units of solar density
#rhostar = 51.14; sig_rho=0.001 # In units of solar density

# Artificially increase density of star to see what it
# would take to make an Earth-like density:
#rhostar *= 1.23

# Mass of star from Mann et al. (2019):  https://github.com/awmann/M_-M_K-
mstar = 0.0898; smstar = 0.0023  # Solar units
#mstar = 0.0898; smstar = 0.0001  # Solar units

# Latest masses of planets assuming a stellar mass of 0.09 M_sun.
# Take from plot_masses.jl in
# /Users/ericagol/Observing/Spitzer/Cycle13/TRAPPIST1/campaign03/April2019
# Includes delta-omega prior, includes mixture model for handling outliers.

#mp= [1.37024,1.31417,0.364033,0.652497,0.997529,1.28074,0.340946]
#sp= [0.0848279,0.0602858,0.022669,0.0230492,0.0202153,0.022538,0.0355592]
#mp = [1.3785, 1.3201, 0.3911, 0.7042, 1.0463, 1.3237, 0.3283]
#sp = [0.0530, 0.0397, 0.0076, 0.0116, 0.0128, 0.0146, 0.0190]

#mp = [1.37927,1.31014,0.388479,0.693422,1.04092,1.32341,0.325985]
#sp = [0.0596173,0.0465207,0.00748253,0.0127809,0.0154444,0.0168958,0.0185088]

# Load in the posterior from the transit-timing analysis:
#@load "/Users/ericagol/Observing/Spitzer/DDT2019/Campaign04/v09/Hyak/T1_hmc_total_02132020.jld2" state_total
@load "../../../data/T1_hmc_total_02212020.jld2" state_total

# Load in the posterior from the photodynamic analysis:
@load "../../../data/T1_pd_MCMC_run_000_001.jld2" flatchain
flatchain1 = flatchain[:,35001:200000]
#@load "/Users/ericagol/Observing/Spitzer/DDT2019/Campaign04/v06/Hyak/T1_pd_MCMC_run_000_002.jld2" flatchain
@load "../../../data/T1_pd_MCMC_run_000_002.jld2" flatchain
flatchain12 = hcat(flatchain1,flatchain[:,35001:200000])
#@load "/Users/ericagol/Observing/Spitzer/DDT2019/Campaign04/v06/Hyak/T1_pd_MCMC_run_000_003.jld2" flatchain
@load "../../../data/T1_pd_MCMC_run_000_003.jld2" flatchain
flatchain = hcat(flatchain12,flatchain[:,35001:200000])

# Make an array to hold posterior of mass-ratios:
nplanet = 7
mass_ratio_total = zeros(nplanet,size(state_total)[2])
for ip=1:nplanet
  mass_ratio_total[ip,:] .= state_total[(ip-1)*5+1,:]
end

periods = [1.51087637,2.42180746,4.049959,6.099043,9.205585,12.354473,18.767953]
t0 = [7322.51654,7282.80879,7670.14227,7660.37910,7671.39470,7665.35084,7662.55467]

logplot = false
include("mass_radius2.jl")

#merror = mass_radius2(mass_ratio_total,mstar,smstar,flatchain,false)

# Now, plot mass-radius relations from
# Zeng https://www.cfa.harvard.edu/~lzeng/tables/mrtable3.txt

dmr = readdlm("../../../data/zeng_mrtable3.txt",skipstart=2)

#plot(dmr[:,1],dmr[:,16],label="30% Fe",linestyle="--")  # 30% Fe
#plot(dmr[:,1],dmr[:,20],label="10% Fe",linestyle="--")  # 10% Fe
#plot(dmr[:,1],dmr[:,23],label=" 0% Fe",linestyle="--")  #  0% Fe
scatter([1.0,0.815],[1.0,0.95])

# Plot new Zeng rocky model:
# (from https://www.cfa.harvard.edu/~lzeng/planetmodels.html#mrtables)
dmr1 = readdlm("../../../data/zeng_rock.txt")
#plot(dmr1[1:18,1],dmr1[1:18,2],label="Rocky (Zeng et al. 2019)",linestyle="--")

# Plot new Zeng Earth-like mode:
dmr2 = readdlm("../../../data/massradiusEarthlikeRocky.txt")
plot(dmr2[1:23,1],dmr2[1:23,2],label="Earth-like (Zeng et al. 2019)",linestyle="--")

# Read in Carbon-rich planet model (for kicks):

#dmr_carbon = readdlm("../../../data/M_R_carbon-rich_Eric.ddat")
#plot(dmr_carbon[:,1],dmr_carbon[:,2],label="Carbon-rich (Miozzi et al. 2018)")

# Interpolate between the two Zeng models:
nm = 100
mm = logarithmspace(-2.5,0.2,nm)

function interpolate_models(iron,dmr1,dmr2,mm,nm)
  # Loop over each model and interpolate/extrapolate:
  r1 = zeros(nm)
  r2 = zeros(nm)
  r = zeros(nm)
  for i=1:nm
    if mm[i] < dmr1[1,1]
      i1 = 1; i2 = 2
    else
      for j=1:size(dmr1)[1]
        if mm[i] > dmr1[j,1] && mm[i] <= dmr1[j+1,1]
           i1 = j; i2 = j+1
        end
      end
    end
    r1[i] = dmr1[i1,2]*(dmr1[i2,2]/dmr1[i1,2])^((log(mm[i]/dmr1[i1,1]))/(log(dmr1[i2,1]/dmr1[i1,1])))
    if mm[i] < dmr2[1,1]
      i1 = 1; i2 = 2
    else
      for j=1:size(dmr2)[1]
        if mm[i] > dmr2[j,1] && mm[i] <= dmr2[j+1,1]
           i1 = j; i2 = j+1
        end
      end
    end
    r2[i] = dmr2[i1,2]*(dmr2[i2,2]/dmr2[i1,2])^(log(mm[i]/dmr2[i1,1])/log(dmr2[i2,1]/dmr2[i1,1]))
    # Now interpolate to the right iron content:
    r[i] = r1[i]+(r2[i]-r1[i])*iron/0.325
  end
  return r
end

#plot(mm,interpolate_models(0.10,dmr1,dmr2,mm,nm),label="10% iron (Zeng et al. 2019)",linestyle="--")
#plot(mm,interpolate_models(0.25,dmr1,dmr2,mm,nm),label="25% iron (Zeng et al. 2019)",linestyle="--")

legend(loc = "lower right",fontsize=10)

tight_layout()
#savefig("T1_mass_radius_Earth_carbon.pdf",bbox_inches="tight")
#read(stdin,Char)

## Semi-amplitudes:
fac = (2*pi*GRAV/(24*3600)/MSUN^2)^(1//3)*MEARTH
Kp = fac.*mavg./(mstar.^2 .*periods).^(1//3)
# Plot RV amplitude:
duration = 300
nrv = duration * 24
RV = zeros(nplanet+1,nrv)
time = zeros(nrv)
t0 .-= 7600.0
for i=1:nrv
  time[i] = i/24.0
  for j=1:nplanet
    RV[j,i] = Kp[j].*sin.((time[i].-t0[j])./periods[j])
  end
  RV[nplanet+1,i] = sum(RV[1:nplanet,i])
#    RV[j,i] = Kp[j].*sin.((time[i].-t0[j])./periods[j]))
end

clf()
plot(time,RV[nplanet+1,:]./100,label="Sum",color="k")
#plot(time,RV[nplanet+1,:],label="Sum",color="k")
xlabel("Time [days]")
#ylabel("RV [cm/sec]")
ylabel("RV [m/sec]")
#axis([0,duration,-1200,1200])
axis([0,duration,-12.00,12.00])
for i=1:7
#  plot(time,RV[i,:],label=planet_name[i],color=cname[i])
  plot(time,RV[i,:]./100,label=planet_name[i],color=cname[i])
#  errorbar([120+10*i],[-1000],Kp[i]*merror[i]/mavg[i],color=cname[i])
  errorbar([90+10*i],[-10],Kp[i]*merror[i]/mavg[i]/100,color=cname[i])
end
#text(40,-1020,"Equivalent RV precision:",fontsize=10)
text(10,-10.20,"Equivalent RV precision:",fontsize=10)
text(180,-10.2,"Measured RV precision: ")
errorbar([280],[-9.0],2.5,label="Hirano et al. (2020)",fmt="o")
legend(ncol=3)
tight_layout()
savefig("../Equivalent_RV_precision.pdf",bbox_inches="tight")

println("Semi-amplitudes: ",Kp," cm/s")
println("Precision: ",merror./mavg.*Kp," cm/s")
Kpsun = fac/(365.25)^(1//3)
println("Earth-Sun amplitude: ",Kpsun)
println("Equivalent Solar precision: ",Kpsun.*merror)

#include("/Users/ericagol/JWST/Exomoon/plot_moons.jl")

#@save "T1_mass_radius_cmf_000.jld2" nsamp msamp rsamp cmf rho_samp mstar_samp rstar_samp nr nm mrgrid rgrid mgrid level_mr conf m1 m2 r1 r2 mavg ravg merror rerror


text:
"""
## geom_rect!

This subroutine forms the coordinates and connectivity for a
rectangular mesh of right angled triangular elements (3, 6, 10 or 15-node)
or quadrilateral elements (4, 8 or 9-node) counting in the
x- or y-dir.

### Method
```julia
geom_rect!(fin_el, iel, x_coords, y_coords, coord, num, dir)
```

### Arguments
```julia
* fin_el::FiniteElement            : Shape of finite element
                                     (Trangle or Quadrilateral)
* iel::Int                       : Element number
* x_coords::FloatRange{Float64}    : x coordinates
* y_coords::FloatRange{Float64}    : y coordinates
* coord::Matrix{Float64}           : Nodal coordinates (Updated)
* num::Vector{Int}               : Node numbers (Updated)
* dir::Symbol                      : Node numbering direction
```
"""
function geom_rect!(fin_el::Triangle, iel::Int, x_coords::Array{Float64, 1},
  y_coords::Array{Float64, 1}, coord::Matrix{Float64}, num::Vector{Int}, dir::Symbol)

  nxe = size(x_coords, 1) - 1
  nye = (size(y_coords, 1) - 1) * 2
  nod = size(num, 1)
  #println([nxe nye nod])
  if dir == :x || dir == :r
    jel = round.(Int, 2 * nxe * floor((iel - 1) ÷ (2 * nxe)))
    ip = round.(Int, floor((iel - jel + 1) ÷ 2))
    iq = round.(Int, floor(2 * floor(floor((iel - 1) ÷ (2 * nxe)) + 1) - 1 + floor(floor(iel ÷ 2) * 2) / iel))
  else
    #println("Direction != :x or :r.")
    jel = round.(Int, floor((iel - 1) ÷ nye))
    ip = jel + 1
    iq = iel - nye * jel
  end
  #println([jel ip iq])
  if nod == 3
    if mod(iq, 2) != 0
      if dir == :x || dir == :r
        num[1] = (nxe + 1) * round.(Int, (iq - 1) / 2) + ip
        num[2] = num[1] + 1
        num[3] = (nxe + 1) * round.(Int, (iq + 1) / 2) + ip
      else
        num[1] = (ip - 1) * round.(Int, (nye + 2) / 2) + (iq + 1) / 2
        num[2] = num[1] + round.(Int, (nye + 2) / 2)
        num[3] = num[1] + 1
      end
      coord[1, 1] = x_coords[ip]
      coord[1, 2] = y_coords[round.(Int, (iq + 1) / 2)]
      coord[2, 1] = x_coords[ip + 1]
      coord[2, 2] = y_coords[round.(Int, (iq + 1) / 2)]
      coord[3, 1] = x_coords[ip]
      coord[3, 2] = y_coords[round.(Int, (iq + 3) / 2)]
    else
      if dir == :x || dir == :r
        num[1] = (nxe + 1) * round.(Int, iq / 2) + ip + 1
        num[2] = num[1] - 1
        num[3] = (nxe + 1) * round.(Int, (iq - 2) / 2) + ip + 1
      else
        num[1] = ip * round.(Int, (nye + 2) / 2) + round.(Int, (iq + 2) / 2)
        num[2] = (ip - 1) * round.(Int, (nye + 2) / 2) + round.(Int, (iq + 1) / 2 + 1)
        num[3] = num[1] - 1
      end
      coord[1, 1] = x_coords[ip+ 1]
      coord[1, 2] = y_coords[round.(Int, (iq + 2) / 2)]
      coord[2, 1] = x_coords[ip]
      coord[2, 2] = y_coords[round.(Int, (iq + 2) / 2)]
      coord[3, 1] = x_coords[ip + 1]
      coord[3, 2] = y_coords[round.(Int, iq / 2)]
    end
  elseif nod == 6
    if mod(iq, 2) !== 0
     if dir == 'x' || dir == 'r'
       num[1] = (iq-1)*(2*nxe+1)+2*ip-1
       num[2] = num[1]+1
       num[3] = num[1]+2
       num[4] = (iq-1)*(2*nxe+1)+2*nxe+2*ip+1
       num[5] = (iq+1)*(2*nxe+1)+2*ip-1
       num[6] = num[4]-1
     else
       num[1] = 2*(nye+1)*(ip-1)+iq
       num[2] = 2*(nye+1)*(ip-1)+nye+1+iq
       num[3] = 2*(nye+1)*ip+iq
       num[4] = num[2]+1
       num[5] = num[1]+2
       num[6] = num[1]+1
     end
     coord[1,1] = x_coords[ip]
     coord[1,2] = y_coords[(iq+1)÷2]
     coord[3,1] = x_coords[ip+1]
     coord[3,2] = y_coords[(iq+1)÷2]
     coord[5,1] = x_coords[ip]
     coord[5,2] = y_coords[(iq+3)÷2]
    else
      if dir == 'x' || dir == 'r'
       num[1] = iq*(2*nxe+1)+2*ip+1
       num[2] = num[1]-1
       num[3] = num[1]-2
       num[4] = (iq-2)*(2*nxe+1)+2*nxe+2*ip+1
       num[5] = (iq-2)*(2*nxe+1)+2*ip+1
       num[6] = num[4]+1
      else
       num[1] = 2*(nye+1)*ip+iq+1
       num[2] = 2*(nye+1)*(ip-1)+nye+iq+2
       num[3] = 2*(nye+1)*(ip-1)+iq+1
       num[4] = num[2]-1
       num[5] = num[1]-2
       num[6] = num[1]-1
      end
      coord[1,1] = x_coords[ip+1]
      coord[1,2] = y_coords[(iq+2)÷2]
      coord[3,1] = x_coords[ip]
      coord[3,2] = y_coords[(iq+2)÷2]
      coord[5,1] = x_coords[ip+1]
      coord[5,2] = y_coords[iq÷2]
    end
    coord[2,:] = 0.5*(coord[1,:]+coord[3,:])
    coord[4,:] = 0.5*(coord[3,:]+coord[5,:])
    coord[6,:] = 0.5*(coord[5,:]+coord[1,:])
  elseif nod == 10
    if mod(iq, 2) !== 0
     if dir == 'x' || dir == 'r'
       num[1] = (iq-1)/2*(3*nxe+1)*3+3*ip-2
       num[2] = num[1]+1
       num[3] = num[1]+2
       num[4] = num[1]+3
       num[5] = (iq-1)/2*(3*nxe+1)*3+3*nxe+1+3*ip
       num[6] = (iq-1)/2*(3*nxe+1)*3+6*nxe+2+3*ip-1
       num[7] = (iq-1)/2*(3*nxe+1)*3+9*nxe+3+3*ip-2
       num[8] = num[6]-1
       num[9] = num[5]-2
       num[10] = num[9]+1
     else
       num[1] = (9*(nye-2)/2+12)*(ip-1)+3*(iq-1)/2+1
       num[2] = (9*(nye-2)/2+12)*(ip-1)+3*(nye-2)/2+4+3*(iq-1)/2+1
       num[3] = (9*(nye-2)/2+12)*(ip-1)+3*(nye-2)+8+3*(iq-1)/2+1
       num[4] = (9*(nye-2)/2+12)*(ip-1)+9*(nye-2)/2+12+3*(iq-1)/2+1
       num[5] = num[3]+1
       num[6] = num[2]+2
       num[7] = num[1]+3
       num[8] = num[1]+2
       num[9] = num[1]+1
       num[10] = num[2]+1
     end
     coord[1,1] = x_coords[ip]
     coord[2,1] = x_coords[ip]+(x_coords[ip+1]-x_coords[ip])/3.0
     coord[3,1] = x_coords[ip]+2.0*(x_coords[ip+1]-x_coords[ip])/3.0
     coord[4,1] = x_coords[ip+1]
     coord[4,2] = y_coords[(iq+1)/2]
     coord[5,2] = y_coords[(iq+1)/2]+(y_coords[(iq+3)/2]-y_coords[(iq+1)/2])/3.0
     coord[6,2] = y_coords[(iq+1)/2]+2.0*(y_coords[(iq+3)/2]-y_coords[(iq+1)/2])/3.0
     coord[7,2] = y_coords[(iq+3)/2]
    else
     if dir == 'x' || dir == 'r'
       num[1] = (iq-2)/2*(3*nxe+1)*3+9*nxe+3+3*ip+1
       num[2] = num[1]-1
       num[3] = num[1]-2
       num[4] = num[1]-3
       num[5] = (iq-2)/2*(3*nxe+1)*3+6*nxe+2+3*ip-1
       num[6] = (iq-2)/2*(3*nxe+1)*3+3*nxe+1+3*ip
       num[7] = (iq-2)/2*(3*nxe+1)*3+3*ip+1
       num[8] = num[6]+1
       num[9] = num[5]+2
       num[10] = num[9]-1
     else
       num[1] = (9*(nye-2)/2+12)*(ip-1)+9*(nye-2)/2+12+3*iq/2+1
       num[2] = (9*(nye-2)/2+12)*(ip-1)+3*(nye-2)+8+3*iq/2+1
       num[3] = (9*(nye-2)/2+12)*(ip-1)+3*(nye-2)/2+4+3*iq/2+1
       num[4] = (9*(nye-2)/2+12)*(ip-1)+3*iq/2+1
       num[5] = num[3]-1
       num[6] = num[2]-2
       num[7] = num[1]-3
       num[8] = num[1]-2
       num[9] = num[1]-1
       num[10] = num[2]-1
     end
     coord[1,1] = x_coords[ip+1]
     coord[2,1] = x_coords[ip+1]-(x_coords[ip+1]-x_coords[ip])/3.0
     coord[3,1] = x_coords[ip+1]-2.0*(x_coords[ip+1]-x_coords[ip])/3.0
     coord[4,1] = x_coords[ip]
     coord[4,2] = y_coords[(iq+2)/2]
     coord[5,2] = y_coords[(iq+2)/2]-(y_coords[(iq+2)/2]-y_coords[iq/2])/3.0
     coord[6,2] = y_coords[(iq+2)/2]-2.0*(y_coords[(iq+2)/2]-y_coords[iq/2])/3.0
     coord[7,2] =y_coords[iq/2]
    end
    coord[5,1] = coord[3,1]
    coord[6,1] = coord[2,1]
    coord[7,1] = coord[1,1]
    coord[8,1] = coord[1,1]
    coord[9,1] = coord[1,1]
    coord[10,1] = coord[2,1]
    coord[1,2] = coord[4,2]
    coord[2,2] = coord[4,2]
    coord[3,2] = coord[4,2]
    coord[8,2] = coord[6,2]
    coord[9,2] = coord[5,2]
    coord[10,2] = coord[5,2]
  elseif nod == 15
   if mod(iq, 2) !== 0
     if dir == 'x' || dir == 'r'
       fac1 = 4*(4*nxe+1)*(iq-1)/2
       num[1] = fac1+4*ip-3
       num[2] = num[1]+1
       num[3] = num[1]+2
       num[4] = num[1]+3
       num[5] = num[1]+4
       num[6] = fac1+ 4*nxe+1+4*ip
       num[7] = fac1+ 8*nxe+1+4*ip
       num[8] = fac1+12*nxe+1+4*ip
       num[9] = fac1+16*nxe+1+4*ip
       num[10] = num[8]-1
       num[11] = num[7]-2
       num[12] = num[6]-3
       num[13] = num[12]+1
       num[14] = num[12]+2
       num[15] = num[11]+1
     else
       fac1=4*(2*nye+1)*(ip-1)+2*iq-1
       num[1] = fac1
       num[2] = fac1+2*nye+1
       num[3] = fac1+4*nye+2
       num[4] = fac1+6*nye+3
       num[5] = fac1+8*nye+4
       num[6] = fac1+6*nye+4
       num[7] = fac1+4*nye+4
       num[8] = fac1+2*nye+4
       num[9] = fac1+4
       num[10] = fac1+3
       num[11] = fac1+2
       num[12] = fac1+1
       num[13] = fac1+2*nye+2
       num[14] = fac1+4*nye+3
       num[15] = fac1+2*nye+3
     end
     coord[1,1] = x_coords[ip]
     coord[1,2] = y_coords[round.(Int, (iq+1)/2)]
     coord[5,1] = x_coords[ip+1]
     coord[5,2] = y_coords[round.(Int, (iq+1)/2)]
     coord[9,1] = x_coords[ip]
     coord[9,2] = y_coords[round.(Int, (iq+3)/2)]
   else
     if dir == 'x' || dir == 'r'
       fac1=4*(4*nxe+1)*(iq-2)/2
       num[1] = fac1+16*nxe+5+4*ip
       num[2] = num[1]-1
       num[3] = num[1]-2
       num[4] = num[1]-3
       num[5] = num[1]-4
       num[6] = fac1+12*nxe+1+4*ip
       num[7] = fac1+8*nxe+1+4*ip
       num[8] = fac1+4*nxe+1+4*ip
       num[9] = fac1+4*ip+1
       num[10] = num[8]+1
       num[11] = num[7]+2
       num[12] = num[6]+3
       num[13] = num[12]-1
       num[14] = num[12]-2
       num[15] = num[11]-1
     else
       fac1=4*(2*nye+1)*(ip-1)+2*iq+8*nye+5
       num[1] = fac1
       num[2] = fac1-2*nye-1
       num[3] = fac1-4*nye-2
       num[4] = fac1-6*nye-3
       num[5] = fac1-8*nye-4
       num[6] = fac1-6*nye-4
       num[7] = fac1-4*nye-4
       num[8] = fac1-2*nye-4
       num[9] = fac1-4
       num[10] = fac1-3
       num[11] = fac1-2
       num[12] = fac1-1
       num[13] = fac1-2*nye-2
       num[14] = fac1-4*nye-3
       num[15] = fac1-2*nye-3
     end
     coord[1,1] = x_coords[ip+1]
     coord[1,2] = y_coords[round.(Int, (iq+2)/2)]
     coord[5,1] = x_coords[ip]
     coord[5,2] = y_coords[round.(Int, (iq+2)/2)]
     coord[9,1] = x_coords[ip+1]
     coord[9,2] = y_coords[round.(Int, iq/2)]
   end
   coord[3,:] = 0.5*(coord[1,:]+coord[5,:])
   coord[7,:] = 0.5*(coord[5,:]+coord[9,:])
   coord[11,:] = 0.5*(coord[9,:]+coord[1,:])
   coord[2,:] = 0.5*(coord[1,:]+coord[3,:])
   coord[4,:] = 0.5*(coord[3,:]+coord[5,:])
   coord[6,:] = 0.5*(coord[5,:]+coord[7,:])
   coord[8,:] = 0.5*(coord[7,:]+coord[9,:])
   coord[10,:] = 0.5*(coord[9,:]+coord[11,:])
   coord[12,:] = 0.5*(coord[11,:]+coord[1,:])
   coord[15,:] = 0.5*(coord[7,:]+coord[11,:])
   coord[14,:] = 0.5*(coord[3,:]+coord[7,:])
   coord[13,:] = 0.5*(coord[2,:]+coord[15,:])
  else
   println("Wrong number of nodes for triangular fin_el.")
   exit(1)
  end
end

function geom_rect!(fin_el::Quadrilateral, iel::Int, x_coords::Array{Float64, 1},
  y_coords::Array{Float64, 1}, coord::Matrix{Float64}, num::Vector{Int}, dir::Symbol)

  nxe = size(x_coords,1)-1
  nye = size(y_coords,1)-1
  nod = size(num, 1)
  if dir == :x || dir == :r
    iq  = round.(Int, floor((iel-1)/nxe))+1
    ip = iel-(iq-1)*nxe
  else
    ip = round.(Int, floor((iel-1)/nye))+1
    iq = iel-(ip-1)*nye
  end
  if nod == 4
    if dir == :x || dir == :r
     num[1] = iq*(nxe+1)+ip
     num[2] = (iq-1)*(nxe+1)+ip
     num[3] = num[2]+1
     num[4] = num[1]+1
    else
     num[1] = (ip-1)*(nye+1)+iq+1
     num[2] = num[1]-1
     num[3] = ip*(nye+1)+iq
     num[4] = num[3]+1
    end
    coord[1:2,1] .= x_coords[ip]
    coord[3:4,1] .= x_coords[ip+1]
    coord[1,2] = y_coords[iq+1]
    coord[2:3,2] .= y_coords[iq]
    coord[4,2] = coord[1,2]
  elseif nod == 5
    if dir == :x || dir == :r
     num[1] = iq*(2*nxe+1)+ip
     num[2] = (iq-1)*(2*nxe+1)+ip
     num[3] = num[2]+1
     num[4] = num[1]+1
     num[5] = iq*(2*nxe+1)+ip-nxe
    else
     num[1] = (ip-1)*(2*nye+1)+iq+1
     num[2] = num[1]-1
     num[3] = ip*(2*nye+1)+iq
     num[4] = num[3]+1
     num[5] = ip*(2*nye+1)+iq-nye
    end
    coord[1:2,1] = x_coords[ip]
    coord[3:4,1] = x_coords[ip+1]
    coord[1,2] = y_coords[iq+1]
    coord[2:3,2] = y_coords[iq]
    coord[4,2] = coord[1,2]
    coord[5,:] = 0.25_iwp*(coord[1,:]+coord[2,:]+coord[3,:]+coord[4,:])
  elseif nod == 8
    if dir == :x || dir == :r
     num[1] = iq*(3*nxe+2)+2*ip-1
     num[2] = iq*(3*nxe+2)+ip-nxe-1
     num[3] = (iq-1)*(3*nxe+2)+2*ip-1
     num[4] = num[3]+1
     num[5] = num[4]+1
     num[6] = num[2]+1
     num[7] = num[1]+2
     num[8] = num[1]+1
    else
     num[1] = (ip-1)*(3*nye+2)+2*iq+1
     num[2] = num[1]-1
     num[3] = num[1]-2
     num[4] = (ip-1)*(3*nye+2)+2*nye+iq+1
     num[5] = ip*(3*nye+2)+2*iq-1
     num[6] = num[5]+1
     num[7] = num[5]+2
     num[8] = num[4]+1
    end
    coord[1:3,1] .= x_coords[ip]
    coord[5:7,1] .= x_coords[ip+1]
    coord[4,1] = 0.5*(coord[3,1]+coord[5,1])
    coord[8,1] = 0.5*(coord[7,1]+coord[1,1])
    coord[1,2] = y_coords[iq+1]
    coord[7:8,2] .= y_coords[iq+1]
    coord[3:5,2] .= y_coords[iq]
    coord[2,2] = 0.5*(coord[1,2]+coord[3,2])
    coord[6,2] = 0.5*(coord[5,2]+coord[7,2])
  elseif nod == 9
    if dir == :x || dir == :r
     num[1] = iq*(4*nxe+2)+2*ip-1
     num[2] = iq*(4*nxe+2)+2*ip-nxe-4
     num[3] =  (iq-1)*(4*nxe+2)+2*ip-1
     num[4] = num[3]+1
     num[5] = num[4]+1
     num[6] = num[2]+2
     num[7] = num[1]+2
     num[8] = num[1]+1
     num[9] = num[2]+1
    else
     num[1] = (ip-1)*2*(2*nye+1)+2*iq+1
     num[2] = num[1]-1
     num[3] = num[1]-2
     num[4] = (ip-1)*2*(2*nye+1)+2*nye+2*iq
     num[5] = ip*2*(2*nye+1)+2*iq-1
     num[6] = num[5]+1
     num[7] = num[5]+2
     num[8] = num[4]+2
     num[9] = num[4]+1
    end
    coord[1:3,1] = x_coords[ip]
    coord[5:7,1] = x_coords[ip+1]
    coord[4,1] = 0.5*(coord[3,1]+coord[5,1])
    coord[8,1] = 0.5*(coord[7,1]+coord[1,1])
    coord[1,2] = y_coords[iq+1]
    coord[7:8,2] = y_coords[iq+1]
    coord[3:5,2] = y_coords[iq]
    coord[2,2] = 0.5*(coord[1,2]+coord[3,2])
    coord[6,2] = 0.5*(coord[5,2]+coord[7,2])
    coord[9,:] = 0.5*(coord[4,:]+coord[8,:])
  else
    println("Wrong number of nodes for quadrilateral fin_el.")
    exit(1)
  end
end


text:
import DataStructures: Deque

@testset "PriceLevel" begin
    pl = PriceLevel()

    # Test creation of price level
    @test size(orders(pl)) == (0,)
    @test size(pl) == 0

    # Test adding a limit order
    oid = UInt128(0)
    o = LimitOrder(UInt128(10), UInt128(1), UInt128(1), BUY, oid)
    add!(pl, o)
    @test size(orders(pl)) == (1,)
    @test size(pl) == 1

    # Test removing a limit order
    remove!(pl, oid)
    @test size(orders(pl)) == (0,)
    @test size(pl) == 0

    # Test that orders are in the correct order
    o2 = LimitOrder(UInt128(10), UInt128(1), UInt128(1), BUY, UInt128(2))
    add!(pl, o)
    add!(pl, o2)

    first_order = first(pl)
    @test first_order == o

    remove!(pl, id(o))

    first_order = first(pl)
    @test first_order == o2

    # Test to input many orders and remove them
    pl = PriceLevel()
    for i in 1:100
        add!(pl, LimitOrder(UInt128(10),
                            UInt128(1),
                            UInt128(1),
                            BUY, UInt128(i)))
    end

    @test size(pl) == 100
    @test size(orders(pl)) == (100,)

    for i in 1:100
        remove!(pl, UInt128(i))
    end

    @test size(pl) == 0
    @test size(orders(pl)) == (0,)

    @testset "Matching" begin
        pl = PriceLevel()
        # Adds three buy orders of size 10
        o1 = LimitOrder(UInt128(10), UInt128(10), UInt128(1), BUY, UInt128(1))
        add!(pl,o1)
        o2 = LimitOrder(UInt128(10), UInt128(10), UInt128(1), BUY, UInt128(2))
        add!(pl,o2)
        o3 = LimitOrder(UInt128(10), UInt128(10), UInt128(1), BUY, UInt128(3))
        add!(pl,o3)

        trades = Deque{LimitOrder}()
        size = UInt128(10)
        trades, size = match!(pl, size, trades)

        #expected_trades = Deque{LimitOrder}()
        #push!(expected_trades, o1)
        #@test trades === expected_trades
        # The entire sized was matched

        @test size == 0
        @test pl.size == 20

        trades, size = match!(pl, UInt128(15), trades)

        @test size == 0
        @test pl.size == 5

        trades, size = match!(pl, UInt128(10), trades)

        @test size == 5
        @test pl.size == 0

    end

end


text:
using Base.Test

addprocs(2)
using PmapProgressMeter
using ProgressMeter

# just make sure it runs
vals = 1:10
p = Progress(length(vals))
@test pmap(x->begin sleep(1); x*2 end, p, vals)[1] == vals[1]*2

# test the do-block syntax
pmap(p, vals) do x
    sleep(1)
    x*2
end

# make sure it runs it kwargs
vals = 1:10
p = Progress(length(vals))
@test pmap(x->begin sleep(.1); x*2 end, p, vals, err_stop=true)[1] == vals[1]*2

# try with multiple lists
vals2 = 10:-1:1
p = Progress(length(vals))
@test pmap(+, p, vals, vals2) == 11*ones(Int,length(vals))

# make sure callback passing works
vals = 1:10
@test pmap((cb, x) -> begin sleep(.1); cb(1); x*2 end, Progress(length(vals)),vals,passcallback=true)[1] == vals[1]*2

max_stars_repo_name: slundberg/PmapProgressMeter

max_issues_repo_name: slundberg/PmapProgressMeter

max_forks_repo_name: slundberg/PmapProgressMeter


text:
@test Canvas <: GraphicsArea

@testset "Nan / Inf" begin
    c = BrailleCanvas(40, 15)
    lines!(c, [0, NaN], [0, 0])
    lines!(c, [0, 0], [+Inf, 0])
    lines!(c, [0, 0], [0, -Inf])
end

@testset "interface" begin
    for (T, xres, yres) in [
        (BrailleCanvas, 2, 4),
        (DensityCanvas, 1, 2),
        (BlockCanvas, 2, 2),
        (AsciiCanvas, 3, 3),
        (DotCanvas, 1, 2),
        (HeatmapCanvas, 1, 2),
    ]
        @testset "$(nameof(T))" begin
            @test T <: Canvas
            c = @inferred T(30, 15, origin_x = -1, origin_y = -1.5, width = 2, height = 3)
            @test @inferred(size(c)) === (2.0, 3.0)
            @test @inferred(width(c)) === 2.0
            @test @inferred(height(c)) === 3.0
            @test @inferred(origin(c)) === (-1.0, -1.5)
            @test @inferred(origin_x(c)) === -1.0
            @test @inferred(origin_y(c)) === -1.5
            @test @inferred(ncols(c)) === 30
            @test @inferred(nrows(c)) === (T == HeatmapCanvas ? 8 : 15)
            @test @inferred(UnicodePlots.x_pixel_per_char(typeof(c))) === xres
            @test @inferred(UnicodePlots.y_pixel_per_char(typeof(c))) === yres
            @test @inferred(pixel_width(c)) === 30 * xres
            @test @inferred(pixel_height(c)) === 15 * yres
            @test @inferred(pixel_size(c)) === (pixel_width(c), pixel_height(c))
            if T <: UnicodePlots.LookupCanvas  # coverage
                @test length(UnicodePlots.lookup_encode(c)) > 0
                @test length(UnicodePlots.lookup_decode(c)) > 0
            end
        end
    end
end

@testset "print and show" begin
    seed!(RNG, 1337)
    x1, y1 = rand(RNG, 20), rand(RNG, 20)
    x2, y2 = rand(RNG, 50), rand(RNG, 50)
    for (T, str) in [
        (BrailleCanvas, "braille"),
        (DensityCanvas, "density"),
        (BlockCanvas, "block"),
        (AsciiCanvas, "ascii"),
        (DotCanvas, "dot"),
        (HeatmapCanvas, "heatmap"),
    ]
        @testset "$(nameof(T))" begin
            c = T(40, 10, origin_x = 0.0, origin_y = 0.0, width = 1.0, height = 1.0)
            if T == BrailleCanvas
                test_ref("references/canvas/empty_braille_show.txt", @show_col(c))
            elseif T == HeatmapCanvas
                test_ref("references/canvas/empty_heatmap_show.txt", @show_col(c))
            else
                test_ref("references/canvas/empty_show.txt", @show_col(c))
            end
            @test @inferred(lines!(c, 0.0, 0.0, 1.0, 1.0, :blue)) === c
            @test @inferred(points!(c, x1, y1, :white)) === c
            @test @inferred(pixel!(c, 2, 4, color = :cyan)) === c
            points!(c, x2, y2, color = :red)
            lines!(c, 0.0, 1.0, 0.5, 0.0, :green)
            points!(c, 0.05, 0.3, color = :cyan)
            lines!(c, [1.0, 2], [2.0, 1])
            lines!(c, 0.0, 0.0, 0.9, 9999.0, :yellow)
            lines!(c, 0, 0, 1, 1, color = :blue)
            lines!(c, 0.1, 0.7, 0.9, 0.6, :red)
            test_ref(
                "references/canvas/$(str)_printrow.txt",
                @io2str(printrow(IOContext(::IO, :color => true), c, 3))
            )
            test_ref("references/canvas/$(str)_print.txt", @print_col(c))
            test_ref("references/canvas/$(str)_print_nocolor.txt", @print_nocol(c))
            test_ref("references/canvas/$(str)_show.txt", @show_col(c))
            test_ref("references/canvas/$(str)_show_nocolor.txt", @show_nocol(c))
        end
    end
end

@testset "color mixing" begin
    c = BrailleCanvas(40, 15, origin_x = 0.0, origin_y = 0.0, width = 1.2, height = 1.2)
    @inferred lines!(c, 0.0, 0.0, 1.2, 0.0, :blue)
    lines!(c, 0.0, 0.2, 1.2, 0.2, :red)
    lines!(c, 0.0, 0.4, 1.2, 0.4, :green)
    lines!(c, 0.0, 0.6, 1.2, 0.6, :yellow)
    lines!(c, 0.0, 0.8, 1.2, 0.8, :magenta)
    lines!(c, 0.0, 1.0, 1.2, 1.0, :cyan)
    lines!(c, 0.0, 1.2, 1.2, 1.2, :white)

    lines!(c, 0.0, 0.0, 0.0, 1.2, :blue)
    lines!(c, 0.2, 0.0, 0.2, 1.2, :red)
    lines!(c, 0.4, 0.0, 0.4, 1.2, :green)
    lines!(c, 0.6, 0.0, 0.6, 1.2, :yellow)
    lines!(c, 0.8, 0.0, 0.8, 1.2, :magenta)
    lines!(c, 1.0, 0.0, 1.0, 1.2, :cyan)
    lines!(c, 1.2, 0.0, 1.2, 1.2, :normal)
    test_ref("references/canvas/color_mixing.txt", @print_col(c))
end

@testset "colors" begin
    c = BrailleCanvas(40, 15, origin_x = 0.0, origin_y = 0.0, width = 1.2, height = 1.2)
    @inferred lines!(c, 0.0, 0.0, 1.2, 0.0, (0, 0, 255))
    lines!(c, 0.0, 0.2, 1.2, 0.2, (255, 0, 0))
    lines!(c, 0.0, 0.4, 1.2, 0.4, (0, 255, 0))
    lines!(c, 0.0, 0.6, 1.2, 0.6, (255, 255, 0))
    lines!(c, 0.0, 0.8, 1.2, 0.8, 127)
    lines!(c, 0.0, 1.0, 1.2, 1.0, 51)
    lines!(c, 0.0, 1.2, 1.2, 1.2, 231)

    lines!(c, 0.0, 0.0, 0.0, 1.2, :blue)
    lines!(c, 0.2, 0.0, 0.2, 1.2, :red)
    lines!(c, 0.4, 0.0, 0.4, 1.2, :green)
    lines!(c, 0.6, 0.0, 0.6, 1.2, :yellow)
    lines!(c, 0.8, 0.0, 0.8, 1.2, :magenta)
    lines!(c, 1.0, 0.0, 1.0, 1.2, :cyan)
    lines!(c, 1.2, 0.0, 1.2, 1.2, :normal)
    test_ref("references/canvas/colors.txt", @print_col(c))
end


text:
using BenchmarkTools
using Yao.Intrisics

const nbit = 16
const V = randn(ComplexF64, 1<<nbit)
const res4 = copy(V)

it = itercontrol(nbit, [3],[1])
@benchmark controldo(x->swaprows!($res4, x+1, x-3), $it)
@benchmark for i in $it swaprows!($res4, i+1, i-3) end
@benchmark controldo($(x->x), $it)


text:
args = Dict{Any,Any}()
args["file"] = String(ARGS[1])
args["h_or_u"] = String(ARGS[2])
args["hcase"] = String(ARGS[3])
branches_h = Dict{Any,Any}()
branches = Dict{Any,Any}()

for k in 2:15
    if k <= 10
        f = open("$(args["file"])_$k")
        lines = readlines(f)
        branches[k] = Dict{AbstractString, Any}("soc" => Set(sort([parse(Int, branch) for branch in split(split(split(lines[4], "[")[2], "]")[1], ",")])))
    else
        f = open("$(args["file"])_$(k)_nonf")
        lines = readlines(f)
        branches[k] = Dict{AbstractString, Any}("soc" => Set(sort([parse(Int, branch) for branch in split(split(split(lines[3], "[")[2], "]")[1], ",")])))
    end
    f_h = open("$(args["file"])_$(k)_soc_$(args["h_or_u"])$(args["hcase"])")
    lines_h = readlines(f_h)
    branches_h[k] = Dict{AbstractString, Any}("soc" => Set(sort([parse(Int, branch) for branch in split(split(split(lines_h[2], "[")[2], "]")[1], ",")])))
end

hamming_distance = zeros(14)

for k in 2:15
    hamming_distance[k-1] = length(intersect(branches_h[k]["soc"], branches[k]["soc"]))
end

for i in hamming_distance
    println(i)
end

max_stars_repo_name: kaarthiksundar/acopf-n-k

max_issues_repo_name: kaarthiksundar/acopf-n-k

max_forks_repo_name: kaarthiksundar/acopf-n-k


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

## core string functions ##

endof(s::AbstractString) = error("you must implement endof(", typeof(s), ")")
next(s::AbstractString, i::Int) = error("you must implement next(", typeof(s), ",Int)")
next(s::DirectIndexString, i::Int) = (s[i],i+1)
next(s::AbstractString, i::Integer) = next(s,Int(i))

string() = ""
string(s::AbstractString) = s

"""
    String(s::AbstractString)

Convert a string to a contiguous byte array representation encoded as UTF-8 bytes.
This representation is often appropriate for passing strings to C.
"""
String(s::AbstractString) = print_to_string(s)
String(s::String) = s

# String constructor docstring from boot.jl, workaround for #16730
# and the unavailability of @doc in boot.jl context.
"""
    String(v::Vector{UInt8})

Create a new `String` from a vector `v` of bytes containing
UTF-8 encoded characters.   This function takes "ownership" of
the array, which means that you should not subsequently modify
`v` (since strings are supposed to be immutable in Julia) for
as long as the string exists.

If you need to subsequently modify `v`, use `String(copy(v))` instead.
"""
String(v::Array{UInt8,1})

"""
    unsafe_string(p::Ptr{UInt8}, [length::Integer])

Copy a string from the address of a C-style (NUL-terminated) string encoded as UTF-8.
(The pointer can be safely freed afterwards.) If `length` is specified
(the length of the data in bytes), the string does not have to be NUL-terminated.

This function is labelled "unsafe" because it will crash if `p` is not
a valid memory address to data of the requested length.

See also [`unsafe_wrap(String, p, [length])`](:func:`unsafe_wrap`), which takes a pointer
and wraps a string object around it without making a copy.
"""
function unsafe_string(p::Union{Ptr{UInt8},Ptr{Int8}}, len::Integer)
    p == C_NULL && throw(ArgumentError("cannot convert NULL to string"))
    ccall(:jl_pchar_to_string, Ref{String}, (Ptr{UInt8},Int), p, len)
end
function unsafe_string(p::Union{Ptr{UInt8},Ptr{Int8}})
    p == C_NULL && throw(ArgumentError("cannot convert NULL to string"))
    ccall(:jl_cstr_to_string, Ref{String}, (Ptr{UInt8},), p)
end

convert(::Type{Vector{UInt8}}, s::AbstractString) = String(s).data
convert(::Type{Array{UInt8}}, s::AbstractString) = String(s).data
convert(::Type{String}, s::AbstractString) = String(s)
convert(::Type{Vector{Char}}, s::AbstractString) = collect(s)
convert(::Type{Symbol}, s::AbstractString) = Symbol(s)
convert(::Type{String}, s::Symbol) = unsafe_string(Cstring(s))

## generic supplied functions ##

start(s::AbstractString) = 1
done(s::AbstractString,i) = (i > endof(s))
getindex(s::AbstractString, i::Int) = next(s,i)[1]
getindex(s::AbstractString, i::Integer) = s[Int(i)]
getindex(s::AbstractString, i::Colon) = s
getindex{T<:Integer}(s::AbstractString, r::UnitRange{T}) = s[Int(first(r)):Int(last(r))]
# TODO: handle other ranges with stride ±1 specially?
getindex(s::AbstractString, v::AbstractVector) =
    sprint(length(v), io->(for i in v; write(io,s[i]) end))

Symbol(s::AbstractString) = Symbol(String(s))

"""
    sizeof(s::AbstractString)

The number of bytes in string `s`.

```jldoctest
julia> sizeof("❤")
3
```
"""
sizeof(s::AbstractString) = error("type $(typeof(s)) has no canonical binary representation")

eltype{T<:AbstractString}(::Type{T}) = Char

"""
```
*(s::AbstractString, t::AbstractString)
```

Concatenate strings. The `*` operator is an alias to this function.

```jldoctest
julia> "Hello " * "world"
"Hello world"
```
"""
(*)(s1::AbstractString, ss::AbstractString...) = string(s1, ss...)
(.*){T<:AbstractString}(v::Vector{T},s::AbstractString) = [i*s for i in v]
(.*){T<:AbstractString}(s::AbstractString,v::Vector{T}) = [s*i for i in v]

length(s::DirectIndexString) = endof(s)

"""
    length(s::AbstractString)

The number of characters in string `s`.
"""
function length(s::AbstractString)
    i = start(s)
    if done(s,i)
        return 0
    end
    n = 1
    while true
        c, j = next(s,i)
        if done(s,j)
            return n
        end
        n += 1
        i = j
    end
end

## string comparison functions ##

function cmp(a::AbstractString, b::AbstractString)
    if a === b
        return 0
    end
    i = start(a)
    j = start(b)
    while !done(a,i) && !done(b,i)
        c, i = next(a,i)
        d, j = next(b,j)
        if c != d
            return c < d ? -1 : +1
        end
    end
    done(a,i) && !done(b,j) ? -1 :
    !done(a,i) && done(b,j) ? +1 : 0
end

==(a::AbstractString, b::AbstractString) = cmp(a,b) == 0
isless(a::AbstractString, b::AbstractString) = cmp(a,b) < 0

# faster comparisons for byte strings and symbols

cmp(a::String, b::String) = lexcmp(a.data, b.data)
cmp(a::Symbol, b::Symbol) = Int(sign(ccall(:strcmp, Int32, (Cstring, Cstring), a, b)))

==(a::String, b::String) = a.data == b.data
isless(a::Symbol, b::Symbol) = cmp(a,b) < 0

## Generic validation functions ##

isvalid(s::DirectIndexString, i::Integer) = (start(s) <= i <= endof(s))

"""
    isvalid(str::AbstractString, i::Integer)

Tells whether index `i` is valid for the given string.
"""
function isvalid(s::AbstractString, i::Integer)
    i < 1 && return false
    done(s,i) && return false
    try
        next(s,i)
        true
    catch
        false
    end
end

## Generic indexing functions ##

prevind(s::DirectIndexString, i::Integer) = Int(i)-1
prevind(s::AbstractArray    , i::Integer) = Int(i)-1
nextind(s::DirectIndexString, i::Integer) = Int(i)+1
nextind(s::AbstractArray    , i::Integer) = Int(i)+1

function prevind(s::String, i::Integer)
    j = Int(i)
    e = endof(s.data)
    if j > e
        return endof(s)
    end
    j -= 1
    while j > 0 && is_valid_continuation(s.data[j])
        j -= 1
    end
    j
end

function nextind(s::String, i::Integer)
    j = Int(i)
    if j < 1
        return 1
    end
    e = endof(s.data)
    j += 1
    while j <= e && is_valid_continuation(s.data[j])
        j += 1
    end
    j
end

"""
    prevind(str::AbstractString, i::Integer)

Get the previous valid string index before `i`.
Returns a value less than `1` at the beginning of the string.
"""
function prevind(s::AbstractString, i::Integer)
    e = endof(s)
    if i > e
        return e
    end
    j = Int(i)-1
    while j >= 1
        if isvalid(s,j)
            return j
        end
        j -= 1
    end
    return 0 # out of range
end

"""
    nextind(str::AbstractString, i::Integer)

Get the next valid string index after `i`.
Returns a value greater than `endof(str)` at or after the end of the string.
"""
function nextind(s::AbstractString, i::Integer)
    e = endof(s)
    if i < 1
        return 1
    end
    if i > e
        return Int(i)+1
    end
    for j = Int(i)+1:e
        if isvalid(s,j)
            return j
        end
    end
    next(s,e)[2] # out of range
end

checkbounds(s::AbstractString, i::Integer) = start(s) <= i <= endof(s) || throw(BoundsError(s, i))
checkbounds{T<:Integer}(s::AbstractString, r::Range{T}) = isempty(r) || (minimum(r) >= start(s) && maximum(r) <= endof(s)) || throw(BoundsError(s, r))
# The following will end up using a deprecated checkbounds, when T is not Integer
checkbounds{T<:Real}(s::AbstractString, I::AbstractArray{T}) = all(i -> checkbounds(s, i), I)
checkbounds{T<:Integer}(s::AbstractString, I::AbstractArray{T}) = all(i -> checkbounds(s, i), I)

ind2chr(s::DirectIndexString, i::Integer) = begin checkbounds(s,i); i end
chr2ind(s::DirectIndexString, i::Integer) = begin checkbounds(s,i); i end

"""
    ind2chr(s::AbstractString, i::Integer)

Convert a byte index `i` to a character index with
respect to string `s`.
"""
function ind2chr(s::AbstractString, i::Integer)
    s[i] # throws error if invalid
    j = 1
    k = start(s)
    while true
        c, l = next(s,k)
        if i <= k
            return j
        end
        j += 1
        k = l
    end
end

"""
    chr2ind(s::AbstractString, i::Integer)

Convert a character index `i` to a byte index.
"""
function chr2ind(s::AbstractString, i::Integer)
    i < start(s) && throw(BoundsError(s, i))
    j = 1
    k = start(s)
    while true
        c, l = next(s,k)
        if i == j
            return k
        end
        j += 1
        k = l
    end
end

immutable EachStringIndex{T<:AbstractString}
    s::T
end
eachindex(s::AbstractString) = EachStringIndex(s)

length(e::EachStringIndex) = length(e.s)
start(e::EachStringIndex) = start(e.s)
next(e::EachStringIndex, state) = (state, nextind(e.s, state))
done(e::EachStringIndex, state) = done(e.s, state)
eltype(::Type{EachStringIndex}) = Int

typealias Chars Union{Char,Tuple{Vararg{Char}},AbstractVector{Char},Set{Char}}

typealias ByteArray Union{Vector{UInt8},Vector{Int8}}

## character column width function ##

"""
    strwidth(s::AbstractString)

Gives the number of columns needed to print a string.

```jldoctest
julia> strwidth("March")
5
```
"""
strwidth(s::AbstractString) = (w=0; for c in s; w += charwidth(c); end; w)

"""
    isascii(c::Union{Char,AbstractString}) -> Bool

Tests whether a character belongs to the ASCII character set, or whether this is true for
all elements of a string.
"""
isascii(c::Char) = c < Char(0x80)
isascii(s::AbstractString) = all(isascii, s)

## string promotion rules ##

promote_rule{S<:AbstractString,T<:AbstractString}(::Type{S}, ::Type{T}) = String

"""
    isxdigit(c::Union{Char,AbstractString}) -> Bool

Tests whether a character is a valid hexadecimal digit, or whether this is true for all elements of a string.

```jldoctest
julia> isxdigit("abc")
true

julia> isxdigit("0x9")
false
```
"""
isxdigit(c::Char) = '0'<=c<='9' || 'a'<=c<='f' || 'A'<=c<='F'
isxdigit(s::AbstractString) = all(isxdigit, s)

## checking UTF-8 & ACSII validity ##

byte_string_classify(data::Vector{UInt8}) =
    ccall(:u8_isvalid, Int32, (Ptr{UInt8}, Int), data, length(data))
byte_string_classify(s::String) = byte_string_classify(s.data)
    # 0: neither valid ASCII nor UTF-8
    # 1: valid ASCII
    # 2: valid UTF-8

isvalid(::Type{String}, s::Union{Vector{UInt8},String}) = byte_string_classify(s) != 0
isvalid(s::String) = isvalid(String, s)

## uppercase and lowercase transformations ##

"""
    uppercase(s::AbstractString)

Returns `s` with all characters converted to uppercase.

```jldoctest
julia> uppercase("Julia")
"JULIA"
```
"""
uppercase(s::AbstractString) = map(uppercase, s)

"""
    lowercase(s::AbstractString)

Returns `s` with all characters converted to lowercase.

```jldoctest
julia> lowercase("STRINGS AND THINGS")
"strings and things"
```
"""
lowercase(s::AbstractString) = map(lowercase, s)

"""
    ucfirst(s::AbstractString)

Returns `string` with the first character converted to uppercase.

```jldoctest
julia> ucfirst("python")
"Python"
```
"""
function ucfirst(s::AbstractString)
    isempty(s) || isupper(s[1]) ? s : string(uppercase(s[1]),s[nextind(s,1):end])
end

"""
    lcfirst(s::AbstractString)

Returns `string` with the first character converted to lowercase.

```jldoctest
julia> lcfirst("Julia")
"julia"
```
"""
function lcfirst(s::AbstractString)
    isempty(s) || islower(s[1]) ? s : string(lowercase(s[1]),s[nextind(s,1):end])
end

## string map, filter, has ##

function map(f, s::AbstractString)
    out = IOBuffer(Array{UInt8}(endof(s)),true,true)
    truncate(out,0)
    for c in s
        c2 = f(c)
        if !isa(c2,Char)
            throw(ArgumentError("map(f,s::AbstractString) requires f to return Char; try map(f,collect(s)) or a comprehension instead"))
        end
        write(out, c2::Char)
    end
    String(take!(out))
end

function filter(f, s::AbstractString)
    out = IOBuffer(Array{UInt8}(endof(s)),true,true)
    truncate(out,0)
    for c in s
        if f(c)
            write(out, c)
        end
    end
    String(take!(out))
end

max_stars_repo_name: habemus-papadum/julia

max_issues_repo_name: habemus-papadum/julia

max_forks_repo_name: habemus-papadum/julia


text:
using SnoopCompile

### Log the compiles
# This only needs to be run once (to generate "/tmp/images_compiles.csv")

SnoopCompile.@snoop "/tmp/cudart_compiles.csv" begin
    include(joinpath(dirname(@__FILE__), "..", "test","runtests.jl"))
end

### Parse the compiles and generate precompilation scripts
# This can be run repeatedly

# IMPORTANT: we must have the module(s) defined for the parcelation
# step, otherwise we will get no precompiles for the Images module
using CUDArt

data = SnoopCompile.read("/tmp/cudart_compiles.csv")

# The Images tests are run inside a module ImagesTest, so all
# the precompiles get credited to ImagesTest. Credit them to Images instead.
subst = Dict() #Dict("ImagesTests"=>"Images")

# Blacklist helps fix problems:
# - MIME uses type-parameters with symbols like :image/png, which is
#   not parseable
blacklist = ["MIME", "__init__", "stream"]

# Use these two lines if you want to create precompile functions for
# individual packages
pc, discards = SnoopCompile.parcel(data[end:-1:1,2], subst=subst, blacklist=blacklist)
SnoopCompile.write("/tmp/precompile", pc)


text: test_gage_principal_axes([1. 0. 0.;], eye(3), "unit vector along X")


text:
include("multpack.jl")
include("multtest.jl")

module MultExperimental

export enkst, sqsum
export sm, sm1, smn, smn1, sv, svref

import MultiPoly
using MultiPoly
import Polynomials
using Polynomials
import MultRoot
using MultRoot

n, k, s, t = generators(MPoly{Float64}, :n, :k, :s, :t)

A = [n (n+1)*n/2; -(n-(k-1)/2)*k -(n*k-(k-1)*(4k+1)/6)*k]

function enkst(A,n,k,s,t)
    enk(a) = evaluate(a, n, k)
	ank = map(enk, A)
	ank \ [s; t]
end

function sqsum(a::AbstractFloat, b::AbstractFloat, xi, yi, k::Int)

	s1 = sum( (a - b*xi[1:k] -yi[1:k])^2 )
	s2 = sum( (a - b*k -yi[1:k])^2 )
	s = (s1+s2) / 2
end

function sqsum(xi, yi, k::Int)
    n = length(xi)
	if n != length(yi) throw(ArgumentError("length mismatch")) end

    ysum = sum(yi)
	yxsum = sum(xi .* yi)
	ab = enkst(A, n, k, ysum, yxsum)
	a, b = ab[1], ab[2]
	sqsum(a, b, xi, yi, k)
end

sm(p::Poly, m::Int64) = sylves(reverse(p.a), reverse((p').a)/degree(p), m)
sm1(p::Poly, m::Int64) = sylves1(reverse(p.a), reverse((p').a)/degree(p), m)

function scalerow(s::Array{Float64,2})
	m = size(s,1)
	b = Float64[1.0 / norm(s[i,:]) for i = 1:m]
	scale(b, s)
end

smn(p::Poly, m::Int64) = scalerow(sm(p, m))
smn1(p::Poly, m::Int64) = scalerow(sm1(p, m))

sv(p::Poly, m::Int) = zminsv(sm(p,m), 1e-10)

sv(p::Poly) = m::Int -> sv(p,m)[1]

function svref(p::Poly, m::Int64)
    A = sm(p, m)
	s0, x0 = zminsv(A, 1e-10)
	f = reverse(p.a)
	g = reverse((p').a)/(length(p.a)-2)
	h, u, v, res0, res, s1 = MultRoot.gcd_refinement(x0, m, f, g, A)
	s0, s1, res0, res
end

svref(p::Poly) = m::Int -> svref(p,m)

function qrref(p::Poly, m::Int64)
    A = smn1(p, m)
	Q, R = qr(A)
	D = abs(diag(R))
	minimum(D)
end

qrref(p::Poly) = m::Int64 -> qrref(p, m)

end

import MultExperimental
using MultExperimental


text:
using DataFrames, ReusePatterns

struct DataFrameMeta <: AbstractDataFrame
    p::DataFrame
    meta::Dict{String, Any}
    DataFrameMeta(args...; kw...) = new(DataFrame(args...; kw...), Dict{Symbol, Any}())
end
meta(d::DataFrameMeta) = getfield(d,:meta)  # <-- new functionality added to DataFrameMeta
@forward((DataFrameMeta, :p), DataFrame)    # <-- reuse all existing functionalities

# Use a `DataFrameMeta` object as if it was a common DataFrame object ...
v = ["x","y","z"][rand(1:3, 10)]
df1 = DataFrameMeta(Any[collect(1:10), v, rand(10)], [:A, :B, :C])
df2 = DataFrameMeta(A = 1:10, B = v, C = rand(10))
dump(df1)
dump(df2)
describe(df2)
first(df1, 10)
df1[:, :A] .+ df2[:, :C]
df1[1:4, 1:2]
df1[:, [:A,:C]]
df1[1:2, [:A,:C]]
df1[:, [1,3]]
df1[1:4, :]
df1[1:4, :C]
df1[1:4, :C] = 40. * df1[1:4, :C]
[df1; df2]
size(df1)
meta(df1)["key"] = :value
meta(df1)["key"]