text:
__precompile__()

module Syslogs

export Syslog

using Sockets
using Printf

# Default UDP and TCP ports
const UDP_PORT = 514
const TCP_PORT = 514

const LOG_EMERG   = 0  # system is unusable
const LOG_ALERT   = 1  # action must be taken immediately
const LOG_CRIT    = 2  # critical conditions
const LOG_ERR     = 3  # error conditions
const LOG_WARNING = 4  # warning conditions
const LOG_NOTICE  = 5  # normal but significant condition
const LOG_INFO	  = 6  # informational
const LOG_DEBUG   = 7  # debug-level messages

#  facility codes
const LOG_KERN     = 0  # kernel messages
const LOG_USER     = 1  # random user-level messages
const LOG_MAIL     = 2  # mail system
const LOG_DAEMON   = 3  # system daemons
const LOG_AUTH     = 4  # security/authorization messages
const LOG_SYSLOG   = 5  # messages generated internally by syslogd
const LOG_LPR      = 6  # line printer subsystem
const LOG_NEWS     = 7  # network news subsystem
const LOG_UUCP     = 8  # UUCP subsystem
const LOG_CRON     = 9  # clock daemon
const LOG_AUTHPRIV = 10	# security/authorization messages (private)

#  other codes through 15 reserved for local use
const LOG_LOCAL0 = 16
const LOG_LOCAL1 = 17
const LOG_LOCAL2 = 18
const LOG_LOCAL3 = 19
const LOG_LOCAL4 = 20
const LOG_LOCAL5 = 21
const LOG_LOCAL6 = 22
const LOG_LOCAL7 = 23

# Define dictionary mappings from keywords to codes
# for levels and facilities.
const LEVELS = Dict(
    :emergency => LOG_EMERG,
    :emerg => LOG_EMERG,
    :alert => LOG_ALERT,
    :critical => LOG_CRIT,
    :crit => LOG_CRIT,
    :error => LOG_ERR,
    :err => LOG_ERR,
    :warning => LOG_WARNING,
    :warn => LOG_WARNING,
    :notice => LOG_NOTICE,
    :info => LOG_INFO,
    :debug => LOG_DEBUG,
)

# Define valid syslog facilities (even kern, which isn't accessible from userspace).
const FACILITIES = Dict(
    :kern => LOG_KERN,
    :user => LOG_USER,
    :mail => LOG_MAIL,
    :daemon => LOG_DAEMON,
    :auth => LOG_AUTH,
    :security => LOG_AUTH,
    :syslog => LOG_SYSLOG,
    :lpr => LOG_LPR,
    :news => LOG_NEWS,
    :uucp => LOG_UUCP,
    :cron => LOG_CRON,
    :authpriv => LOG_AUTHPRIV,
    :local0 => LOG_LOCAL0,
    :local1 => LOG_LOCAL1,
    :local2 => LOG_LOCAL2,
    :local3 => LOG_LOCAL3,
    :local4 => LOG_LOCAL4,
    :local5 => LOG_LOCAL5,
    :local6 => LOG_LOCAL6,
    :local7 => LOG_LOCAL7,
)

openlog(ident::String, logopt::Integer, facility::Integer) =
    ccall(:openlog, Cvoid, (Ptr{UInt8}, UInt, UInt), ident, logopt, facility)

syslog(priority::Integer, msg::String) =
    ccall(:syslog, Cvoid, (UInt, Ptr{UInt8}), priority, msg)

closelog() = ccall(:closelog, Cvoid, ())

makepri(facility::Integer, priority::Integer) = (UInt(facility) << 3) | UInt(priority)

"""
    Syslog <: IO

`Syslog` handles writing logs to local (libc inteface) and remote (UDP / TCP sockets) syslog servers.

# Fields
* `address::Union{Tuple{IPAddr, Int}, Nothing}`: The host and ip for the remote servers. Logs locally if `nothing`.
* `facility::Symbol`: The syslog facility to write to (e.g., :local0, :ft, :daemon, etc) (defaults to :user)
* `socket::Union{Base.LibuvStream, Nothing}`: A UDP or TCP socket for logging messages. Logs locally if `nothing`.
"""
mutable struct Syslog <: IO
    address::Union{Tuple{IPAddr, Int}, Nothing}
    facility::Symbol
    socket::Union{Base.LibuvStream, Nothing}
end

Syslog(host::IPAddr, port::Int, facility::Symbol=:user; tcp::Bool=false) =
    Syslog((host, port), facility, tcp ? connect(host, port) : UDPSocket())

Syslog(host::IPAddr, facility::Symbol=:user; tcp::Bool=false) =
    Syslog(host, tcp ? TCP_PORT : UDP_PORT, facility; tcp=tcp)

Syslog(facility::Symbol=:user) = Syslog(nothing, facility, nothing)

"""
    println(::Syslog, ::AbstractString, ::AbstractString)

Converts the first AbstractString to a Symbol and call
`println(::Syslog, ::Symbol, ::AbstractString)`
"""
function Base.println(io::Syslog, level::AbstractString, msg::AbstractString)
    println(io, Symbol(lowercase(level)), msg)
end

function Base.println(io::Syslog, level::Symbol, msg::String)
    haskey(LEVELS, level) || throw(ArgumentError("Invalid logging level: $level."))
    haskey(FACILITIES, io.facility) || throw(ArgumentError("Invalid logging facility: $(io.facility)."))

    if io.socket !== nothing
        sock = io.socket

        if !(isa(sock, UDPSocket) || isa(sock, TCPSocket))
            throw(ArgumentError("Syslog only supports UDP or TCP sockets."))
        end
    end

    pri = makepri(FACILITIES[io.facility], LEVELS[level])

    if io.socket === nothing && io.address === nothing
        syslog(pri, msg)
    else
        sock = io.socket
        addr = io.address
        content = @sprintf("<%d>%s\0", pri, msg)

        if isa(sock, UDPSocket)
            # info("Sending ($(typeof(sock))): $content")
            send(sock, addr[1], addr[2], content)
        else
            # info("Sending ($(typeof(sock))): $content")
            try
                write(sock, content)
            catch e
                close(sock)
                s = connect(addr[1], addr[2])
                io.socket = s
                write(s, content)
            end
        end
    end
end

Base.log(io::Syslog, args...) = println(io, args...)

function Base.flush(io::Syslog)
    if io.socket !== nothing && isa(io.socket, TCPSocket)
        flush(io.socket)
    end
end

function Base.close(io::Syslog)
    if io.socket !== nothing
        close(io.socket)
    end
end

end  # module


text:
function run(s::String)::Int
    # Your code here
    numbers = Dict{Int,Int}()

    start = parse.(Int, split(s, ","))

    for idx in 1:length(start)
        numbers[start[idx]] = idx
    end

    last_seen_number::Int = start[end]
    current_number::Int = 0
    for idx::Int in length(start):30000000 - 1
        current_number = idx - get(numbers, last_seen_number, idx)
        numbers[last_seen_number] = idx
        last_seen_number = current_number
    end
    return last_seen_number
end

#########################################

function main()
    res, time, memory = @timed run(ARGS[1])
    println("_duration:$(time * 1000)")
    println(res)
end

main()

ext: jl


text:
mutable struct NLOpt
    opt
    x0
end
"""
    nlopt(lower, upper;
          x0 = rand(length(lower)) .* (upper .- lower) .+ lower,
          method = :GN_ESCH, kwargs...)

Resonable methods are `:GN_CRS2_LM`, `:GN_ESCH`, `:GN_ISRES`, `:LN_SBPLX`.
Other options like `maxtime` etc. can be set as keyword arguments
(see https://github.com/JuliaOpt/NLopt.jl).
"""
nlopt(prior; kwargs...) = nlopt(lower(prior), upper(prior); kwargs...)
function nlopt(lower, upper;
               x0 = rand(length(lower)) .* (upper .- lower) .+ lower,
               method = :GN_ESCH, kwargs...)
    D = length(lower)
    opt = NLopt.Opt(method, D)
    for (option, value) in pairs(kwargs)
        setproperty!(opt, option, value)
    end
    NLopt.lower_bounds!(opt, lower)
    NLopt.upper_bounds!(opt, upper)
    NLOpt(opt, x0)
end
export nlopt
function _optimize(o::NLOpt, f; maxfevals, verbose)
    o.opt.min_objective = (x, g) -> f(x)
    o.opt.maxeval = maxfevals
    res = NLopt.optimize(o.opt, o.x0)
    (x = res[2], f = res[1], res = res)
end


text:
module pushinterface_displayinterface_test_framebuffer

using Test

import PushInterface: DisplayInterface
import .DisplayInterface: rgb_to_pixel

@test rgb_to_pixel(1, 0xff, 0xff, 0xff, 0x00) == (0x18, 0xf4)
@test rgb_to_pixel(2, 0xff, 0xff, 0xff, 0x00) == (0x18, 0xf8)

@test rgb_to_pixel(1, 0xff, 0xd0, 0x02, 0x1b) == (0xfd, 0xeb)
@test rgb_to_pixel(2, 0xff, 0xd0, 0x02, 0x1b) == (0xfd, 0xe7)
@test rgb_to_pixel(3, 0xff, 0xff, 0xff, 0x00) == (0x18, 0xf4)
@test rgb_to_pixel(4, 0xff, 0xff, 0xff, 0x00) == (0x18, 0xf8)

import .DisplayInterface: Canvas, fill!, draw!
import .DisplayInterface: framebuffer
import Colors: @colorant_str
using Pkg

buf = reshape(read(joinpath(@__DIR__, "sample.bytes")), (20480, 16))
expected = vec(buf')

path = joinpath(@__DIR__, "line.png")
canvas = Canvas()
fill!(canvas, colorant"yellow") # @colorant_str
draw!(canvas, path)

a = framebuffer(canvas)
@test expected[1:16]             == UInt8[0xfd, 0xeb, 0xfd, 0xe7, 0x18, 0xf4, 0x18, 0xf8, 0x18, 0xf4, 0xfd, 0xe7, 0xfd, 0xeb, 0x18, 0xf8]
@test expected[16+1:16+16]       == UInt8[0x18, 0xf4, 0x18, 0xf8, 0x18, 0xf4, 0x18, 0xf8, 0x18, 0xf4, 0x18, 0xf8, 0x18, 0xf4, 0x18, 0xf8]
@test expected[2048+1 : 2048+16] == a[2048+1 : 2048+16]
@test expected                   == a

end # module pushinterface_displayinterface_test_framebuffer


text:
module TestExamples3DMHD

using Test
using Trixi

include("test_trixi.jl")

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "tree_3d_dgsem")

@testset "MHD" begin
  @trixi_testset "elixir_mhd_ec.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_ec.jl"),
      l2   = [0.017590099293094203, 0.017695875823827714, 0.017695875823827686, 0.017698038279620777, 0.07495006099352074, 0.010391801950005755, 0.010391801950005759, 0.010393502246627087, 2.524766553484067e-16],
      linf = [0.28173002819718196, 0.3297583616136297, 0.32975836161363004, 0.356862935505337, 1.2893514981209626, 0.10950981489747313, 0.10950981489747136, 0.11517234329681891, 2.0816911067714202e-15])
  end

  @trixi_testset "elixir_mhd_ec.jl with initial_condition=initial_condition_constant" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_ec.jl"),
      l2   = [4.270231310667203e-16, 2.4381208042014784e-15, 5.345107673575357e-15, 3.00313882171883e-15, 1.7772703118758417e-14, 1.0340110783830874e-15, 1.1779095371939702e-15, 9.961878521814573e-16, 8.1201730630719145e-16],
      linf = [2.4424906541753444e-15, 2.881028748902281e-14, 2.4646951146678475e-14, 2.3092638912203256e-14, 2.3447910280083306e-13, 1.7763568394002505e-14, 1.0436096431476471e-14, 2.042810365310288e-14, 7.057203733035201e-15],
      atol = 1000*eps(),
      initial_condition=initial_condition_constant)
  end

  @trixi_testset "elixir_mhd_alfven_wave.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave.jl"),
      l2   = [0.0032216851061739474, 0.009667756811390567, 0.00426372333971404, 0.011792748908209512, 0.009561034484102497, 0.009417462436849494, 0.005103363929916489, 0.01169498223549228, 0.003998265503740793],
      linf = [0.011829165032604894, 0.03710399824480111, 0.015496687949451582, 0.04700041585617602, 0.03190969001885602, 0.03490056585762913, 0.011667706905677644, 0.049549546108596104, 0.018477076688890922])
  end

  @trixi_testset "elixir_mhd_alfven_wave.jl with flux_derigs_etal" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave.jl"),
      l2   = [0.003872168982949962, 0.009234032031394519, 0.004166855440339939, 0.01160194799648245, 0.006989139550604536, 0.009258319391439123, 0.004972194549900642, 0.011701012705405973, 0.003796916988492686],
      linf = [0.012897588621946787, 0.03317614952204696, 0.011594778657069232, 0.04437495747712844, 0.025695368105128846, 0.03336443221122043, 0.01000776455308272, 0.044820519469304224, 0.01606741448289395],
      volume_flux = (flux_derigs_etal, flux_nonconservative_powell))
  end

  @trixi_testset "elixir_mhd_alfven_wave_mortar.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave_mortar.jl"),
      l2   = [0.001879021634926363, 0.007032724521848316, 0.0032793932234187325, 0.009056594733320348, 0.007514150120617965, 0.007328739509868727, 0.00309794018112387, 0.009026356949274878, 0.0035732583778049776],
      linf = [0.013734346970999622, 0.06173467158736011, 0.02183946452704291, 0.06258216169457917, 0.03672304497348122, 0.055120532123884625, 0.018202716205672487, 0.06133688282205586, 0.019888161885935608],
      tspan = (0.0, 0.25))
  end

  @trixi_testset "elixir_mhd_alfven_wave.jl with Orszag-Tang setup + flux_hll" begin
    # OBS! This setup does not make much sense and is only used to exercise all components of the
    # flux_hll implementation
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave.jl"),
      l2   = [0.004391143689111404, 0.04144737547475548, 0.041501307637678286, 0.04150353006408862, 0.03693135855995625, 0.021125605214031118, 0.03295607553556973, 0.03296235755245784, 7.16035229384135e-6],
      linf = [0.017894703320895378, 0.08486850681397005, 0.0891044523165206, 0.08492024792056754, 0.10448301878352373, 0.05381260695579509, 0.0884774018719996, 0.07784546966765199, 7.71609149516089e-5],
      initial_condition = initial_condition_orszag_tang,
      surface_flux = (flux_hll, flux_nonconservative_powell),
      volume_flux  = (flux_central, flux_nonconservative_powell),
      coordinates_min = (0, 0, 0),
      coordinates_max = (1, 1, 1),
      initial_refinement_level=3,
      cfl = 1.1,
      tspan = (0.0, 0.06))
  end
end

end # module


text:
module IndirectArrays

export IndirectArray

"""
    IndirectArray(index, values)

creates an array `A` where the values are looked up in the value table,
`values`, using the `index`.  Concretely, `A[i,j] =
values[index[i,j]]`.
"""
struct IndirectArray{T,N,A<:AbstractArray{<:Integer,N},V<:AbstractVector{T}} <: AbstractArray{T,N}
    index::A
    values::V

    @inline function IndirectArray{T,N,A,V}(index, values) where {T,N,A,V}
        # The typical logic for testing bounds and then using
        # @inbounds will not check whether index is inbounds for
        # values. So we had better check this on construction.
        @boundscheck checkbounds(values, index)
        new{T,N,A,V}(index, values)
    end
end
Base.@propagate_inbounds IndirectArray(index::AbstractArray{<:Integer,N}, values::AbstractVector{T}) where {T,N} =
    IndirectArray{T,N,typeof(index),typeof(values)}(index, values)
function IndirectArray{T}(A::AbstractArray) where {T}
    values = unique(A)
    index = convert(Array{T}, indexin(A, values))
    return IndirectArray(index, values)
end
IndirectArray(A::AbstractArray) = IndirectArray{UInt8}(A)

Base.size(A::IndirectArray) = size(A.index)
Base.axes(A::IndirectArray) = axes(A.index)
Base.IndexStyle(::Type{IndirectArray{T,N,A,V}}) where {T,N,A,V} = IndexStyle(A)

Base.copy(A::IndirectArray) = IndirectArray(copy(A.index), copy(A.values))

@inline function Base.getindex(A::IndirectArray, i::Int)
    @boundscheck checkbounds(A.index, i)
    @inbounds idx = A.index[i]
    @boundscheck checkbounds(A.values, idx)
    @inbounds ret = A.values[idx]
    ret
end

@inline function Base.getindex(A::IndirectArray{T,N}, I::Vararg{Int,N}) where {T,N}
    @boundscheck checkbounds(A.index, I...)
    @inbounds idx = A.index[I...]
    @boundscheck checkbounds(A.values, idx)
    @inbounds ret = A.values[idx]
    ret
end

@inline function Base.setindex!(A::IndirectArray, x, i::Int)
    @boundscheck checkbounds(A.index, i)
    idx = findfirst(isequal(x), A.values)
    if idx == nothing
        push!(A.values, x)
        A.index[i] = length(A.values)
    else
        A.index[i] = idx
    end
    return A
end

@inline function Base.push!(A::IndirectArray{T,1} where T, x)
    idx = findfirst(isequal(x), A.values)
    if idx == nothing
        push!(A.values, x)
        push!(A.index, length(A.values))
    else
        push!(A.index, idx)
    end
    return A
end

function Base.append!(A::IndirectArray{T,1}, B::IndirectArray{T,1}) where T
    if A.values == B.values
        append!(A.index, B.index)
    else # pretty inefficient but let's get something going
        for b in B
            push!(A, b)
        end
    end
    return A
end

function Base.append!(A::IndirectArray{<:Any,1}, B::AbstractVector)
    for b in B
        push!(A, b)
    end
    return A
end

end # module


text:
module SkinModHelperSkinSwapTrigger
using ..Ahorn, Maple

@mapdef Trigger "SkinModHelper/SkinSwapTrigger" SkinSwapTrigger(x::Integer, y::Integer, width::Integer=16, height::Integer=16,
	skinId::String="Default", revertOnLeave::Bool=false)

const placements = Ahorn.PlacementDict(
	"Skin Swap Trigger (Skin Mod Helper)" => Ahorn.EntityPlacement(
		SkinSwapTrigger,
		"rectangle"
	)
)

end


text:
function _pretty_print(io::IO, k, color)
    if get(io, :color, false)
        printstyled(io, k; color=color)
    else
        print(io, k)
    end
end

function _get_key(axis, i)
    if i isa StaticInt
        return Int(i)
    else
        return i
    end
end
_get_key(axis::Axis, i) = keys(axis)[findfirst(==(i), eachindex(axis))]
_get_key(axis::StructAxis, i) = keys(axis)[findfirst(==(i), eachindex(axis))]

# this makes it so we don't print StaticInt

function Base.alignment(
    io::IO,
    X::AxisVecOrMat,
    rows::AbstractVector,
    cols::AbstractVector,
    cols_if_complete::Integer,
    cols_otherwise::Integer,
    sep::Integer
)

    a = Tuple{Int, Int}[]
    for j in cols # need to go down each column one at a time
        l, r = alignment(io, _get_key(axes(X, 2), j))
        #=
        l = 0
        r = length(sprint(show, "text/plain", , context=io, sizehint=0))
        =#
        for i in rows # plumb down and see what largest element sizes are
            if isassigned(X,i,j)
                aij = alignment(io, X[i,j])
            else
                aij = undef_ref_alignment
            end
            l = max(l, aij[1]) # left characters
            r = max(r, aij[2]) # right characters
        end
        push!(a, (l, r)) # one tuple per column of X, pruned to screen width
        if length(a) > 1 && sum(map(sum,a)) + sep*length(a) >= cols_if_complete
            pop!(a) # remove this latest tuple if we're already beyond screen width
            break
        end
    end
    if 1 < length(a) < length(axes(X,2))
        while sum(map(sum,a)) + sep * length(a) >= cols_otherwise
            pop!(a)
        end
    end
    return a
end

function print_col_keys(
    io::IO,
    X::AbstractVecOrMat,
    A::Vector,
    cols::AbstractVector,
    sep::AbstractString
)

    for (k, j) = enumerate(cols)
        k > length(A) && break
        x = _get_key(axes(X, 2), j)
        a = alignment(io, x)::Tuple{Int,Int}

        # First try 3-arg show
        sx = sprint(show, "text/plain", x, context=io, sizehint=0)

        # If the output contains line breaks, try 2-arg show instead.
        if occursin('\n', sx)
            sx = sprint(show, x, context=io, sizehint=0)
        end

        _pretty_print(io, sx, 63)
        if j == axes(X, 2)[end]
            print(io, repeat(" ", A[k][1] - a[1]))
        else
            print(io, repeat(" ", (A[k][1]-a[1]) + (A[k][2]-a[2])))
        end

        if k < length(A)
            print(io, sep)
        end
    end
end

function print_row_key(io, axis, i, a, pre, post)
    sx = sprint(show, "text/plain", _get_key(axis, i), context=io, sizehint=0)
    print(io, pre)
    _pretty_print(io, sx, 63)
    print(io, repeat(" ", a + length(post) - length(sx)))
end

function row_alignment(io, row)
    a = 0
    for i in row
        a = max(a, length(sprint(show, "text/plain", _get_key(row, i), context=io, sizehint=0)))
    end
    return a
end

function Base.print_matrix(
    io::IO,
    X::AxisVecOrMat,
    pre::AbstractString = " ",  # pre-matrix string
    sep::AbstractString = "  ", # separator between elements
    post::AbstractString = "  ",  # post-matrix string
    hdots::AbstractString = "  \u2026  ",
    vdots::AbstractString = "\u22ee",
    ddots::AbstractString = "  \u22f1  ",
    hmod::Integer = 5, vmod::Integer = 5
)

    hmod, vmod = Int(hmod)::Int, Int(vmod)::Int
    if !(get(io, :limit, false)::Bool)
        screenheight = screenwidth = typemax(Int)
    else
        sz = displaysize(io)::Tuple{Int,Int}
        screenheight, screenwidth = sz[1] - 4, sz[2]
    end
    screenwidth -= length(pre)::Int + length(post)::Int
    presp = repeat(" ", length(pre)::Int)  # indent each row to match pre string
    postsp = ""
    @assert Base.textwidth(hdots) == Base.textwidth(ddots)
    sepsize = length(sep)::Int
    rowsA = UnitRange{Int}(axes(X,1))
    colsA = UnitRange{Int}(axes(X,2))
    m, n = length(rowsA), length(colsA)
    # To figure out alignments, only need to look at as many rows as could
    # fit down screen. If screen has at least as many rows as A, look at A.
    # If not, then we only need to look at the first and last chunks of A,
    # each half a screen height in size.
    halfheight = div(screenheight, 2)
    if m > screenheight
        rowsA = [rowsA[(0:halfheight - 1) .+ firstindex(rowsA)];
                 rowsA[(end-div(screenheight - 1, 2) + 1):end]]
    end
    # Similarly for columns, only necessary to get alignments for as many
    # columns as could conceivably fit across the screen
    maxpossiblecols = div(screenwidth, 1+sepsize)
    if n > maxpossiblecols
        colsA = [colsA[(0:maxpossiblecols - 1) .+ firstindex(colsA)];
                 colsA[(end - maxpossiblecols + 1):end]]
    end
    A = alignment(io, X, rowsA, colsA, screenwidth, screenwidth, sepsize)
    # Nine-slicing is accomplished using print_matrix_row repeatedly
    if m <= screenheight # rows fit vertically on screen
        if n <= length(A) # rows and cols fit so just print whole matrix in one piece
            rA = row_alignment(io, axes(X, 1))
            print(io, " " * repeat(" ", rA + length(pre) + length(post)))  # print space before column keys
            print_col_keys(io, X, A, colsA, sep)
            println(io)
            for i in rowsA
                print(io, i == first(rowsA) ? pre : presp)
                print_row_key(io, axes(X, 1), i, rA, pre, post)
                print_matrix_row(io, X, A, i, colsA, sep)
                print(io, i == last(rowsA) ? post : postsp)
                if i != last(rowsA)
                    println(io)
                end
            end
        else # rows fit down screen but cols don't, so need horizontal ellipsis
            c = div(screenwidth-length(hdots)::Int+1,2)+1  # what goes to right of ellipsis
            Ralign = reverse(alignment(io, X, rowsA, reverse(colsA), c, c, sepsize)) # alignments for right
            c = screenwidth - sum(map(sum, Ralign)) - (length(Ralign)-1)*sepsize - length(hdots)::Int
            Lalign = alignment(io, X, rowsA, colsA, c, c, sepsize) # alignments for left of ellipsis
            rA = row_alignment(io, axes(X, 1))
            print(io, " " * repeat(" ", rA + length(pre) + length(post)))  # print space before column keys
            print_col_keys(io, X, Lalign, colsA[1:length(Lalign)], sep)
            print(io, repeat(" ", length(hdots)::Int))
            print_col_keys(io, X, Ralign, (n - length(Ralign)) .+ colsA, sep)
            println(io)
            for i in rowsA
                print(io, i == first(rowsA) ? pre : presp)
                print_row_key(io, axes(X, 1), i, rA, pre, post)
                print_matrix_row(io, X, Lalign, i, colsA[1:length(Lalign)], sep)
                print(io, (i - first(rowsA)) % hmod == 0 ? hdots : repeat(" ", length(hdots)::Int))
                print_matrix_row(io, X, Ralign, i, (n - length(Ralign)) .+ colsA, sep)
                print(io, i == last(rowsA) ? post : postsp)
                if i != last(rowsA)
                    println(io)
                end
            end
        end
    else # rows don't fit so will need vertical ellipsis
        if n <= length(A) # rows don't fit, cols do, so only vertical ellipsis
            rA = row_alignment(io, axes(X, 1))
            print(io, " " * repeat(" ", rA + length(pre) + length(post)))  # print space before column keys
            print_col_keys(io, X, A, colsA, sep)
            println(io)
            for i in rowsA
                print_row_key(io, axes(X, 1), i, rA, pre, post)
                print(io, i == first(rowsA) ? pre : presp)
                print_matrix_row(io, X, A, i, colsA, sep)
                print(io, i == last(rowsA) ? post : postsp)
                if i != rowsA[end] || i == rowsA[halfheight]; println(io); end
                if i == rowsA[halfheight]
                    print(io, i == first(rowsA) ? pre : presp)
                    Base.print_matrix_vdots(io, vdots, A, sep, vmod, 1, false)
                    print(io, i == last(rowsA) ? post : postsp * '\n')
                end
            end
        else # neither rows nor cols fit, so use all 3 kinds of dots
            c = div(screenwidth - length(hdots)::Int + 1,2) + 1
            Ralign = reverse(alignment(io, X, rowsA, reverse(colsA), c, c, sepsize))
            c = screenwidth - sum(map(sum,Ralign)) - (length(Ralign)-1)*sepsize - length(hdots)::Int
            Lalign = alignment(io, X, rowsA, colsA, c, c, sepsize)
            r = mod((length(Ralign) - n + 1),vmod) # where to put dots on right half
            rA = row_alignment(io, axes(X, 1))
            print(io, " " * repeat(" ", rA + length(pre) + length(post)))  # print space before column keys
            print_col_keys(io, X, Lalign, colsA[1:length(Lalign)], sep)
            print(io, repeat(" ", length(hdots)::Int))
            print_col_keys(io, X, Ralign, (n - length(Ralign)) .+ colsA, sep)
            println(io)
            for i in rowsA
                print_row_key(io, axes(X, 1), i, rA, pre, post)
                print(io, i == first(rowsA) ? pre : presp)
                print_matrix_row(io, X, Lalign, i, colsA[1:length(Lalign)], sep)
                print(io, (i - first(rowsA)) % hmod == 0 ? hdots : repeat(" ", length(hdots)::Int))
                print_matrix_row(io, X, Ralign, i, (n-length(Ralign)).+colsA, sep)
                print(io, i == last(rowsA) ? post : postsp)
                if i != rowsA[end] || i == rowsA[halfheight]; println(io); end
                if i == rowsA[halfheight]
                    print(io, i == first(rowsA) ? pre : presp)
                    Base.print_matrix_vdots(io, vdots, Lalign, sep, vmod, 1, true)
                    print(io, " " * repeat(" ", rA + length(pre) + length(post)))  # print space before column keys
                    print(io, ddots)
                    Base.print_matrix_vdots(io, vdots, Ralign, sep, vmod, r, false)
                    print(io, i == last(rowsA) ? post : postsp * '\n')
                end
            end
        end
        if isempty(rowsA)
            print(io, pre)
            print(io, vdots)
            length(colsA) > 1 && print(io, "    ", ddots)
            print(io, post)
        end
    end
end

function print_axes_summary(io::IO, axs::NamedTuple{L}) where {L}
    compact_io = IOContext(io, :compact => true)
    lft_pad =lpad(' ', 5)
    print(io, lpad("$(lpad(Char(0x2022), 3)) axes:", 0))
    for i in OneTo(length(axs))
        println(compact_io)
        print(compact_io, lft_pad)
        _pretty_print(compact_io, "$(L[i]) = ", 129)
        _pretty_print(compact_io, axs[i], 63)
    end
end

function print_axes_summary(io::IO, axs::Tuple)
    compact_io = IOContext(io, :compact => true)
    lft_pad =lpad(' ', 5)
    print(io, lpad("$(lpad(Char(0x2022), 3)) axes:", 0))
    for i in OneTo(length(axs))
        println(compact_io)
        print(compact_io, lft_pad)
        _pretty_print(compact_io, "$i = ", 129)
        _pretty_print(compact_io, axs[i], 63)
    end
end

function Base.show_nd(io::IO, a::AxisArray, print_matrix::Function, label_slices::Bool)
    limit::Bool = get(io, :limit, false)
    if isempty(a)
        return
    end
    tailinds = tail(tail(axes(a)))
    nd = ndims(a)-2
    for I in CartesianIndices(tailinds)
        idxs = I.I
        if limit
            for i = 1:nd
                ii = idxs[i]
                ind = tailinds[i]
                if length(ind) > 10
                    if ii == ind[firstindex(ind)+3] && all(d->idxs[d]==first(tailinds[d]),1:i-1)
                        for j=i+1:nd
                            szj = length(axes(a, j+2))
                            indj = tailinds[j]
                            if szj>10 && first(indj)+2 < idxs[j] <= last(indj)-3
                                @goto skip
                            end
                        end
                        #println(io, idxs)
                        print(io, "...\n\n")
                        @goto skip
                    end
                    if ind[firstindex(ind)+2] < ii <= ind[end-3]
                        @goto skip
                    end
                end
            end
        end
        if label_slices
            _pretty_print(io, "[:, :, ", 129)
            for i = 1:(nd-1)
                _pretty_print(io, "$(_get_key(axes(a, i), idxs[i])), ", 63)
            end
            _pretty_print(io, "$(_get_key(axes(a, ndims(a)), idxs[end]))", 63)
            _pretty_print(io, "] =", 129)
            println(io)
        end
        slice = view(a, axes(a,1), axes(a,2), idxs...)
        Base.print_matrix(io, slice)
        print(io, idxs == map(last,tailinds) ? "" : "\n\n")
        @label skip
    end
end

function Base.show(io::IO, ::MIME"text/plain", X::AxisArray)
    if isempty(X) && (get(io, :compact, false) || X isa Vector)
        return show(io, X)
    end
    # 0) show summary before setting :compact
    summary(io, X)
    isempty(X) && return
    Base.show_circular(io, X) && return

    # 1) compute new IOContext
    if !haskey(io, :compact) && length(axes(X, 2)) > 1
        io = IOContext(io, :compact => true)
    end
    if get(io, :limit, false) && eltype(X) === Method
        # override usual show method for Vector{Method}: don't abbreviate long lists
        io = IOContext(io, :limit => false)
    end

    if get(io, :limit, false) && displaysize(io)[1]-4 <= 0
        return print(io, " …")
    else
        println(io)
    end

    # 2) update typeinfo
    #
    # it must come after printing the summary, which can exploit :typeinfo itself
    # (e.g. views)
    # we assume this function is always called from top-level, i.e. that it's not nested
    # within another "show" method; hence we always print the summary, without
    # checking for current :typeinfo (this could be changed in the future)
    io = IOContext(io, :typeinfo => eltype(X))

    # 2) show actual content
    recur_io = IOContext(io, :SHOWN_SET => X)
    Base.print_array(recur_io, X)
end

function Base.showarg(io::IO, x::AxisArray, toplevel)
    if toplevel
        print(io, Base.dims2string(length.(axes(x))), " ")
    end

    print(io, "AxisArray(")
    Base.showarg(io, parent(x), false)
    print(io, "\n")
    print_axes_summary(io, axes(x))
    print(io, "\n)")
end

Base.summary(io::IO, x::AxisArray) = Base.showarg(io, x, true)


text:
module QDLDL

export qdldl, \, solve, solve!

using AMD

const QDLDL_UNKNOWN = -1;
const QDLDL_USED   = true;
const QDLDL_UNUSED = false;

struct QDLDLFactorisation{Tf<:AbstractFloat,Ti<:Integer}

    #contains factors L, D^{-1}
    L::SparseMatrixCSC{Tf,Ti}
    Dinv::Array{Tf}
    #permutation vector (nothing if no permutation)
    p
end

#constructor for enforcing no permutation
function qdldl{Tv<:AbstractFloat,Ti<:Integer}(A::SparseMatrixCSC{Tv,Ti},p::Void)
    #no permutation at all
    (L, Dinv) = factor(triu(A));
    return QDLDLFactorisation(L,Dinv,nothing)
end

#constructor for user-specified permutation, including no permutation
function qdldl{Tv<:AbstractFloat,Ti<:Integer}(A::SparseMatrixCSC{Tv,Ti},p::Array{Ti})
    (L, Dinv) = factor(triu(A[p,p]));
    return QDLDLFactorisation(L,Dinv,p)
end

#constructor for the default (AMD) permutation
function qdldl{Tv<:AbstractFloat,Ti<:Integer}(A::SparseMatrixCSC{Tv,Ti})
    #use AMD ordering as a default
    p = amd(A)::Array{Ti};
    return qdldl(A,p);
end

function Base.:\(QDLDL::QDLDLFactorisation,b)
    return solve(QDLDL,b)
end

# Solves Ax = b using LDL factors for A.
# Returns x, preserving b
function solve(QDLDL::QDLDLFactorisation,b)
    x = copy(b)
    solve!(QDLDL,x)
    return x
end

# Solves Ax = b using LDL factors for A.
# Solves in place (x replaces b)
function solve!(QDLDL::QDLDLFactorisation,b)
    #permute b
    if QDLDL.p != nothing
        permute!(b,QDLDL.p)
    end
    QDLDL_solve!(QDLDL.L.n,
                 QDLDL.L.colptr,
                 QDLDL.L.rowval,
                 QDLDL.L.nzval,
                 QDLDL.Dinv,b)
    #inverse permutation
    if QDLDL.p != nothing
        ipermute!(b,QDLDL.p)
    end
    return nothing
end

function factor{Tv<:AbstractFloat,Ti<:Integer}(A::SparseMatrixCSC{Tv,Ti})

    A = triu(A)

    etree  = Array{Ti}(A.n)
    Lnz    = Array{Ti}(A.n)
    iwork  = Array{Ti}(A.n*3)
    bwork  = Array{Bool}(A.n)
    fwork  = Array{Tv}(A.n)

    #compute elimination gree using QDLDL converted code
    sumLnz = QDLDL_etree!(A.n,A.colptr,A.rowval,iwork,Lnz,etree)

    if(sumLnz < 0)
        error("A matrix is not upper triangular or has an empty column")
    end

    #allocate space for the L matrix row indices and data
    Lp = Array{Ti}(A.n + 1);
    Li = Array{Ti}(sumLnz);
    Lx = Array{Tv}(sumLnz);

    #allocate for and D D inverse
    D  = Array{Tv}(A.n);
    Dinv = Array{Tv}(A.n);

    #factor using QDLDL converted code
    posDCount = QDLDL_factor(A.n,A.colptr,A.rowval,A.nzval,Lp,Li,Lx,D,Dinv,Lnz,etree,bwork,iwork,fwork)

    if(posDCount < 0)
        error("Zero entry in D (matrix is not quasidefinite)")
    end

    L = SparseMatrixCSC(A.n,A.n,Lp,Li,Lx);

    return L, Dinv

end

# Compute the elimination tree for a quasidefinite matrix
# in compressed sparse column form.

function QDLDL_etree!(n,Ap,Ai,work,Lnz,etree)

    @inbounds for i = 1:n
      # zero out Lnz and work.  Set all etree values to unknown
        work[i]  = 0
        Lnz[i]   = 0i
        etree[i] = QDLDL_UNKNOWN

        #Abort if A doesn't have at least one entry
        #one entry in every column
        if(Ap[i] == Ap[i+1])
          return -1
        end
    end

    @inbounds for j = 1:n
        work[j] = j
        @inbounds for p = Ap[j]:(Ap[j+1]-1)
          i = Ai[p]
            if(i > j)
                return -1
            end
            @inbounds while(work[i] != j)
                if(etree[i] == QDLDL_UNKNOWN)
                    etree[i] = j
                end
                Lnz[i] += 1        #nonzeros in this column
                work[i] = j
                i = etree[i]
            end
        end #end for p
    end

  #tally the total nonzeros
  sumLnz = sum(Lnz);

  return sumLnz
end

function QDLDL_factor(n,Ap,Ai,Ax,Lp,Li,Lx,D,Dinv,Lnz,etree,bwork,iwork,fwork)

  positiveValuesInD = 0

  #partition working memory into pieces
  yMarkers        = bwork
  yIdx            = view(iwork,      1:n)
  elimBuffer      = view(iwork,  (n+1):2*n)
  LNextSpaceInCol = view(iwork,(2*n+1):3*n)
  yVals           = fwork;

  Lp[1] = 1 #first column starts at index one / Julia is 1 indexed

  @inbounds for i = 1:n

    #compute L column indices
    Lp[i+1] = Lp[i] + Lnz[i]   #cumsum, total at the end

    # set all Yidx to be 'unused' initially
    #in each column of L, the next available space
    #to start is just the first space in the column
    yMarkers[i]  = QDLDL_UNUSED
    yVals[i]     = 0.0
    D[i]         = 0.0
    LNextSpaceInCol[i] = Lp[i]
end

  # First element of the diagonal D.
  D[1]     = Ax[1]
  if(D[1] == 0.0) return -1 end
  if(D[1]  > 0.0) positiveValuesInD += 1 end
  Dinv[1] = 1/D[1];

  #Start from 1 here. The upper LH corner is trivially 0
  #in L b/c we are only computing the subdiagonal elements
  @inbounds for k = 2:n

    #NB : For each k, we compute a solution to
    #y = L(0:(k-1),0:k-1))\b, where b is the kth
    #column of A that sits above the diagonal.
    #The solution y is then the kth row of L,
    #with an implied '1' at the diagonal entry.

    #number of nonzeros in this row of L
    nnzY = 0  #number of elements in this row

    #This loop determines where nonzeros
    #will go in the kth row of L, but doesn't
    #compute the actual values
    @inbounds for i = Ap[k]:(Ap[k+1]-1)

      bidx = Ai[i]   # we are working on this element of b

      #Initialize D[k] as the element of this column
      #corresponding to the diagonal place.  Don't use
      #this element as part of the elimination step
      #that computes the k^th row of L
      if(bidx == k)
        D[k] = Ax[i];
        continue
       end

      yVals[bidx] = Ax[i]   # initialise y(bidx) = b(bidx)

      # use the forward elimination tree to figure
      # out which elements must be eliminated after
      # this element of b
      nextIdx = bidx

      if(yMarkers[nextIdx] == QDLDL_UNUSED)  #this y term not already visited

        yMarkers[nextIdx] = QDLDL_USED     #I touched this one
        elimBuffer[1]     = nextIdx  # It goes at the start of the current list
        nnzE              = 1         #length of unvisited elimination path from here

        nextIdx = etree[bidx];

        @inbounds while(nextIdx != QDLDL_UNKNOWN && nextIdx < k)
          if(yMarkers[nextIdx] == QDLDL_USED) break; end

          yMarkers[nextIdx] = QDLDL_USED;   #I touched this one
          #NB: Julia is 1-indexed, so I increment nnzE first here,
          #no after writing into elimBuffer as in the C version
          nnzE += 1                   #the list is one longer than before
          elimBuffer[nnzE] = nextIdx; #It goes in the current list
          nextIdx = etree[nextIdx];   #one step further along tree

        end #end while

        # now I put the buffered elimination list into
        # my current ordering in reverse order
        @inbounds while(nnzE != 0)
         #NB: inc/dec reordered relative to C because
         #the arrays are 1 indexed
         nnzY += 1;
         yIdx[nnzY] = elimBuffer[nnzE];
         nnzE -= 1;
       end #end while
   end #end if

end #end for i

    #This for loop places nonzeros values in the k^th row
    @inbounds for i = nnzY:-1:1

      #which column are we working on?
      cidx = yIdx[i]

      # loop along the elements in this
      # column of L and subtract to solve to y
      tmpIdx = LNextSpaceInCol[cidx];
      yVals_cidx = yVals[cidx]
      @inbounds for j = Lp[cidx]:(tmpIdx-1)
        yVals[Li[j]] -= Lx[j]*yVals_cidx
      end

      #Now I have the cidx^th element of y = L\b.
      #so compute the corresponding element of
      #this row of L and put it into the right place
      Li[tmpIdx] = k
      Lx[tmpIdx] = yVals_cidx *Dinv[cidx]

      #D[k] -= yVals[cidx]*yVals[cidx]*Dinv[cidx];
      D[k] -= yVals_cidx*Lx[tmpIdx]
      LNextSpaceInCol[cidx] += 1

      #reset the yvalues and indices back to zero and QDLDL_UNUSED
      #once I'm done with them
      yVals[cidx]     = 0.0
      yMarkers[cidx]  = QDLDL_UNUSED;

    end #end for i

    #Maintain a count of the positive entries
    #in D.  If we hit a zero, we can't factor
    #this matrix, so abort
    if(D[k] == 0.0) return -1 end
    if(D[k]  > 0.0) positiveValuesInD += 1 end

    #compute the inverse of the diagonal
    Dinv[k]= 1/D[k]

end #end for k

  return positiveValuesInD

end

# Solves (L+I)x = b, with x replacing b
function QDLDL_Lsolve!(n,Lp,Li,Lx,x)

  @inbounds for i = 1:n
      @inbounds for j = Lp[i]: (Lp[i+1]-1)
          x[Li[j]] -= Lx[j]*x[i];
      end
  end
  return nothing
end

# Solves (L+I)'x = b, with x replacing b
function QDLDL_Ltsolve!(n,Lp,Li,Lx,x)

  @inbounds for i = n:-1:1
      @inbounds for j = Lp[i]:(Lp[i+1]-1)
          x[i] -= Lx[j]*x[Li[j]]
      end
  end
  return nothing
end

# Solves Ax = b where A has given LDL factors,
# with x replacing b
function QDLDL_solve!(n,Lp,Li,Lx,Dinv,b)

    QDLDL_Lsolve!(n,Lp,Li,Lx,b)
    b .*= Dinv;
    QDLDL_Ltsolve!(n,Lp,Li,Lx,b)

end

end #end module

max_stars_repo_name: ANYbotics/qdldl

max_issues_repo_name: ANYbotics/qdldl

max_forks_repo_name: ANYbotics/qdldl


text:
using NewPackage
using Test

include("foo_test.jl")


text:
@testset "Candidate" begin

    struct DummyFitness <: AbstractFitness end
    NaiveGAflux._fitness(::DummyFitness, f::AbstractCandidate) = 17
    using NaiveGAflux: FileCandidate, AbstractWrappingCandidate, FittedCandidate
    using Functors: fmap
    import MemPool
    @testset "$ctype" for (ctype, candfun) in (
        (CandidateModel, CandidateModel),
        (CandidateOptModel, g -> CandidateOptModel(Descent(0.01), g))
    )

        @testset " $lbl" for (lbl, wrp) in (
            ("", identity),
            (FileCandidate, FileCandidate),
            (FittedCandidate, c -> FittedCandidate(3, 0.34, c)),
            ("FittedCandidate ∘ FileCandidate", c -> FittedCandidate(1, 0.78, FileCandidate(c)))
            )

            invertex = denseinputvertex("in", 3)
            hlayer = fluxvertex("hlayer", Dense(3,4), invertex)
            outlayer = fluxvertex("outlayer", Dense(4, 2), hlayer)
            graph = CompGraph(invertex, outlayer)

            try

                cand = wrp(candfun(graph))

                @test fitness(DummyFitness(), cand) == 17
                # Just to make sure we can access the graph
                @test fitness(SizeFitness(), cand) == nparams(graph)

                graphmutation = VertexMutation(MutationFilter(v -> name(v)=="hlayer", AddVertexMutation(ArchSpace(DenseSpace([1], [relu])))))
                optmutation = OptimizerMutation((Momentum, Nesterov, ADAM))
                evofun = evolvemodel(graphmutation, optmutation)
                newcand = evofun(cand)

                @test NaiveGAflux.model(nvertices, newcand) == 4
                @test NaiveGAflux.model(nvertices, cand) == 3

                optimizer(c) = typeof(opt(c))

                if ctype === CandidateOptModel
                    @test optimizer(newcand) !== optimizer(cand) !== Nothing
                    fmapped = fmap(identity, newcand)
                    @test opt(fmapped) !== opt(newcand)
                    @test optimizer(fmapped) === optimizer(newcand)
                else
                    @test optimizer(newcand) === optimizer(cand) === Nothing
                end

                teststrat() = NaiveGAflux.default_crossoverswap_strategy(v -> 1)
                graphcrossover = VertexCrossover(CrossoverSwap(;pairgen = (v1,v2) -> (1,1), strategy=teststrat); pairgen = (v1,v2;ind1) -> ind1==1 ? (2,3) : nothing)
                optcrossover = OptimizerCrossover()
                crossfun = evolvemodel(graphcrossover, optcrossover)

                newcand1, newcand2 = crossfun((cand, newcand))

                @test optimizer(newcand1) === optimizer(newcand)
                @test optimizer(newcand2) === optimizer(cand)

                @test NaiveGAflux.model(nvertices, newcand1) == 4
                @test NaiveGAflux.model(nvertices, newcand2) == 3
            finally
                MemPool.cleanup()
            end
        end
    end

    @testset "FileCandidate" begin
        try
           @testset "FileCandidate cleanup" begin
                fc = FileCandidate([1,2,3], 1.0)
                # Waiting for timer to end does not seem to be reliable, so we'll just stop the timer and call the timeout function manually
                close(fc.movetimer)
                NaiveGAflux.candtodisk(fc.c, fc.writelock)

                fname = MemPool.default_path(fc.c)

                @test isfile(fname)

                finalize(fc)
                @test !isfile(fname)
            end

            @testset "Functor" begin
                v1 = fluxvertex("v1", Dense(3,3;init=Flux.identity_init), denseinputvertex("in", 3))
                cand1 = FileCandidate(CandidateModel(CompGraph(inputs(v1)[], v1)))

                mul(x::AbstractArray) = 2 .* x
                mul(x) = x
                cand2 = Flux.fmap(mul, cand1)
                @test cand2 isa FileCandidate

                @test cand1.hold == cand2.hold == false

                indata = collect(Float32, reshape(1:6,3,2))
                @test model(cand2)(indata) == 2 .* indata
                @test cand2.hold == true

                # We accessed wrappedcand when calling graph, so now we hold the model in memory until release!:d
                # Make sure this is the case for fmapped candidates too
                cand3 = Flux.fmap(identity, cand2)
                @test cand3.hold == cand2.hold == true
            end

            @testset "Serialization" begin
                using Serialization

                v = fluxvertex("v1", Dense(3,3), denseinputvertex("in", 3))
                secand = FileCandidate(CandidateModel(CompGraph(inputs(v)[], v)))

                indata = randn(3, 2)
                expected = v(indata)

                io = PipeBuffer()
                serialize(io, secand)

                MemPool.cleanup()

                # Make sure the data is gone
                @test_throws KeyError NaiveGAflux.model(secand)

                decand = deserialize(io)

                @test model(decand)(indata) == expected
            end

            @testset "Hold in mem" begin
                import NaiveGAflux: wrappedcand, callcand, candinmem
                struct BoolCand <: AbstractCandidate
                    x::Ref{Bool}
                end
                testref(c::BoolCand, f=identity) = f(c.x)
                testref(c::AbstractWrappingCandidate) = testref(wrappedcand(c))
                testref(c::FileCandidate, f) =  callcand(testref, c, f)

                fc = FileCandidate(BoolCand(Ref(true)), 0.1)

                #1 test the testcase: Change true to false after move to disk
                x = testref(fc, identity) # And neither should this now
                @test isopen(fc.movetimer) == true

                t0 = time()
                while candinmem(fc) && time() - t0 < 5
                    sleep(0.01)
                end
                @test time() - t0 < 5 # Or else we timed out
                x[] = false

                @test false == x[] != testref(fc)[] # Note model(fc) should not start movetimer
                @test !isopen(fc.movetimer)

                x = testref(fc, identity) # And neither should this now
                @test !isopen(fc.movetimer)
                NaiveGAflux.release!(fc)
                @test isopen(fc.movetimer)
            end

            @testset "Move to disk collision" begin
                using NaiveGAflux: AbstractCandidate, wrappedcand, release!
                using Serialization

                struct TakesLongToSerialize <: AbstractCandidate end

                finish_serialize = Condition()
                function Serialization.serialize(s::AbstractSerializer, ::TakesLongToSerialize)
                    wait(finish_serialize)
                    Serialization.writetag(s.io, Serialization.OBJECT_TAG)
                    serialize(s, TakesLongToSerialize)
                end

                fc = FileCandidate(TakesLongToSerialize(), 10.0)
                # Waiting for timer to end does not seem to be reliable, so we'll just stop the timer and call the timeout function manually
                close(fc.movetimer)
                # Async as the serialization will block
                movetask = @async NaiveGAflux.candtodisk(fc.c, fc.writelock)

                # Also blocked due to locking
                gettask = @async @test_logs (:warn, r"Try to access FileCandidate which is being moved to disk") wrappedcand(fc)

                # Hammer the condition to improve robustness towards deadlocks
                cnt = 0
                while !istaskdone(gettask)
                    sleep(0.05)
                    notify(finish_serialize)
                    cnt += 1
                    if cnt == 20
                        @warn "Task did not complete after $cnt attempts. Testcase might be deadlocked!"
                    end
                end

                @test fetch(gettask) == TakesLongToSerialize()
                # Calling wrappedcand causes FileCandidate to hold the candidate in memory until released
                @test isopen(fc.movetimer) == false

                release!(fc)

                @test isopen(fc.movetimer) == true
                close(fc.movetimer) # So we don't create garbage after fc.movedelay seconds
            end

        finally
            MemPool.cleanup()
        end
    end

    @testset "Global optimizer mutation" begin
        import NaiveGAflux.Flux.Optimise: Optimiser
        import NaiveGAflux: sameopt, learningrate, BoundedRandomWalk, global_optimizer_mutation, randomlrscale

        @testset "Random learning rate scale" begin
            using Random
            so = ShieldedOpt

            omf = randomlrscale();

            om1 = omf()
            @test learningrate(om1(Descent(0.1))) ≈ learningrate(om1(Momentum(0.1)))

            opt = Optimiser(so(Descent(0.1)), Momentum(0.1), so(Descent(1.0)), ADAM(1.0), Descent(1.0))
            @test length(om1(opt).os) == 4
            @test learningrate(om1(opt)) ≈ learningrate(om1(Descent(0.01)))

            om2 = omf()
            @test learningrate(om2(Descent(0.1))) ≈ learningrate(om2(Momentum(0.1)))

            # Differnt iterations shall yield different results ofc
            @test learningrate(om1(Descent(0.1))) != learningrate(om2(Momentum(0.1)))

            # Make sure learning rate stays within bounds when using BoundedRandomWalk
            rng = MersenneTwister(0);
            brw = BoundedRandomWalk(-1.0, 1.0, () -> randn(rng))
            @test collect(extrema(cumprod([10^brw() for i in 1:10000]))) ≈ [0.1, 10] atol = 1e-10
        end

        @testset "Global learning rate scaling" begin
            v1 = denseinputvertex("in", 3)
            pop = CandidateOptModel.(Descent.(0.1:0.1:1.0), Ref(CompGraph(v1, v1)))

            lr(c) = c.opt.eta
            @test lr.(pop) == 0.1:0.1:1.0

            popscal = global_optimizer_mutation(pop, pp -> OptimizerMutation(o -> sameopt(o, 10learningrate(o))))

            @test lr.(popscal) == 1:10
        end
    end
end


text:
@testset "Densities" begin
    @testset "Single functions" begin
        rmin = 0.01
        rmax = 10.0
        δr = 0.05
        N = ceil(Int, rmax/δr)
        k = 7
        Nn = 71

        δrmin=rmin
        δrmax=0.5
        α=0.01

        @testset "$R" for (R,kind,rtol) in [
            (FiniteDifferences(N, δr), :orthogonal_uniform, 1e-14),
            (StaggeredFiniteDifferences(δrmin, δrmax, α, rmax), :orthogonal_non_uniform, 1e-14),
            (FEDVR(range(0, stop=rmax, length=Nn), k), :orthogonal_non_uniform, 1e-14),
            (BSpline(LinearKnotSet(k, 0, rmax, Nn)), :non_orthogonal, 1e-5),
        ]
            f = x -> sin(2π*x)
            g = x -> x*exp(-x)
            h = x -> sin(2π*x)*x*exp(-x)

            r = axes(R,1)

            cf = R \ f.(r)
            cg = R \ g.(r)
            ch = R \ h.(r)

            ρ = Density(applied(*,R,cf), applied(*,R,cg))
            ch2 = ρ.ρ

            @test eltype(ρ) == eltype(cf)

            if kind == :orthogonal_uniform
                @test ρ.LV == I
                @test ρ.RV == I
                @test ρ.C == I
            elseif kind == :orthogonal_non_uniform
                @test ρ.LV == I
                @test ρ.RV == I
            end

            @test ch2 ≈ ch atol=1e-14 rtol=rtol
        end
    end

    @testset "Multiple functions" begin
        L = 6.0
        N = 1000
        ρ = L/(N+1)
        @testset "$R" for R in [
            FiniteDifferences(N, ρ),
            StaggeredFiniteDifferences(0.01, 0.2, 0.1, L),
            FEDVR(range(0, stop=L, length=20), 7)[:,2:end-1],
            BSpline(LinearKnotSet(7, 0, L, 200))
        ]
            r = axes(R,1)
            f = exp.(-(r .- L/2).^2/0.5)
            gs = [
                sin.(1π/L*r),
                sin.(2π/L*r),
                sin.(3π/L*r),
                sin.(4π/L*r),
                sin.(5π/L*r),
                sin.(6π/L*r)
            ]
            w = r -> 1 / r

            cf = R \ f
            cgs = reduce(hcat, [R \ g for g in gs])
            wr = w.(r)

            u = applied(*, R, cf)
            v = applied(*, R, cgs)

            # One times many functions
            ρuv = Density(u, v, w=w)
            uv_refs = reduce(hcat, [R \ (f .* g .* wr) for g in gs])
            @test ρuv.ρ ≈ uv_refs

            # Many times many functions
            ρvv = Density(v, v, w=w)
            vv_refs = reduce(hcat, [R \ (g .* g .* wr) for g in gs])
            @test ρvv.ρ ≈ vv_refs

            @test_throws DimensionMismatch Density(applied(*, R, rand(size(R, 2), 2)), applied(*, R, rand(size(R, 2), 3)))
        end
    end

    @testset "Restricted bases" begin
        rmin = 0.0
        rmax = 10.0
        δr = 0.05
        N = ceil(Int, rmax/δr)
        k = 7
        Nn = 71

        δrmin=0.01
        δrmax=0.5
        α=0.01

        f = x -> sinpi(x/rmax)
        g = x -> cospi(x/rmax)
        @testset "$(label)" for (w,label) in [(one,"Unweighted"),(exp,"Weighted")]
            h = x -> f(x)*w(x)*g(x)

            @testset "$R" for (R,rtol) in [
                (FiniteDifferences(N, δr), 1e-15),
                (StaggeredFiniteDifferences(δrmin, δrmax, α, rmax), 1e-15),
                (FEDVR(range(0, stop=rmax, length=Nn), k), 1e-15),
                (BSpline(LinearKnotSet(k, 0, rmax, Nn)), 1e-14),
            ]
                R̃ = R[:,1:ceil(Int, size(R,2))]

                r = axes(R̃,1)

                u = R̃ \ f.(r)
                v = R̃ \ g.(r)

                ρ = Density(applied(*, R̃, u), applied(*, R̃, v), w=w)
                @test ρ.ρ ≈ (R̃ \ h.(r)) rtol=rtol
            end
        end
    end

    @testset "Pretty-printing" begin
        R = FiniteDifferences(99, 0.1)
        u = R*rand(size(R,2))
        v = R*rand(size(R,2))
        ρ = Density(u, v)
        fp = FunctionProduct{false}(u, v)

        @test string(ρ) == "99 .* 99 -> 99 FunctionProduct Float64, conjugated (<=> Density)"
        @test string(fp) == "99 .* 99 -> 99 FunctionProduct Float64"

        buf = IOBuffer()
        show(buf, MIME"text/plain"(), ρ)
        @test String(take!(buf)) == """
            99 .* 99 -> 99 FunctionProduct Float64, conjugated (<=> Density); L .* R -> R, with
              L: $(string(R))
              R: $(string(R))"""

        buf = IOBuffer()
        show(buf, MIME"text/plain"(), fp)
        @test String(take!(buf)) == """
            99 .* 99 -> 99 FunctionProduct Float64; L .* R -> R, with
              L: $(string(R))
              R: $(string(R))"""

        u = R*rand(size(R,2), 1)
        v = R*rand(size(R,2), 5)
        ρ = Density(u, v)
        @test string(ρ) == "(99, :) .* (99, :) -> (99, 5) FunctionProduct Float64, conjugated (<=> Density)"

        u = R*rand(size(R,2), 5)
        v = R*rand(size(R,2), 5)
        ρ = Density(u, v)
        @test string(ρ) == "(99, :) .* (99, :) -> (99, 5) FunctionProduct Float64, conjugated (<=> Density)"

        R = BSpline(LinearKnotSet(7, 0, 6.0, 200))

        u = R*rand(size(R,2))
        v = R*rand(size(R,2), 5)
        ρ = Density(u, v)
        @test string(ρ) == "206 .* (206, 5) -> (206, 5) FunctionProduct Float64, conjugated (<=> Density)"

        u = R*rand(size(R,2), 5)
        v = R*rand(size(R,2))
        ρ = Density(u, v)
        @test string(ρ) == "(206, 5) .* 206 -> (206, 5) FunctionProduct Float64, conjugated (<=> Density)"

        u = R*rand(size(R,2), 1)
        v = R*rand(size(R,2), 5)
        ρ = Density(u, v)
        @test string(ρ) == "(206, 1) .* (206, 5) -> (206, 5) FunctionProduct Float64, conjugated (<=> Density)"

        u = R*rand(size(R,2), 5)
        v = R*rand(size(R,2), 1)
        ρ = Density(u, v)
        @test string(ρ) == "(206, 5) .* (206, 1) -> (206, 5) FunctionProduct Float64, conjugated (<=> Density)"

        u = R*rand(size(R,2), 5)
        v = R*rand(size(R,2), 5)
        ρ = Density(u, v)
        @test string(ρ) == "(206, 5) .* (206, 5) -> (206, 5) FunctionProduct Float64, conjugated (<=> Density)"
    end
end


text:
using ProgressMeter

function main()
    out_dir="result"
    PARAMDAT="PARAMDAT.json"

    xmax, ymax, zmax, nodes, vecAx, vecAy, vecAz = read_allgrid()

    out_file_front, out_ext, restartnum, restart_file, init_small, norm_ok,
    time_integ, nt, dt, every_outnum, in_nt, dtau,
    init_rho, init_u, init_v, init_w, init_p, specific_heat_ratio, Rd, bdcon = input_para(PARAMDAT)

    Qbase, cellxmax, cellymax, cellzmax, restartnum = set_initQbase(xmax, ymax, zmax, restart_file, init_rho,
                                                                    init_u, init_v, init_w, init_p,
                                                                    specific_heat_ratio, out_file_front,
                                                                    out_ext, out_dir, restartnum, Rd)

    # init Delta_Qcon_hat
    volume         = set_volume(nodes, vecAx, vecAy, vecAz, cellxmax, cellymax, cellzmax)
    cellcenter     = set_cellcenter(nodes, cellxmax, cellymax, cellzmax)

    wallpoint, nop = set_wallpoint(nodes, bdcon, cellxmax, cellymax, cellzmax)
    wally          = set_wally(cellcenter, wallpoint, nop, cellxmax, cellymax, cellzmax)

    # main loop
    print("threads num : ")
    println(Threads.nthreads())

    #throw(UndefVarError(:x))
    # main loop
    prog = Progress(nt,1)
    for t in 1:nt
        next!(prog)

        evalnum = t+restartnum
        if time_integ == "1"
            # exlicit scheme
            Qbase    = set_boundary(Qbase, cellxmax, cellymax, cellzmax, vecAx, vecAy, vecAz, bdcon)
            Qcon     = base_to_conservative(Qbase, cellxmax, cellymax, cellzmax, specific_heat_ratio)
            Qcon_hat = setup_Qcon_hat(Qcon, cellxmax, cellymax, cellzmax, volume)

            # initial_setup
            mu     = set_mu(Qbase, cellxmax, cellymax, cellzmax, specific_heat_ratio, Rd)
            lambda = set_lambda(Qbase, cellxmax, cellymax, cellzmax, mu, specific_heat_ratio, Rd)

            # RHS
            # advection_term
            # E_adv_hat, F_adv_hat = AUSM(Qbase,Qcon,cellxmax, cellymax, cellzmax,vecAx,vecAy,specific_heat_ratio)
            cell_Ahat_plas, cell_Ahat_minus, cell_Bhat_plas, cell_Bhat_minus, cell_Chat_plas, cell_Chat_minus =
                cal_jacobi(Qbase, Qcon, cellxmax, cellymax, cellzmax, specific_heat_ratio, vecAx, vecAy, vecAz, volume)
            cell_E_hat_plas, cell_E_hat_minus, cell_F_hat_plas, cell_F_hat_minus, cell_G_hat_plas, cell_G_hat_minus =
                setup_cell_flux_hat(Qcon, cellxmax, cellymax, cellzmax, cell_Ahat_plas, cell_Ahat_minus, cell_Bhat_plas, cell_Bhat_minus, cell_Chat_plas, cell_Chat_minus)
            E_adv_hat, F_adv_hat, G_adv_hat =
                FVS(cell_E_hat_plas, cell_E_hat_minus, cell_F_hat_plas, cell_F_hat_minus, cell_G_hat_plas, cell_G_hat_minus, cellxmax, cellymax, cellzmax)

            # viscos_term
            E_vis_hat, F_vis_hat, G_vis_hat = central_diff(Qbase, Qcon, cellxmax, cellymax, cellzmax, mu, lambda, vecAx, vecAy, vecAz, specific_heat_ratio, volume, Rd)
            #=
            E_vis_hat = zeros(cellxmax+1, cellymax+1, cellzmax+1, 5)
            F_vis_hat = zeros(cellxmax+1, cellymax+1, cellzmax+1, 5)
            G_vis_hat = zeros(cellxmax+1, cellymax+1, cellzmax+1, 5)
            =#

            # turblence_term
            yplus = cal_yplus(Qbase, wally, mu, cellxmax, cellymax, cellzmax)

            # RHS
            RHS = setup_RHS(cellxmax, cellymax, cellzmax, E_adv_hat, F_adv_hat, G_adv_hat, E_vis_hat, F_vis_hat, G_vis_hat)

            # time integral
            Qcon_hat = time_integration_explicit(dt, Qcon_hat, RHS, cellxmax, cellymax, cellzmax)

            # reset
            Qcon  = Qhat_to_Q(Qcon_hat, cellxmax, cellymax, cellzmax, volume)
            Qbase = conservative_to_base(Qcon, cellxmax, cellymax, cellzmax, specific_heat_ratio)

        elseif time_integ == "2"
            # Qcn
            Qbasen    = copy(Qbase)
            Qconn     = base_to_conservative(Qbasen, cellxmax, cellymax, cellzmax, specific_heat_ratio)
            Qconn_hat = setup_Qcon_hat(Qconn, cellxmax, cellymax, cellzmax, volume)

            Qbasem = copy(Qbase)

            for tau in 1:in_nt
                # LHS (A_adv_hat=jacobian)
                Qbasem   = set_boundary(Qbasem, cellxmax, cellymax, cellzmax, vecAx, vecAy, vecAz, bdcon)
                Qcon     = base_to_conservative(Qbasem, cellxmax, cellymax, cellzmax, specific_heat_ratio)
                Qcon_hat = setup_Qcon_hat(Qcon, cellxmax, cellymax, cellzmax, volume)

                # initial_setup
                mu     = set_mu(Qbasem, cellxmax, cellymax, cellzmax, specific_heat_ratio, Rd)
                lambda = set_lambda(Qbasem, cellxmax, cellymax, cellzmax, mu, specific_heat_ratio, Rd)

                # turb
                # tauwall = set_tauwall(Qbasem,wally,mu,cellxmax, cellymax, cellzmax)

                # RHS
                #advection_term
                #E_adv_hat, F_adv_hat = AUSM(Qbasem,Qcon,cellxmax, cellymax, cellzmax,vecAx,vecAy,specific_heat_ratio)
                cell_Ahat_plas, cell_Ahat_minus, cell_Bhat_plas, cell_Bhat_minus, cell_Chat_plas, cell_Chat_minus =
                    cal_jacobi(Qbasem, Qcon, cellxmax, cellymax, cellzmax, specific_heat_ratio, vecAx, vecAy, vecAz, volume)
                cell_E_hat_plas, cell_E_hat_minus, cell_F_hat_plas, cell_F_hat_minus, cell_G_hat_plas, cell_G_hat_minus =
                    setup_cell_flux_hat(Qcon, cellxmax, cellymax, cellzmax, cell_Ahat_plas, cell_Ahat_minus, cell_Bhat_plas, cell_Bhat_minus, cell_Chat_plas, cell_Chat_minus)
                E_adv_hat, F_adv_hat, G_adv_hat =
                    FVS(cell_E_hat_plas, cell_E_hat_minus, cell_F_hat_plas, cell_F_hat_minus, cell_G_hat_plas, cell_G_hat_minus, cellxmax, cellymax, cellzmax)

                # viscos_term
                E_vis_hat, F_vis_hat, G_vis_hat = central_diff(Qbasem, Qcon, cellxmax, cellymax, cellzmax, mu, lambda, vecAx, vecAy, vecAz, specific_heat_ratio, volume, Rd)

                RHS = setup_RHS(cellxmax, cellymax, cellzmax, E_adv_hat, F_adv_hat, G_adv_hat, E_vis_hat, F_vis_hat, G_vis_hat)

                #println(E_adv_hat[:,50,:])
                #println(E_adv_hat[:,50,:])

                #println(RHS[:,50,:])

                #throw(UndefVarError(:x))

                # lusgs_advection_term
                A_adv_hat_p, A_adv_hat_m, B_adv_hat_p, B_adv_hat_m, C_adv_hat_p, C_adv_hat_m, A_beta_shig, B_beta_shig, C_beta_shig =
                    one_wave(Qbase, Qcon, cellxmax, cellymax, cellzmax, vecAx, vecAy, vecAz, specific_heat_ratio, volume)
                # lusgs_viscos_term
                jalphaP, jbetaP, jgammaP =
                    central_diff_jacobian(Qbasem, Qcon, cellxmax, cellymax, cellzmax, mu, lambda, vecAx, vecAy, vecAz, specific_heat_ratio, volume)
                #jalphaP = zeros(cellxmax, cellymax)
                #jbetaP = zeros(cellxmax, cellymax)

                # LUSGS
                delta_Q = zeros(cellxmax, cellymax, cellzmax,5)

                ite=0
                while true
                    delta_Q_temp = copy(delta_Q)
                    delta_Q = lusgs(dt, Qcon_hat, Qconn_hat, delta_Q, A_adv_hat_p, A_adv_hat_m, B_adv_hat_p, B_adv_hat_m, C_adv_hat_p, C_adv_hat_m,
                                    A_beta_shig, B_beta_shig, C_beta_shig, jalphaP, jbetaP, jgammaP, RHS, cellxmax, cellymax,cellzmax, volume)

                    #println(RHS[2,:,:])
                    #println(delta_Q[2,:,:])

                    res   = set_res(delta_Q, delta_Q_temp, cellxmax, cellymax, cellzmax)
                    norm2 = check_converge(res, RHS, cellxmax, cellymax,cellzmax, init_small)

                    ite += 1
                    if norm2[1] < norm_ok && norm2[2] < norm_ok && norm2[3] < norm_ok && norm2[4] < norm_ok && norm2[5] < norm_ok
                        break
                    end
                    if ite % 100 ==0
                        println(norm2)
                    end
                end
                #println(ite)
                #println(delta_Q[:,50,:])

                for i in 2:cellxmax-1
                    for j in 2:cellymax-1
                        for k in 2:cellzmax-1
                            for l in 1:5
                                Qcon_hat[i,j,k,l] = Qcon_hat[i,j,k,l] + delta_Q[i,j,k,l]
                            end
                        end
                    end
                end

                Qcon   = Qhat_to_Q(Qcon_hat, cellxmax, cellymax, cellzmax, volume)
                Qbasem = conservative_to_base(Qcon, cellxmax, cellymax, cellzmax, specific_heat_ratio)

                #println("xxx")
            end
            Qbase = copy(Qbasem)
        end

        if round(evalnum) % every_outnum == 0
            println("nt_______________________________"*string(round(evalnum)))
            output_result(evalnum, Qbase, cellxmax, cellymax, cellzmax, out_file_front, out_ext, out_dir, Rd, yplus)
        end

        if isequal(Qbase[2,2,2,1], NaN) == true
            println("\n NaN !! \n")
            println(" check !!")

            throw(UndefVarError(:x))
            break
        end
        if isequal(yplus[2,2,2], NaN) == true
            println("\n NaN yplus!! \n")
            println(" check yplus!!")

            throw(UndefVarError(:x))
            break
        end

        for i in 2:cellxmax-1
            for j in 2:cellymax-1
                for k in 2:cellzmax-1
                    if Qbase[i,j,k,2] < 1.0e-12
                        Qbase[i,j,k,2] = 0.0
                    end
                end
            end
        end
    end
end

# -- main --
main()

max_stars_repo_name: hide-dog/test

max_issues_repo_name: hide-dog/test

max_forks_repo_name: hide-dog/test


text:
using CUDA
using Enzyme
using Test

function mul_kernel(A)
    i = threadIdx().x
    if i <= length(A)
        A[i] *= A[i]
    end
    return nothing
end

function grad_mul_kernel(A, dA)
    Enzyme.autodiff_deferred(mul_kernel, Const, Duplicated(A, dA))
    return nothing
end

@testset "mul_kernel" begin
    A = CUDA.ones(64,)
    @cuda threads=length(A) mul_kernel(A)
    dA = similar(A)
    dA .= 1
    @cuda threads=length(A) grad_mul_kernel(A, dA)
    @test all(dA .== 2)
end


text:
using Printf

abstract type NonlinearSolver end

config(s::NonlinearSolver) = error("config not implemented for $(typeof(s))")
status(s::NonlinearSolver) = error("status not implemented for $(typeof(s))")
initialize!(s::NonlinearSolver, x₀::AbstractArray) = error("initialize! not implemented for $(typeof(s))")
solver_step!(s::NonlinearSolver) = error("solver_step! not implemented for $(typeof(s))")

function solve!(x, s::NonlinearSolver)
    initialize!(s, x)

    while !meets_stopping_criteria(status(s), config(s))
        next_iteration!(status(s))
        solver_step!(s)
        # residual!(status(s))
    end

    warn_iteration_number(status(s), config(s))

    copyto!(x, solution(status(s)))

    return x
end

struct NonlinearSolverException <: Exception
    msg::String
end

Base.showerror(io::IO, e::NonlinearSolverException) = print(io, "Nonlinear Solver Exception: ", e.msg, "!")

# get_solver_status!(solver::NonlinearSolver{T}, status_dict::Dict) where {T} =
#             get_solver_status!(status(solver), params(solver), status_dict)

# get_solver_status(solver::NonlinearSolver{T}) where {T} = get_solver_status!(solver,
#             Dict(:nls_niter => 0,
#                  :nls_atol => zero(T),
#                  :nls_rtol => zero(T),
#                  :nls_stol => zero(T),
#                  :nls_converged => false)
#             )


text:
module SparseTools

using SparseArrays
using SuiteSparse
using LinearAlgebra

# SPQR begin
include("options.jl")

include("support.jl")
include("wrapper.jl")
include("basic.jl")
include("null.jl")
include("pinv.jl")
include("rank_defaults.jl")
include("rank_deflation.jl")
include("rank_form_basis.jl")
include("ssi.jl")
include("ssp.jl")
# SPQR end

end # module


text:
__precompile__()
module Notes2Graph

import ArgParse; ap = ArgParse
import TextAnalysis: TokenDocument, stem!
import SQLite; sqlite = SQLite
import DataFrames: DataFrame
import Combinatorics: combinations
include("databases.jl")
include("parse.jl")
include("BibTeX.jl")

export initialize_database, find_hashtags, word_stems, add_proposition!, add_t1!, add_t2!, add_t3!, add_t4!, add_t5!, add_t6!, table_length, update!, DataFrame, sqlite, load_database, find_nodeid, related_concepts, find_descrid, descriptions, derivatives, relationid, relation_descr

# main()

end # module


text:
@testset "OpenCL.Device" begin
    @testset "Device Type" begin
        for p in cl.platforms()
            for (t, k) in zip((cl.CL_DEVICE_TYPE_GPU, cl.CL_DEVICE_TYPE_CPU,
                               cl.CL_DEVICE_TYPE_ACCELERATOR, cl.CL_DEVICE_TYPE_ALL),
                              (:gpu, :cpu, :accelerator, :all))

                #for (dk, dt) in zip(cl.devices(p, k), cl.devices(p, t))
                #    @fact dk == dt --> true
                #end
                #devices = cl.devices(p, k)
                #for d in devices
                #    @fact d[:device_type] == t --> true
                #end
            end
        end
    end

    @testset "Device Equality" begin
        for platform in cl.platforms()
            devices = cl.devices(platform)
            if length(devices) > 1
                test_dev = devices[1]
                for dev in devices[2:end]
                   @test pointer(dev) != pointer(test_dev)
                   @test hash(dev) != hash(test_dev)
                   @test isequal(dev, test_dev) == false
               end
           end
       end

    end

   @testset "Device Info" begin
        device_info_keys = Symbol[
                :driver_version,
                :version,
                :extensions,
                :platform,
                :name,
                :device_type,
                :has_image_support,
                :queue_properties,
                :has_queue_out_of_order_exec,
                :has_queue_profiling,
                :has_native_kernel,
                :vendor_id,
                :max_compute_units,
                :max_work_item_size,
                :max_clock_frequency,
                :address_bits,
                :max_read_image_args,
                :max_write_image_args,
                :global_mem_size,
                :max_mem_alloc_size,
                :max_const_buffer_size,
                :local_mem_size,
                :has_local_mem,
                :host_unified_memory,
                :available,
                :compiler_available,
                :max_work_group_size,
                :max_parameter_size,
                :profiling_timer_resolution,
                :max_image2d_shape,
                :max_image3d_shape,
            ]
        for p in cl.platforms()
            if occursin("Portable", p[:name])
                msg = "Skipping Device Info tests for Portable Computing Language Platform "
                @warn(msg)
                continue
            end
            @test isa(p, cl.Platform)
            @test_throws ArgumentError p[:zjdlkf]
            for d in cl.devices(p)
                @test isa(d, cl.Device)
                @test_throws ArgumentError d[:zjdlkf]
                for k in device_info_keys
                    @test d[k] == cl.info(d, k)
                    if k == :extensions
                        @test isa(d[k], Array)
                        if length(d[k]) > 0
                            @test isa(d[k], Array{String, 1})
                        end
                    elseif k == :platform
                        @test d[k] == p
                    elseif k == :max_work_item_sizes
                        @test length(d[k]) == 3
                    elseif k == :max_image2d_shape
                        @test length(d[k]) == 2
                    elseif k == :max_image3d_shape
                        @test length(d[k]) == 3
                    end
                end
            end
        end
    end
end


text:
__precompile__()

module Cbc

if isfile(joinpath(dirname(@__FILE__),"..","deps","deps.jl"))
    include("../deps/deps.jl")
else
    error("Cbc not properly installed. Please run Pkg.build(\"Cbc\")")
end

include("CbcCInterface.jl")
include("CbcSolverInterface.jl")

using Cbc.CbcMathProgSolverInterface
export CbcSolver

end # module


text:
abstract type PMUType end

pmutype(::T) where {T} = error("`pmutype` not defined for arguments of type $T")

# Defaults
_unitstatus(x::PMUType, i...)  = error("`unitstatus` undefined for $(typeof(x))")
_unitcontrol(x::PMUType, i...) = error("`unitcontrol` undefined for $(typeof(x))")
_counter(x::PMUType, i...)     = error("`counter` undefined for $(typeof(x))")
_control(x::PMUType, i...)     = error("`control` undefined for $(typeof(x))")
_extras(x::PMUType, i...)      = error("`extras` undefined for $(typeof(x))")

writetype(x::PMUType) = UInt32

numcounters(x::PMUType)        = error("`numcounters` undefined for $(typeof(x))")
numbytes(x::PMUType) = sizeof(UInt64) * numcounters(x)

unpack(x) = ()
unitstatus(x, i...)     = _unitstatus(pmutype(x), indexzero.((unpack(x)..., i...))...)
unitcontrol(x, i...)    = _unitcontrol(pmutype(x), indexzero.((unpack(x)..., i...))...)
counter(x, i...)        = _counter(pmutype(x), indexzero.((unpack(x)..., i...))...)
control(x, i...)        = _control(pmutype(x), indexzero.((unpack(x)..., i...))...)
extras(x, i...)         = _extras(pmutype(x), indexzero.((unpack(x)..., i...))...)
writetype(x) = writetype(pmutype(x))

numcounters(x) = numcounters(pmutype(x))
numbytes(x) = numbytes(pmutype(x))

### Integrated Memory Controller
struct IMC <: PMUType end
_unitstatus(::IMC)    = IndexZero(0xF8)
_unitcontrol(::IMC)   = IndexZero(0xF4)
_counter(::IMC, i)    = IndexZero(0xA0 + value(i) * 0x8)
_control(::IMC, i)    = IndexZero(0xD8 + value(i) * 0x4)
numcounters(::IMC) = 4

### CHA Counters
struct CHA <: PMUType end
_unitstatus(::CHA, i)   = IndexZero(0xE07 + value(i) * 0x10)
_unitcontrol(::CHA, i)  = IndexZero(0xE00 + value(i) * 0x10)
_counter(::CHA, cha, i) = IndexZero(0xE08 + value(cha) * 0x10 + value(i))
_control(::CHA, cha, i) = IndexZero(0xE01 + value(cha) * 0x10 + value(i))
_extras(x::CHA, cha, i) = IndexZero(0xE05 + value(cha) * 0x10 + value(i))
writetype(::CHA) = UInt64
numcounters(::CHA) = 4

# Customize for various types Specialize
abstract type AbstractUncorePMU end

##### IMC Uncore PMU
# PMU implementation for monitoring the integrated memory controller
struct IMCUncorePMU <: AbstractUncorePMU
    # A handle to the underlying
    handle::Handle

    # Pre-allocated buffer for reading new counter values.
    # We return counter values as a tuple for even better performance.
    buffer::Vector{UInt8}
end
IMCUncorePMU(handle::Handle) = IMCUncorePMU(handle, zeros(UInt8, numbytes(IMC())))
pmutype(::IMCUncorePMU) = IMC()
Base.close(x::IMCUncorePMU) = close(x.handle)

##### CHA Uncore PMU
# PMU implementation for monitoring the CHA
struct CHAUncorePMU <: AbstractUncorePMU
    # We hold on to a single handle for the MSR path, shared by all PMUs
    handle::Handle
    # The number of this CHA
    cha::IndexZero{Int}
    buffer::Vector{UInt8}

    # Allow passing a buffer, or manually create one
    function CHAUncorePMU(
            handle::Handle,
            cha,
            buffer = zeros(UInt8, numbytes(CHA()))
        )
        resize!(buffer, numbytes(CHA()))

        return new(handle, indexzero(cha), buffer)
    end
end

pmutype(::CHAUncorePMU) = CHA()
unpack(x::CHAUncorePMU) = (x.cha,)
Base.close(x::CHAUncorePMU) = close(x.handle)

#####
##### Low level accessing functions
#####

function setunitstatus!(U::AbstractUncorePMU, v)
    seek(U.handle, unitstatus(U))
    write(U.handle, convert(writetype(U), v))
end

function getunitstatus(U::AbstractUncorePMU)
    seek(U.handle, unitstatus(U))
    return read(U.handle, UInt32)
end

function setunitcontrol!(U::AbstractUncorePMU, v)
    seek(U.handle, unitcontrol(U))
    write(U.handle, convert(writetype(U), v))
end

function getunitcontrol(U::AbstractUncorePMU)
    seek(U.handle, unitcontrol(U))
    return read(U.handle, UInt32)
end

function setcontrol!(U::AbstractUncorePMU, counter, v)
    seek(U.handle, control(U, counter))
    write(U.handle, convert(writetype(U), v))
end

function getcontrol(U::AbstractUncorePMU, i)
    seek(U.handle, control(U, i))
    return read(U.handle, UInt32)
end

function getcounter(U::AbstractUncorePMU, i)
    seek(U.handle, counter(U, i))
    return CounterValue(read(U.handle, UInt64))
end

function setextra!(U::AbstractUncorePMU, i, v)
    seek(U.handle, extras(U, i))
    write(U.handle, convert(writetype(U), v))
end

function getextra(U::AbstractUncorePMU, i)
    seek(U.handle, extras(U, i))
    return read(U.handle, UInt32)
end

#####
##### Some higher level functions
#####

function getallcounters(U::AbstractUncorePMU)
    # Need to seek and read since MSR based registers don't automatically progress
    # the position in the system file.
    a = unsafe_read(U.handle, UInt64, counter(U, 1); buffer = U.buffer)
    b = unsafe_read(U.handle, UInt64, counter(U, 2); buffer = U.buffer)
    c = unsafe_read(U.handle, UInt64, counter(U, 3); buffer = U.buffer)
    d = unsafe_read(U.handle, UInt64, counter(U, 4); buffer = U.buffer)
    return CounterValue.((a, b, c, d))
end

function reset!(U::AbstractUncorePMU)
    # Write to the unit control to clear all counters and control registers
    val = setbits(zero(writetype(U)), (0, 1, 8, 16, 17))
    setunitcontrol!(U, val)
end

function enable!(U::AbstractUncorePMU)
    val = setbits(zero(writetype(U)), (16, 17))
    setunitcontrol!(U, val)
end


text:
function ContinuedFraction(α,β,γ,n::Int,N::Int; tol = 1e-10, Nmin = 20, rN = zero(α(1)))
    num = β(0)
    for i = 1:n
        num = β(i) - α(i-1)*γ(i)/num
    end
    lhs = num

    num = β(N)+α(N)*rN
    for i = (N):-1:(n+1)
       num = β(i-1) - α(i-1)*γ(i)/num
    end
    rhs = num - β(n)

    lhs + rhs
end

function rNCoeffs(D0,D1,D2,D3,D4)
    fp = sqrt(-D0 - D1 - D2)
    u1 = -fp
    u2 = (1/2)*(-4 - 2*D0 - D1)
    u3 = (fp)*(8 + 16*D0 + 8*(D0^2) + 12*D1 + 8*D0*D1 + (D1^2) + 8*D2 +
    4*D0*D2 - 4*D3 - 4*D4)/(8*(D0 + D1 + D2))
    u4 = (1/2)*(4 + 4*D0 + 2*D0^2 + D1 + D0*D1 - D3)
    u1,u2,u3,u4
end

function CF(P)
    _,_,(D₀,D₁,D₂,D₃,D₄) = ParameterTransformations(P)
    n = P.n; N = P.Nmax;
    αₙ(n) = n^2 + (D₀+1)*n + D₀
    βₙ(n) = -2*n^2 + (D₁+2)*n +D₃
    γₙ(n) = n^2 + (D₂-3)*n  + D₄ - D₂ + 2
    u1,u2,u3,u4 = rNCoeffs(D₀,D₁,D₂,D₃,D₄)
    rN = 1 + u1*N^(-0.5) + u2*N^(-1) + u3*N^(1.5) + u4*N^(-2)
    ContinuedFraction(αₙ,βₙ,γₙ,n,N; rN = rN)
end


text:
################################################################################
#                                BASIC METHODS
################################################################################
# Macro error function
# inspired from https://github.com/jump-dev/JuMP.jl/blob/709d41b78e56efb4f2c54414266b30932010bd5a/src/macros.jl#L923-L928
function _macro_error(macroname, args, source, str...)
    error("At $(source.file):$(source.line): `@$macroname($(join(args, ", ")))`: ",
          str...)
end

# Escape when needed
# taken from https://github.com/jump-dev/JuMP.jl/blob/709d41b78e56efb4f2c54414266b30932010bd5a/src/macros.jl#L895-L897
_esc_non_constant(x::Number) = x
_esc_non_constant(x::Expr) = isexpr(x,:quote) ? x : esc(x)
_esc_non_constant(x) = esc(x)

# Given a base_name and idxvars, returns an expression that constructs the name
# of the object.
# Taken from https://github.com/jump-dev/JuMP.jl/blob/709d41b78e56efb4f2c54414266b30932010bd5a/src/macros.jl#L930-L946
function _name_call(base_name, idxvars)
    if isempty(idxvars) || base_name == ""
        return base_name
    end
    ex = Expr(:call, :string, base_name, "[")
    for i in 1:length(idxvars)
        # Converting the arguments to strings before concatenating is faster:
        # https://github.com/JuliaLang/julia/issues/29550.
        esc_idxvar = esc(idxvars[i])
        push!(ex.args, :(string($esc_idxvar)))
        i < length(idxvars) && push!(ex.args, ",")
    end
    push!(ex.args, "]")
    return ex
end

# Process macro arugments
# Taken from https://github.com/jump-dev/JuMP.jl/blob/709d41b78e56efb4f2c54414266b30932010bd5a/src/Containers/macro.jl#L26-L36
function _extract_kw_args(args)
    kw_args = filter(x -> isexpr(x, :(=)) && x.args[1] != :container , collect(args))
    flat_args = filter(x->!isexpr(x, :(=)), collect(args))
    requested_container = :Auto
    for kw in args
        if isexpr(kw, :(=)) && kw.args[1] == :container
            requested_container = kw.args[2]
        end
    end
    return flat_args, kw_args, requested_container
end


text:
function FESCnmSnmReader(arg0::JString, arg1::jdouble)
    return FESCnmSnmReader((JString, jdouble), arg0, arg1)
end

function load_data(obj::FESCnmSnmReader, arg0::InputStream, arg1::JString)
    return jcall(obj, "loadData", void, (InputStream, JString), arg0, arg1)
end


text:
#### write to files

open("lumped_"*filename, "w") do f
  write(f, "{Lumping $(filename) using lumper.jl tool. D.Ellis 2017\n\n")

  write(f,"Add defns to constants file \n")

nlumped = length(collect(Base.flatten(lumplist)))

  write(f,"\n\nSpecies - Lumping combinations are : \n\n")

  for i in enumerate(lumplist)
    write(f, " LMP$i \n")
  end

  write(f,"}\n\n")

  m = open(filename);
  for ln in eachline(m)
    #println("$(length(ln)), $ln")
    write(f, ln*"")
    if contains(ln,"DEFVAR")
      p("stop")
      break
    end
  end
  close(m)

  writespecs = []
  for x in specs
    try
      x = "LMP$(lumpdict[x])"
    end
    push!(writespecs,x)
  end
  writespecs = (unique(writespecs,1))
  sort!(writespecs)

  for s in writespecs
    if s != "DUMMY"
      write(f, "$s = IGNORE;\n")
    end
  end

#greater than 80 characters

  write(f, """#INLINE F90_RCONST
  USE constants
  USE lumpedrates
  ! end of USE statements
  ! start of executable statements

  """)

write(f,"Real(dp)::PA($(nlumped))\n")
write(f,"Real(dp)::LC($(nlumped))\n")
write(f,"Integer::index\n")

  write(f,join(ppair,"\n    "))

    write(f, """\n\n
    do index = 1, $(nlumped)
      if (isnan(LC(index))) LC(index) = 0
    end do
    """)

  write(f, "\n\nRO2 = 0 ")

  ro2 = sort([replace(x,"ind_","") for x in  matchall(r"ind_(\w+)\b", join(lines) ) ])
  #groupspecs = collect(Base.flatten(lumplist))
  for i in 1:length(ro2)

    try
      ro2[i] = "C(IND_LMP$(lumpspec[ro2[i]])) * LC($(lumpspec[ro2[i]]))"
    catch
      ro2[i] = "C(IND_$(ro2[i]))"
    end
    write(f, "&\n+ $(ro2[i]) ")

  end

  write(f, """\n\n    CALL mcm_constants(time, temp, M, N2, O2, RO2, H2O)
  #ENDINLINE """)

  write(f,"\n\n#EQUATIONS\n")

  for i in enumerate(eqns)

    rt = i[2][4]
    #rt = tostr(sympy.N(sympy.expand(replace(i[2][4],r"(\d)D([\+\-\d])",s"\1E\2")),3))

    rt = replace(rt,r"EXP\(([-\d\.]+)/TEMP\)",s"ET(\1)")
    rt = replace(rt,r"LC_(\d)_",s"LC(\1)")
    write(f,"{$(i[1])} $(i[2][2]) = $(i[2][3]) : $(rt);\n")
  end

end

run(`scp ./lumped_methanecomplete.kpp dp626@earth0:~/DSMACC-testing/mechanisms/autolump.kpp`)

"""
module lumpedrates
USE model_Precision, ONLY: dp
use constants
IMPLICIT NONE

contains
 FUNCTION ET(a)
    USE model_global, only: TEMP
    real(dp) :: ET
    INTEGER    :: a
    ET = EXP(a/TEMP)
 END FUNCTION

end module lumpedrates
"""

max_stars_repo_name: wolfiex/MechLump

max_issues_repo_name: wolfiex/MechLump

max_forks_repo_name: wolfiex/MechLump


text:
function get_a1(obj::RotationOrder)
    return jcall(obj, "getA1", Vector3D, ())
end

function get_a2(obj::RotationOrder)
    return jcall(obj, "getA2", Vector3D, ())
end

function get_a3(obj::RotationOrder)
    return jcall(obj, "getA3", Vector3D, ())
end

function get_rotation_order(::Type{RotationOrder}, arg0::JString)
    return jcall(RotationOrder, "getRotationOrder", RotationOrder, (JString,), arg0)
end

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

function value_of(::Type{RotationOrder}, arg0::JString)
    return jcall(RotationOrder, "valueOf", RotationOrder, (JString,), arg0)
end

function values(::Type{RotationOrder})
    return jcall(RotationOrder, "values", Vector{RotationOrder}, ())
end


text:
# TODO: Remove some T's from output type signatures

"""
    DataArray{T,N}(d::Array{T,N}, m::AbstractArray{Bool} = falses(size(d)))

Construct a `DataArray`, an `N`-dimensional array with element type `T` that allows missing
values. The resulting array uses the data in `d` with `m` as a bitmask to signify missingness.
That is, for each index `i` in `d`, if `m[i]` is `true`, the array contains `missing` at index `i`,
otherwise it contains `d[i]`.

    DataArray(T::Type, dims...)

Construct a `DataArray` with element type `T` and dimensions specified by `dims`. All elements
default to `missing`.

# Examples

```jldoctest
julia> DataArray([1, 2, 3], [true, false, true])
3-element DataArrays.DataArray{Int64,1}:
  missing
 2
  missing

julia> DataArray(Float64, 3, 3)
3×3 DataArrays.DataArray{Float64,2}:
 missing  missing  missing
 missing  missing  missing
 missing  missing  missing
```
"""
mutable struct DataArray{T, N} <: AbstractDataArray{T, N}
    data::Array{T, N}
    na::BitArray{N}

    function DataArray{T,N}(d::Array{<:Union{T, Missing}, N}, m::BitArray{N}) where {T, N}
        # Ensure data values and missingness metadata match
        if size(d) != size(m)
            msg = "Data and missingness arrays must be the same size"
            throw(ArgumentError(msg))
        end
        # if input array can contain missing values, we need to mark corresponding entries as missing
        if eltype(d) >: Missing
            # If the original eltype is wider than the target eltype T, conversion may fail
            # in the presence of missings: we need to allocate a copy, leaving entries
            # corresponding to missings uninitialized
            if eltype(d) <: T
                @inbounds for i in eachindex(d)
                    if isassigned(d, i) && ismissing(d, i)
                        m[i] = true
                    end
                end
            else
                d2 = similar(d, T)
                @inbounds for i in eachindex(d)
                    isassigned(d, i) || continue
                    if ismissing(d, i)
                        m[i] = true
                    else
                        d2[i] = d[i]
                    end
                end
                return new(d2, m)
            end
        elseif eltype(d) <: Missing
            m = trues(m)
        end
        new(d, m)
    end
end

DataArray{T}(d::Array{S, N}) where {T, S, N} = DataArray{T, N}(d) # -> DataArray{T}

function DataArray{T, N}(d::Array,
                         m::BitArray{N} = falses(size(d))) where {T, N}  # -> DataArray{T}
    DataArray(convert(Array{Missings.T(T), N}, d), m)
end

function DataArray(d::Array{T, N},
                   m::BitArray{N} = falses(size(d))) where {T, N} # -> DataArray{T}
    return DataArray{Missings.T(T), N}(d, m)
end

function DataArray(d::Array, m::AbstractArray{Bool}) # -> DataArray{T}
    return DataArray(d, BitArray(m))
end

function DataArray(T::Type, dims::Integer...) # -> DataArray{T}
    return DataArray(Array{Missings.T(T)}(uninitialized, dims...), trues(dims...))
end

function DataArray(T::Type, dims::NTuple{N, Int}) where N # -> DataArray{T}
    return DataArray(Array{Missings.T(T)}(dims...), trues(dims...))
end

"""
    DataVector{T}

A 1-dimensional `DataArray` with element type `T`.
"""
const DataVector{T} = DataArray{T, 1}

"""
    DataMatrix{T}

A 2-dimensional `DataArray` with element type `T`.
"""
const DataMatrix{T} = DataArray{T, 2}

Base.copy(d::DataArray) = Base.copy!(similar(d), d) # -> DataArray{T}

function Base.copy!(dest::DataArray, src::DataArray) # -> DataArray{T}
    if isbits(eltype(src)) && isbits(eltype(dest))
        copy!(dest.data, src.data)
    else
        # Elements of src_data are not necessarily initialized, so
        # only copy initialized elements
        dest_data = dest.data
        src_data = src.data
        length(dest_data) >= length(src_data) || throw(BoundsError())
        src_chunks = src.na.chunks
        for i = 1:length(src_data)
            @inbounds if !Base.unsafe_bitgetindex(src_chunks, i)
                dest_data[i] = src_data[i]
            end
        end
    end
    copy!(dest.na, src.na)
    dest
end

function Base.copy!(dest::DataArray, doffs::Integer, src::DataArray) # -> DataArray{T}
    copy!(dest, doffs, src, 1, length(src))
end

# redundant on Julia 0.4
function Base.copy!(dest::DataArray, doffs::Integer, src::DataArray, soffs::Integer) # -> DataArray{T}
    soffs <= length(src) || throw(BoundsError())
    copy!(dest, doffs, src, soffs, length(src)-soffs+1)
end

function Base.copy!(dest::DataArray, doffs::Integer, src::DataArray, soffs::Integer, n::Integer) # -> DataArray{T}
    if n == 0
        return dest
    elseif n < 0
        throw(ArgumentError("tried to copy n=$n elements, but n should be nonnegative"))
    end
    if isbits(eltype(src))
        copy!(dest.data, doffs, src.data, soffs, n)
    else
        # Elements of src_data are not necessarily initialized, so
        # only copy initialized elements
        dest_data = dest.data
        src_data = src.data
        if doffs < 1 || length(dest_data) - doffs + 1 < n ||
           soffs < 1 || length(src_data) - soffs + 1 < n
            throw(BoundsError())
        end
        src_chunks = src.na.chunks
        for i = 0:(n-1)
            di, si = doffs + i, soffs + i

            @inbounds if !Base.unsafe_bitgetindex(src_chunks, si)
                dest_data[di] = src_data[si]
            end
        end
    end
    copy!(dest.na, doffs, src.na, soffs, n)
    dest
end

Base.fill!(A::DataArray, ::Missing) = (fill!(A.na, true); A)
Base.fill!(A::DataArray, v) = (fill!(A.data, v); fill!(A.na, false); A)

function Base.deepcopy(d::DataArray) # -> DataArray{T}
    return DataArray(deepcopy(d.data), deepcopy(d.na))
end

function Base.resize!(da::DataArray{T,1}, n::Int) where T
    resize!(da.data, n)
    oldn = length(da.na)
    resize!(da.na, n)
    da.na[oldn+1:n] = true
    da
end

function Base.similar(da::DataArray, T::Type, dims::Dims) #-> DataArray{T}
    return DataArray(Array{Missings.T(T)}(uninitialized, dims), trues(dims))
end

Base.size(d::DataArray) = size(d.data) # -> (Int...)
Base.ndims(da::DataArray) = ndims(da.data) # -> Int
Base.length(d::DataArray) = length(d.data) # -> Int
Base.endof(da::DataArray) = endof(da.data) # -> Int

function Base.find(da::DataArray{Bool}) # -> Array{Int}
    data = da.data
    ntrue = 0
    @inbounds @bitenumerate da.na i na begin
        ntrue += !na && data[i]
    end
    res = Vector{Int}(ntrue)
    count = 1
    @inbounds @bitenumerate da.na i na begin
        if !na && data[i]
            res[count] = i
            count += 1
        end
    end
    return res
end

function Base.convert(::Type{Array{S, N}},
                      x::DataArray{T, N}) where {S, T, N} # -> Array{S, N}
    return S[v for v in x]
end

function Base.convert(::Type{Array{S}}, da::DataArray{T, N}) where {S, T, N}
    return convert(Array{S, N}, da)
end

function Base.convert(::Type{Vector}, dv::DataVector{T}) where T
    return convert(Array{Union{T, Missing}, 1}, dv)
end

function Base.convert(::Type{Matrix}, dm::DataMatrix{T}) where T
    return convert(Array{Union{T, Missing}, 2}, dm)
end

function Base.convert(::Type{Array}, da::DataArray{T, N}) where {T, N}
    return convert(Array{Union{T, Missing}, N}, da)
end

function Base.convert(
    ::Type{Array{S, N}},
    da::DataArray{T, N},
    replacement::Any
) where {S, T, N} # -> Array{S, N}

    replacementS = convert(S, replacement)
    res = Array{S}(uninitialized, size(da))
    for i in 1:length(da)
        if da.na[i]
            res[i] = replacementS
        else
            res[i] = da.data[i]
        end
    end
    return res
end

function Base.convert(::Type{Vector}, dv::DataVector{T}, replacement::Any) where T
    return convert(Array{T, 1}, dv, replacement)
end

function Base.convert(::Type{Matrix}, dm::DataMatrix{T}, replacement::Any) where T
    return convert(Array{T, 2}, dm, replacement)
end

function Base.convert(::Type{Array}, da::DataArray{T, N}, replacement::Any) where {T, N}
    return convert(Array{T, N}, da, replacement)
end

struct EachFailMissing{T<:DataArray}
    da::T
end
Missings.fail(da::DataArray) = EachFailMissing(da)
Base.length(itr::EachFailMissing) = length(itr.da)
Base.start(itr::EachFailMissing) = 1
Base.done(itr::EachFailMissing, ind::Integer) = ind > length(itr)
Base.eltype(itr::EachFailMissing) = Missings.T(eltype(itr.da))
function Base.next(itr::EachFailMissing, ind::Integer)
    if itr.da.na[ind]
        throw(MissingException("missing value encountered in Missings.fail"))
    else
        (itr.da.data[ind], ind + 1)
    end
end

struct EachDropMissing{T<:DataArray}
    da::T
end
Missings.skipmissing(da::DataArray) = EachDropMissing(da)
function _next_nonna_ind(da::DataArray, ind::Int)
    ind += 1
    @inbounds while ind <= length(da) && da.na[ind]
        ind += 1
    end
    ind
end
Base.length(itr::EachDropMissing) = length(itr.da) - sum(itr.da.na)
Base.start(itr::EachDropMissing) = _next_nonna_ind(itr.da, 0)
Base.done(itr::EachDropMissing, ind::Int) = ind > length(itr.da)
Base.eltype(itr::EachDropMissing) = Missings.T(eltype(itr.da))
function Base.next(itr::EachDropMissing, ind::Int)
    (itr.da.data[ind], _next_nonna_ind(itr.da, ind))
end

struct EachReplaceMissing{S<:DataArray, T}
    da::S
    replacement::T
end
Missings.replace(da::DataArray, replacement::Any) =
    EachReplaceMissing(da, replacement)
Base.length(itr::EachReplaceMissing) = length(itr.da)
Base.start(itr::EachReplaceMissing) = 1
Base.done(itr::EachReplaceMissing, ind::Integer) = ind > length(itr)
Base.eltype(itr::EachReplaceMissing) = Missings.T(eltype(itr.da))
function Base.next(itr::EachReplaceMissing, ind::Integer)
    item = itr.da.na[ind] ? itr.replacement : itr.da.data[ind]
    (item, ind + 1)
end

Base.collect(itr::EachDropMissing{<:DataVector}) = itr.da.data[.!itr.da.na] # -> Vector
Base.collect(itr::EachFailMissing{<:DataVector}) = copy(itr.da.data) # -> Vector

Base.any(::typeof(ismissing), da::DataArray) = any(da.na) # -> Bool
Base.all(::typeof(ismissing), da::DataArray) = all(da.na) # -> Bool

Missings.ismissing(da::DataArray, I::Real, Is::Real...) = getindex(da.na, I, Is...)

function Base.isfinite(da::DataArray) # -> DataArray{Bool}
    n = length(da)
    res = Array{Bool}(size(da))
    for i in 1:n
        if !da.na[i]
            res[i] = isfinite(da.data[i])
        end
    end
    return DataArray(res, copy(da.na))
end

# Promotion rules

# promote_rule{T, T}(::Type{AbstractDataArray{T}},
#                    ::Type{T}) = promote_rule(T, T)
# promote_rule{S, T}(::Type{AbstractDataArray{S}},
#                    ::Type{T}) = promote_rule(S, T)
# promote_rule{T}(::Type{AbstractDataArray{T}}, ::Type{T}) = T

function Base.convert(::Type{DataArray{S, N}},
                      a::AbstractArray{T, N}) where {S, T, N} # -> DataArray{S, N}
    return DataArray(convert(Array{S, N}, a), falses(size(a)))
end

function Base.convert(::Type{DataArray{S, N}},
                      a::AbstractArray{T, N}) where {S, T>:Missing, N} # -> DataArray{S, N}
    return DataArray(convert(Array{Union{S, Missing}, N}, a), falses(size(a)))
end

function Base.convert(::Type{DataArray{S}}, x::AbstractArray{T, N}) where {S, T, N}
    convert(DataArray{Missings.T(S), N}, x)
end

function Base.convert(::Type{DataArray}, x::AbstractArray{T, N}) where {T, N}
    convert(DataArray{Missings.T(T), N}, x)
end

Base.convert(::Type{DataArray{T, N}}, x::DataArray{T, N}) where {T, N} = x

function Base.convert(::Type{DataArray{S, N}},
                      x::DataArray{T, N}) where {S, T, N} # -> DataArray{S, N}
    v = similar(x.data, S)
    @inbounds for i = 1:length(x)
        if !x.na[i]
            v[i] = convert(S, x.data[i])
        end
    end
    return DataArray(v, x.na)
end

"""
    data(a::AbstractArray) -> DataArray

Convert `a` to a `DataArray`.

# Examples

```jldoctest
julia> data([1, 2, 3])
3-element DataArrays.DataArray{Int64,1}:
 1
 2
 3

julia> data(@data [1, 2, missing])
3-element DataArrays.DataArray{Int64,1}:
 1
 2
  missing
```
"""
data(a::AbstractArray) = convert(DataArray, a)

# TODO: Make sure these handle copying correctly
# TODO: Remove these? They have odd behavior, because they convert to Array's.
# TODO: Rethink multi-item documentation approach
for f in (:(Base.float),)
    @eval begin
        function ($f)(da::DataArray) # -> DataArray
            if any(ismissing, da)
                err = "Cannot convert DataArray with missings to desired type"
                throw(MissingException(err))
            else
                ($f)(da.data)
            end
        end
    end
end

"""
    finduniques(da::DataArray) -> (Vector, Int)

Get the unique values in `da` as well as the index of the first `missing` value
in `da` if present, or 0 otherwise.
"""
function finduniques(da::DataArray{T}) where T # -> Vector{T}, Int
    out = Vector{T}(uninitialized, 0)
    seen = Set{T}()
    n = length(da)
    firstmissing = 0
    for i in 1:n
        if ismissing(da, i)
            if firstmissing == 0
                firstmissing = length(out) + 1
            else
                continue
            end
        elseif !in(da.data[i], seen)
            push!(seen, da.data[i])
            push!(out, da.data[i])
        end
    end
    return out, firstmissing
end

function Base.unique(da::DataArray{T}) where T # -> DataVector{T}
    unique_values, firstmissing = finduniques(da)
    n = length(unique_values)
    if firstmissing > 0
        res = DataArray(Vector{T}(uninitialized, n + 1))
        i = 1
        for val in unique_values
            if i == firstmissing
                res.na[i] = true
                i += 1
            end
            res.data[i] = val
            i += 1
        end

        if firstmissing == n + 1
            res.na[n + 1] = true
        end

        return res
    else
        return DataArray(unique_values)
    end
end

function Missings.levels(da::DataArray) # -> DataVector{T}
    unique_values, firstmissing = finduniques(da)
    return DataArray(unique_values)
end


text:
@testset "SimpleCovariance" begin
    v = simple()
    @test sprint(show, v) == "Simple covariance estimator"
end

@testset "RobustCovariance" begin
    v = robust()
    @test sprint(show, v) == "Heteroskedasticity-robust covariance estimator"
end

@testset "ClusterCovariance" begin
    @test_throws MethodError cluster()
    c1 = cluster(:a)
    @test names(c1) == (:a,)
    @test length(c1) == 1
    c2 = cluster(:a, :b)
    @test names(c2) == (:a, :b)
    @test length(c2) == 2

    @test sprint(show, c1) == "Cluster-robust covariance estimator"
    @test sprint(show, MIME("text/plain"), c1) == """
        1-way cluster-robust covariance estimator:
          a"""
    @test sprint(show, MIME("text/plain"), c2) == """
        2-way cluster-robust covariance estimator:
          a
          b"""

    N = 10
    a1 = collect(1:N)
    g1 = group(a1)
    c1 = cluster((:a,), (g1,))
    @test nclusters(c1) == (a=N,)
end