text:
#TODO: Move stdio detection from C to Julia (might require some Clang magic)
include("uv_constants.jl")

## types ##
typealias Callback Union(Function,Bool)

abstract AsyncStream <: IO
abstract UVServer

typealias UVHandle Ptr{Void}
typealias UVStream AsyncStream

# A dict of all libuv handles that are being waited on somewhere in the system
# and should thus not be garbage collected
const uvhandles = ObjectIdDict()

preserve_handle(x) = uvhandles[x] = get(uvhandles,x,0)+1
unpreserve_handle(x) = (v = uvhandles[x]; v == 1 ? pop!(uvhandles,x) : (uvhandles[x] = v-1); nothing)

#Wrapper for an OS file descriptor (on both Unix and Windows)
immutable RawFD
    fd::Int32
    RawFD(fd::Integer) = new(int32(fd))
    RawFD(fd::RawFD) = fd
end

convert(::Type{Int32}, fd::RawFD) = fd.fd

function uv_sizeof_handle(handle)
    if !(UV_UNKNOWN_HANDLE < handle < UV_HANDLE_TYPE_MAX)
        throw(DomainError())
    end
    ccall(:uv_handle_size,Csize_t,(Int32,),handle)
end

function uv_sizeof_req(req)
    if !(UV_UNKNOWN_REQ < req < UV_REQ_TYPE_MAX)
        throw(DomainError())
    end
    ccall(:uv_req_size,Csize_t,(Int32,),req)
end

for h in uv_handle_types
@eval const $(symbol("_sizeof_"*lowercase(string(h)))) = uv_sizeof_handle($h)
end
for r in uv_req_types
@eval const $(symbol("_sizeof_"*lowercase(string(r)))) = uv_sizeof_req($r)
end

function eof(s::AsyncStream)
    wait_readnb(s,1)
    !isopen(s) && nb_available(s.buffer)<=0
end

const StatusUninit      = 0 # handle is allocated, but not initialized
const StatusInit        = 1 # handle is valid, but not connected/active
const StatusConnecting  = 2 # handle is in process of connecting
const StatusOpen        = 3 # handle is usable
const StatusActive      = 4 # handle is listening for read/write/connect events
const StatusClosing     = 5 # handle is closing / being closed
const StatusClosed      = 6 # handle is closed
function uv_status_string(x)
    s = x.status
    if x.handle == C_NULL
        if s == StatusClosed
            return "closed"
        elseif s == StatusUninit
            return "null"
        end
        return "invalid status"
    elseif s == StatusUninit
        return "uninit"
    elseif s == StatusInit
        return "init"
    elseif s == StatusConnecting
        return "connecting"
    elseif s == StatusOpen
        return "open"
    elseif s == StatusActive
        return "active"
    elseif s == StatusClosing
        return "closing"
    elseif s == StatusClosed
        return "closed"
    end
    return "invalid status"
end

uv_handle_data(handle) = ccall(:jl_uv_handle_data,Ptr{Void},(Ptr{Void},),handle)
uv_req_data(handle) = ccall(:jl_uv_req_data,Ptr{Void},(Ptr{Void},),handle)
uv_req_set_data(req,data) = ccall(:jl_uv_req_set_data,Void,(Ptr{Void},Any),req,data)
uv_req_set_data(req,data::Ptr{Void}) = ccall(:jl_uv_req_set_data,Void,(Ptr{Void},Ptr{Void}),req,data)

type Pipe <: AsyncStream
    handle::Ptr{Void}
    status::Int
    buffer::IOBuffer
    line_buffered::Bool
    readcb::Callback
    readnotify::Condition
    ccb::Callback
    connectnotify::Condition
    closecb::Callback
    closenotify::Condition
    Pipe(handle) = new(
        handle,
        StatusUninit,
        PipeBuffer(),
        true,
        false,Condition(),
        false,Condition(),
        false,Condition())
end
function Pipe()
    handle = c_malloc(_sizeof_uv_named_pipe)
    try
        ret = Pipe(handle)
        associate_julia_struct(ret.handle,ret)
        finalizer(ret,uvfinalize)
        return init_pipe!(ret;readable=true)
    catch
        c_free(handle)
        rethrow()
    end
end

type PipeServer <: UVServer
    handle::Ptr{Void}
    status::Int
    ccb::Callback
    connectnotify::Condition
    closecb::Callback
    closenotify::Condition
    PipeServer(handle) = new(
        handle,
        StatusUninit,
        false,Condition(),
        false,Condition())
end

function init_pipe!(pipe::Union(Pipe,PipeServer);readable::Bool=false,writable=false,julia_only=true)
    if pipe.handle == C_NULL
        error("failed to initialize pipe")
    elseif pipe.status != StatusUninit
        error("pipe is already initialized")
    end
    uv_error("init_pipe",ccall(:jl_init_pipe, Cint, (Ptr{Void},Int32,Int32,Int32), pipe.handle, writable,readable,julia_only))
    pipe.status = StatusInit
    pipe
end

function PipeServer()
    handle = c_malloc(_sizeof_uv_named_pipe)
    try
        ret = PipeServer(handle)
        associate_julia_struct(ret.handle,ret)
        finalizer(ret,uvfinalize)
        return init_pipe!(ret;readable=true)
    catch
        c_free(handle)
        c_free(handle-1)
        rethrow()
    end
end

show(io::IO,stream::Pipe) = print(io,"Pipe(",uv_status_string(stream),", ",
    nb_available(stream.buffer)," bytes waiting)")
show(io::IO,stream::PipeServer) = print(io,"PipeServer(",uv_status_string(stream),")")

type TTY <: AsyncStream
    handle::Ptr{Void}
    status::Int
    line_buffered::Bool
    buffer::IOBuffer
    readcb::Callback
    readnotify::Condition
    closecb::Callback
    closenotify::Condition
    TTY(handle) = new(
        handle,
        StatusUninit,
        true,
        PipeBuffer(),
        false,Condition(),
        false,Condition())
end

function TTY(fd::RawFD; readable::Bool = false)
    handle = c_malloc(_sizeof_uv_tty)
    ret = TTY(handle)
    associate_julia_struct(handle,ret)
    finalizer(ret,uvfinalize)
    # This needs to go after associate_julia_struct so that there
    # is no garbage in the ->data field
    uv_error("TTY",ccall(:uv_tty_init,Int32,(Ptr{Void},Ptr{Void},Int32,Int32),eventloop(),handle,fd.fd,readable))
    ret.status = StatusOpen
    ret.line_buffered = false
    ret
end

# note that uv_is_readable/writable work for any subtype of
# uv_stream_t, including uv_tty_t and uv_pipe_t
isreadable(io::Union(Pipe,TTY)) =
    bool(ccall(:uv_is_readable, Cint, (Ptr{Void},), io.handle))
iswritable(io::Union(Pipe,TTY)) =
    bool(ccall(:uv_is_writable, Cint, (Ptr{Void},), io.handle))

nb_available(stream::UVStream) = nb_available(stream.buffer)

show(io::IO,stream::TTY) = print(io,"TTY(",uv_status_string(stream),", ",
    nb_available(stream.buffer)," bytes waiting)")

uvtype(::AsyncStream) = UV_STREAM
uvhandle(stream::AsyncStream) = stream.handle

convert(T::Type{Ptr{Void}}, s::AsyncStream) = convert(T, s.handle)
handle(s::AsyncStream) = s.handle
handle(s::Ptr{Void}) = s

make_stdout_stream() = _uv_tty2tty(ccall(:jl_stdout_stream, Ptr{Void}, ()))

associate_julia_struct(handle::Ptr{Void},jlobj::ANY) =
    ccall(:jl_uv_associate_julia_struct,Void,(Ptr{Void},Any),handle,jlobj)
disassociate_julia_struct(uv) = disassociate_julia_struct(uv.handle)
disassociate_julia_struct(handle::Ptr{Void}) =
    handle != C_NULL && ccall(:jl_uv_disassociate_julia_struct,Void,(Ptr{Void},),handle)

function init_stdio(handle)
    t = ccall(:jl_uv_handle_type,Int32,(Ptr{Void},),handle)
    if t == UV_FILE
        return fdio(ccall(:jl_uv_file_handle,Int32,(Ptr{Void},),handle))
    else
        if t == UV_TTY
            ret = TTY(handle)
        elseif t == UV_TCP
            ret = TcpSocket(handle)
        elseif t == UV_NAMED_PIPE
            ret = Pipe(handle)
        else
            error("FATAL: stdio type ($t) invalid")
        end
        ret.status = StatusOpen
        ret.line_buffered = false
        associate_julia_struct(ret.handle,ret)
        finalizer(ret,uvfinalize)
        return ret
    end
end

function stream_wait(x,c...)
    preserve_handle(x)
    try
        return wait(c...)
    finally
        unpreserve_handle(x)
    end
end

function reinit_stdio()
    global uv_jl_asynccb = cglobal(:jl_uv_asynccb)
    global uv_jl_alloc_buf = cglobal(:jl_uv_alloc_buf)
    global uv_jl_readcb = cglobal(:jl_uv_readcb)
    global uv_jl_connectioncb = cglobal(:jl_uv_connectioncb)
    global uv_jl_connectcb = cglobal(:jl_uv_connectcb)
    global uv_jl_writecb = cglobal(:jl_uv_writecb)
    global uv_jl_writecb_task = cglobal(:jl_uv_writecb_task)
    global uv_eventloop = ccall(:jl_global_event_loop, Ptr{Void}, ())
    global STDIN = init_stdio(ccall(:jl_stdin_stream ,Ptr{Void},()))
    global STDOUT = init_stdio(ccall(:jl_stdout_stream,Ptr{Void},()))
    global STDERR = init_stdio(ccall(:jl_stderr_stream,Ptr{Void},()))
end

flush(::AsyncStream) = nothing

function isopen(x)
    if !(x.status != StatusUninit && x.status != StatusInit)
        error("I/O object not initialized")
    end
    x.status != StatusClosed
end

function check_open(x)
    if !isopen(x)
        error("stream is closed or unusable")
    end
end

function wait_connected(x)
    check_open(x)
    while x.status == StatusConnecting
        stream_wait(x,x.connectnotify)
        check_open(x)
    end
end

function wait_readbyte(x::AsyncStream, c::Uint8)
    while isopen(x) && search(x.buffer,c) <= 0
        start_reading(x)
        stream_wait(x,x.readnotify)
    end
end

wait_readline(x) = wait_readbyte(x, uint8('\n'))

function wait_readnb(x::AsyncStream, nb::Int)
    while isopen(x) && nb_available(x.buffer) < nb
        start_reading(x)
        stream_wait(x,x.readnotify)
    end
end

wait_close(x) = if isopen(x) stream_wait(x,x.closenotify); end

#from `connect`
function _uv_hook_connectcb(sock::AsyncStream, status::Int32)
    @assert sock.status == StatusConnecting
    if status >= 0
        sock.status = StatusOpen
        err = ()
    else
        sock.status = StatusInit
        err = UVError("connect",status)
    end
    if isa(sock.ccb,Function)
        sock.ccb(sock, status)
    end
    err===() ? notify(sock.connectnotify) : notify_error(sock.connectnotify, err)
end

#from `listen`
function _uv_hook_connectioncb(sock::UVServer, status::Int32)
    local err
    if status >= 0
        err = ()
    else
        err = UVError("connection",status)
    end
    if isa(sock.ccb,Function)
        sock.ccb(sock,status)
    end
    err===() ? notify(sock.connectnotify) : notify_error(sock.connectnotify, err)
end

## BUFFER ##
## Allocate a simple buffer
function alloc_request(buffer::IOBuffer, recommended_size::Uint)
    ensureroom(buffer, int(recommended_size))
    ptr = buffer.append ? buffer.size + 1 : buffer.ptr
    return (pointer(buffer.data, ptr), length(buffer.data)-ptr+1)
end
function _uv_hook_alloc_buf(stream::AsyncStream, recommended_size::Uint)
    (buf,size) = alloc_request(stream.buffer, recommended_size)
    @assert size>0 # because libuv requires this (TODO: possibly stop reading too if it fails)
    (buf,uint(size))
end

function notify_filled(buffer::IOBuffer, nread::Int, base::Ptr{Void}, len::Uint)
    if buffer.append
        buffer.size += nread
    else
        buffer.ptr += nread
    end
end
function notify_filled(stream::AsyncStream, nread::Int)
    more = true
    while more
        if isa(stream.readcb,Function)
            nreadable = (stream.line_buffered ? int(search(stream.buffer, '\n')) : nb_available(stream.buffer))
            if nreadable > 0
                more = stream.readcb(stream, nreadable)
            else
                more = false
            end
        else
            more = false
        end
    end
end

const READ_BUFFER_SZ=10485760           # 10 MB
function _uv_hook_readcb(stream::AsyncStream, nread::Int, base::Ptr{Void}, len::Uint)
    if nread < 0
        if nread != UV_EOF
            # This is a fatal connectin error. Shutdown requests as per the usual
            # close function won't work and libuv will fail with an assertion failre
            ccall(:jl_forceclose_uv,Void,(Ptr{Void},),stream.handle)
            notify_error(stream.readnotify, UVError("readcb",nread))
        else
            if isa(stream,TTY)
                notify_error(stream.readnotify, EOFError())
            else
                close(stream)
                notify(stream.readnotify)
            end
        end
    else
        notify_filled(stream.buffer, nread, base, len)
        notify_filled(stream, nread)
        notify(stream.readnotify)
    end

    # Stop reading when
    # 1) when we have an infinite buffer, and we have accumulated a lot of unread data OR
    # 2) we have an alternate buffer that has reached its limit.
    if (is_maxsize_unlimited(stream.buffer) && (nb_available(stream.buffer) > READ_BUFFER_SZ )) ||
       (nb_available(stream.buffer) == stream.buffer.maxsize)
        stop_reading(stream)
    end
 end
##########################################
# Async Workers
##########################################

abstract AsyncWork

type SingleAsyncWork <: AsyncWork
    handle::Ptr{Void}
    cb::Function
    function SingleAsyncWork(cb::Function)
        this = new(c_malloc(_sizeof_uv_async), cb)
        associate_julia_struct(this.handle, this)
        preserve_handle(this)
        err = ccall(:uv_async_init,Cint,(Ptr{Void},Ptr{Void},Ptr{Void}),eventloop(),this.handle,uv_jl_asynccb::Ptr{Void})
        this
    end
end

type Timer <: AsyncWork
    handle::Ptr{Void}
    cb::Function
    function Timer(cb::Function)
        this = new(c_malloc(_sizeof_uv_timer), cb)
        # We don't want to set a julia struct, but we also
        # want to make sure there's no garbage data in the
        # ->data field
        disassociate_julia_struct(this.handle)
        err = ccall(:uv_timer_init,Cint,(Ptr{Void},Ptr{Void}),eventloop(),this.handle)
        if err != 0
            c_free(this.handle)
            this.handle = C_NULL
            error(UVError("uv_make_timer",err))
        end
        finalizer(this,uvfinalize)
        this
    end
end

close(t::Timer) = ccall(:jl_close_uv,Void,(Ptr{Void},),t.handle)

function _uv_hook_close(uv::Union(AsyncStream,UVServer))
    uv.handle = 0
    uv.status = StatusClosed
    if isa(uv.closecb, Function)
        uv.closecb(uv)
    end
    notify(uv.closenotify)
    try notify(uv.readnotify) end
    try notify(uv.connectnotify) end
end
_uv_hook_close(uv::AsyncWork) = (uv.handle = C_NULL; nothing)

# This serves as a common callback for all async classes
function _uv_hook_asynccb(async::AsyncWork, status::Int32)
    if isa(async, Timer)
        if ccall(:uv_timer_get_repeat, Uint64, (Ptr{Void},), async.handle) == 0
            # timer is stopped now
            disassociate_julia_struct(async.handle)
        end
    end
    try
        async.cb(async, status)
    catch
    end
    nothing
end

function start_timer(timer::Timer, timeout::Real, repeat::Real)
    associate_julia_struct(timer.handle, timer)
    preserve_handle(timer)
    ccall(:uv_update_time,Void,(Ptr{Void},),eventloop())
    ccall(:uv_timer_start,Cint,(Ptr{Void},Ptr{Void},Uint64,Uint64),
          timer.handle, uv_jl_asynccb::Ptr{Void}, uint64(round(timeout*1000))+1, uint64(round(repeat*1000)))
end

function stop_timer(timer::Timer)
    # ignore multiple calls to stop_timer
    !haskey(uvhandles, timer) && return
    timer.handle == C_NULL && return

    ccall(:uv_timer_stop,Cint,(Ptr{Void},),timer.handle)
    disassociate_julia_struct(timer.handle)
    unpreserve_handle(timer)
end

function sleep(sec::Real)
    w = Condition()
    timer = Timer(function (tmr,status)
        if status == 0
            notify(w)
        else
            notify_error(UVError("timer",status))
        end
    end)
    start_timer(timer, float(sec), 0)
    try
        stream_wait(timer,w)
    finally
        stop_timer(timer)
    end
    nothing
end

## event loop ##
eventloop() = global uv_eventloop::Ptr{Void}
#mkNewEventLoop() = ccall(:jl_new_event_loop,Ptr{Void},()) # this would probably be fine, but is nowhere supported

function run_event_loop()
    ccall(:jl_run_event_loop,Void,(Ptr{Void},),eventloop())
end
function process_events(block::Bool)
    loop = eventloop()
    if block
        ccall(:jl_run_once,Int32,(Ptr{Void},),loop)
    else
        ccall(:jl_process_events,Int32,(Ptr{Void},),loop)
    end
end

## pipe functions ##
function malloc_julia_pipe(x)
    x.handle = c_malloc(_sizeof_uv_named_pipe)
    associate_julia_struct(x.handle,x)
    finalizer(x,uvfinalize)
end

_link_pipe(read_end::Ptr{Void},write_end::Ptr{Void}) = uv_error("pipe_link",ccall(:uv_pipe_link, Int32, (Ptr{Void}, Ptr{Void}), read_end, write_end))

function link_pipe(read_end::Ptr{Void},readable_julia_only::Bool,write_end::Ptr{Void},writable_julia_only::Bool,readpipe::AsyncStream,writepipe::AsyncStream)
    #make the pipe an unbuffered stream for now
    #TODO: this is probably not freeing memory properly after errors
    uv_error("init_pipe",ccall(:jl_init_pipe, Cint, (Ptr{Void},Int32,Int32,Int32), read_end, 0, 1, readable_julia_only))
    uv_error("init_pipe(2)",ccall(:jl_init_pipe, Cint, (Ptr{Void},Int32,Int32,Int32), write_end, 1, 0, writable_julia_only))
    _link_pipe(read_end,write_end)
end

function link_pipe(read_end::Ptr{Void},readable_julia_only::Bool,write_end::Ptr{Void},writable_julia_only::Bool)
    uv_error("init_pipe",ccall(:jl_init_pipe, Cint, (Ptr{Void},Int32,Int32,Int32), read_end, 0, 1, readable_julia_only))
    uv_error("init_pipe(2)",ccall(:jl_init_pipe, Cint, (Ptr{Void},Int32,Int32,Int32), write_end, 1, 0, writable_julia_only))
    _link_pipe(read_end,write_end)
end

function link_pipe(read_end::Pipe,readable_julia_only::Bool,write_end::Ptr{Void},writable_julia_only::Bool)
    if read_end.handle == C_NULL
        malloc_julia_pipe(read_end)
    end
    init_pipe!(read_end; readable = true, writable = false, julia_only = readable_julia_only)
    uv_error("init_pipe",ccall(:jl_init_pipe, Cint, (Ptr{Void},Int32,Int32,Int32), write_end, 1, 0, writable_julia_only))
    _link_pipe(read_end.handle,write_end)
    read_end.status = StatusOpen
end
function link_pipe(read_end::Ptr{Void},readable_julia_only::Bool,write_end::Pipe,writable_julia_only::Bool)
    if write_end.handle == C_NULL
        malloc_julia_pipe(write_end)
    end
    uv_error("init_pipe",ccall(:jl_init_pipe, Cint, (Ptr{Void},Int32,Int32,Int32), read_end, 0, 1, readable_julia_only))
    init_pipe!(write_end; readable = false, writable = true, julia_only = writable_julia_only)
    _link_pipe(read_end,write_end.handle)
    write_end.status = StatusOpen
end
function link_pipe(read_end::Pipe,readable_julia_only::Bool,write_end::Pipe,writable_julia_only::Bool)
    if write_end.handle == C_NULL
        malloc_julia_pipe(write_end)
    end
    if read_end.handle == C_NULL
        malloc_julia_pipe(read_end)
    end
    init_pipe!(read_end; readable = true, writable = false, julia_only = readable_julia_only)
    init_pipe!(write_end; readable = false, writable = true, julia_only = writable_julia_only)
    _link_pipe(read_end.handle,write_end.handle)
    write_end.status = StatusOpen
    read_end.status = StatusOpen
    nothing
end
close_pipe_sync(p::Pipe) = (ccall(:uv_pipe_close_sync,Void,(Ptr{Void},),p.handle); p.status = StatusClosed)
close_pipe_sync(handle::UVHandle) = ccall(:uv_pipe_close_sync,Void,(UVHandle,),handle)

_uv_hook_isopen(stream) = int32(stream.status != StatusUninit && stream.status != StatusInit && isopen(stream))

function close(stream::Union(AsyncStream,UVServer))
    if isopen(stream) && stream.status != StatusClosing
        ccall(:jl_close_uv,Void,(Ptr{Void},),stream.handle)
        stream.status = StatusClosing
    end
    nothing
end

## stream functions ##
function start_reading(stream::AsyncStream)
    if stream.status == StatusOpen
        if !isreadable(stream)
            error("tried to read a stream that is not readable")
        end
        ret = ccall(:uv_read_start,Cint,(Ptr{Void},Ptr{Void},Ptr{Void}),
                    handle(stream),uv_jl_alloc_buf::Ptr{Void},uv_jl_readcb::Ptr{Void})
        stream.status = StatusActive
        ret
    elseif stream.status == StatusActive
        int32(0)
    else
        int32(-1)
    end
end
function start_reading(stream::AsyncStream, cb::Function)
    start_reading(stream)
    stream.readcb = cb
    nread = nb_available(stream.buffer)
    if nread > 0
        notify_filled(stream, nread)
    end
    nothing
end
start_reading(stream::AsyncStream, cb::Bool) = (start_reading(stream); stream.readcb = cb; nothing)

function stop_reading(stream::AsyncStream)
    if stream.status == StatusActive
        ret = ccall(:uv_read_stop,Cint,(Ptr{Void},),stream.handle)
        stream.status = StatusOpen
        ret
    elseif stream.status == StatusOpen
        int32(0)
    else
        int32(-1)
    end
end

function readall(stream::AsyncStream)
    start_reading(stream)
    wait_close(stream)
    return takebuf_string(stream.buffer)
end

function read!{T}(s::AsyncStream, a::Array{T})
    isbits(T) || error("read from buffer only supports bits types or arrays of bits types")
    nb = length(a) * sizeof(T)
    read!(s, reshape(reinterpret(Uint8, a), nb))
    return a
end

function read!{Uint8}(s::AsyncStream, a::Vector{Uint8})
    nb = length(a)
    sbuf = s.buffer
    @assert sbuf.seekable == false
    @assert sbuf.maxsize >= nb

    if nb_available(sbuf) >= nb
        return read!(sbuf, a)
    end

    if nb <= 65536 # Arbitrary 64K limit under which we are OK with copying the array from the stream's buffer
        wait_readnb(s,nb)
        read!(sbuf, a)
    else
        stop_reading(s) # Just playing it safe, since we are going to switch buffers.
        newbuf = PipeBuffer(a, nb)
        newbuf.size = 0
        s.buffer = newbuf
        write(newbuf, sbuf)
        wait_readnb(s,nb)
        s.buffer = sbuf
    end
    return a
end

function read{T}(s::AsyncStream, ::Type{T}, dims::Dims)
    isbits(T) || error("read from buffer only supports bits types or arrays of bits types")
    nb = prod(dims)*sizeof(T)
    a = read!(s, Array(Uint8, nb))
    reshape(reinterpret(T, a), dims)
end

function read(this::AsyncStream,::Type{Uint8})
    buf = this.buffer
    @assert buf.seekable == false
    wait_readnb(this,1)
    read(buf,Uint8)
end

function readline(this::AsyncStream)
    buf = this.buffer
    @assert buf.seekable == false
    wait_readline(this)
    readline(buf)
end

readline() = readline(STDIN)

function readavailable(this::AsyncStream)
    buf = this.buffer
    @assert buf.seekable == false
    wait_readnb(this,1)
    takebuf_string(buf)
end

function readuntil(this::AsyncStream,c::Uint8)
    buf = this.buffer
    @assert buf.seekable == false
    wait_readbyte(this,c)
    readuntil(buf,c)
end

#function finish_read(pipe::Pipe)
#    close(pipe) #handles to UV and ios will be invalid after this point
#end
#
#function finish_read(state::(Pipe,ByteString))
#    finish_read(state...)
#end

macro uv_write(n,call)
    esc(quote
        check_open(s)
        uvw = c_malloc(_sizeof_uv_write+$(n))
        err = $call
        if err < 0
            c_free(uvw)
            uv_error("write", err)
        end
    end)
end

## low-level calls ##

function write!{T}(s::AsyncStream, a::Array{T})
    if isbits(T)
        n = uint(length(a)*sizeof(T))
        @uv_write n ccall(:jl_write_no_copy, Int32, (Ptr{Void}, Ptr{Void}, Uint, Ptr{Void}, Ptr{Void}), handle(s), a, n, uvw, uv_jl_writecb::Ptr{Void})
        return int(length(a)*sizeof(T))
    else
        throw(MethodError(write,(s,a)))
    end
end
function write!(s::AsyncStream, p::Ptr, nb::Integer)
    @uv_write nb ccall(:jl_write_no_copy, Int32, (Ptr{Void}, Ptr{Void}, Uint, Ptr{Void}, Ptr{Void}), handle(s), p, nb, uvw, uv_jl_writecb::Ptr{Void})
    return nb
end
write!(s::AsyncStream, string::ByteString) = write!(s,string.data)

function _uv_hook_writecb(s::AsyncStream, req::Ptr{Void}, status::Int32)
    status < 0 && close(s)
    nothing
end

function write(s::AsyncStream, b::Uint8)
    @uv_write 1 ccall(:jl_putc_copy, Int32, (Uint8, Ptr{Void}, Ptr{Void}, Ptr{Void}), b, handle(s), uvw, uv_jl_writecb_task::Ptr{Void})
    ct = current_task()
    uv_req_set_data(uvw,ct)
    ct.state = :waiting
    stream_wait(ct)
    return 1
end
function write(s::AsyncStream, c::Char)
    @uv_write utf8sizeof(c) ccall(:jl_pututf8_copy, Int32, (Ptr{Void},Uint32, Ptr{Void}, Ptr{Void}), handle(s), c, uvw, uv_jl_writecb_task::Ptr{Void})
    ct = current_task()
    uv_req_set_data(uvw,ct)
    ct.state = :waiting
    stream_wait(ct)
    return utf8sizeof(c)
end
function write{T}(s::AsyncStream, a::Array{T})
    if isbits(T)
        n = uint(length(a)*sizeof(T))
        @uv_write n ccall(:jl_write_no_copy, Int32, (Ptr{Void}, Ptr{Void}, Uint, Ptr{Void}, Ptr{Void}), handle(s), a, n, uvw, uv_jl_writecb_task::Ptr{Void})
        ct = current_task()
        uv_req_set_data(uvw,ct)
        ct.state = :waiting
        stream_wait(ct)
        return int(length(a)*sizeof(T))
    else
        check_open(s)
        invoke(write,(IO,Array),s,a)
    end
end
function write(s::AsyncStream, p::Ptr, nb::Integer)
    @uv_write nb ccall(:jl_write_no_copy, Int32, (Ptr{Void}, Ptr{Void}, Uint, Ptr{Void}, Ptr{Void}), handle(s), p, nb, uvw, uv_jl_writecb_task::Ptr{Void})
    ct = current_task()
    uv_req_set_data(uvw,ct)
    ct.state = :waiting
    stream_wait(ct)
    return int(nb)
end

function _uv_hook_writecb_task(s::AsyncStream,req::Ptr{Void},status::Int32)
    d = uv_req_data(req)
    if status < 0
        err = UVError("write",status)
        close(s)
        if d != C_NULL
            schedule(unsafe_pointer_to_objref(d)::Task,err,error=true)
        end
    elseif d != C_NULL
        schedule(unsafe_pointer_to_objref(d)::Task)
    end
end

## Libuv error handling ##
type UVError <: Exception
    prefix::String
    code::Int32
    UVError(p::String,code::Integer)=new(p,int32(code))
end

struverror(err::UVError) = bytestring(ccall(:uv_strerror,Ptr{Uint8},(Int32,),err.code))
uverrorname(err::UVError) = bytestring(ccall(:uv_err_name,Ptr{Uint8},(Int32,),err.code))

uv_error(prefix::Symbol, c::Integer) = uv_error(string(prefix),c)
uv_error(prefix::String, c::Integer) = c < 0 ? throw(UVError(prefix,c)) : nothing
show(io::IO, e::UVError) = print(io, e.prefix*": "*struverror(e)*" ("*uverrorname(e)*")")

## server functions ##

function accept_nonblock(server::PipeServer,client::Pipe)
    if client.status != StatusInit
        error(client.status == StatusUninit ? "client is not initialized" :
              "client is already in use or has been closed")
    end
    err = ccall(:uv_accept,Int32,(Ptr{Void},Ptr{Void}),server.handle,client.handle)
    if err == 0
        client.status = StatusOpen
    end
    err
end
function accept_nonblock(server::PipeServer)
    client = Pipe()
    uv_error("accept", accept_nonblock(server,client) != 0)
    client
end

function accept(server::UVServer, client::AsyncStream)
    if server.status != StatusActive
        error("server not connected; make sure \"listen\" has been called")
    end
    while isopen(server)
        err = accept_nonblock(server,client)
        if err == 0
            return client
        elseif err != UV_EAGAIN
            uv_error("accept",err)
        end
        stream_wait(server,server.connectnotify)
    end
    error("server was closed while attempting to accept a client")
end

const BACKLOG_DEFAULT = 511

function _listen(sock::UVServer; backlog::Integer=BACKLOG_DEFAULT)
    err = ccall(:uv_listen, Cint, (Ptr{Void}, Cint, Ptr{Void}),
                sock.handle, backlog, uv_jl_connectioncb::Ptr{Void})
    sock.status = StatusActive
    err
end

function bind(server::PipeServer, name::ASCIIString)
    @assert server.status == StatusInit
    err = ccall(:uv_pipe_bind, Int32, (Ptr{Void}, Ptr{Uint8}),
                server.handle, name)
    if err != 0
        if err != UV_EADDRINUSE && err != UV_EACCES
            error(UVError("bind",err))
        else
            return false
        end
    end
    server.status = StatusOpen
    true
end

function listen(path::ByteString)
    sock = PipeServer()
    bind(sock, path) || error("could not listen on path $path")
    uv_error("listen", _listen(sock))
    sock
end

function connect!(sock::Pipe, path::ByteString)
    @assert sock.status == StatusInit
    req = c_malloc(_sizeof_uv_connect)
    uv_req_set_data(req,C_NULL)
    ccall(:uv_pipe_connect, Void, (Ptr{Void}, Ptr{Void}, Ptr{Uint8}, Ptr{Void}), req, sock.handle, path, uv_jl_connectcb::Ptr{Void})
    sock.status = StatusConnecting
    sock
end
connect!(sock::Pipe, path::String) = connect(sock,bytestring(path))

function connect(sock::AsyncStream, args...)
    connect!(sock,args...)
    wait_connected(sock)
    sock
end

connect(path::String) = connect(Pipe(),path)

dup(x::RawFD) = RawFD(ccall((@windows? :_dup : :dup),Int32,(Int32,),x.fd))
dup(src::RawFD,target::RawFD) = systemerror("dup",-1==
    ccall((@windows? :_dup2 : :dup2),Int32,
    (Int32,Int32),src.fd,target.fd))

_fd(x::IOStream) = RawFD(fd(x))
@unix_only _fd(x::AsyncStream) = RawFD(ccall(:jl_uv_handle,Int32,(Ptr{Void},),x.handle))
@windows_only _fd(x::AsyncStream) = WindowsRawSocket(
    ccall(:jl_uv_handle,Ptr{Void},(Ptr{Void},),x.handle))

for (x,writable,unix_fd,c_symbol) in ((:STDIN,false,0,:jl_uv_stdin),(:STDOUT,true,1,:jl_uv_stdout),(:STDERR,true,2,:jl_uv_stderr))
    f = symbol("redirect_"*lowercase(string(x)))
    _f = symbol(string("_",f))
    @eval begin
        function ($_f)(stream)
            global $x
            @windows? (
                ccall(:SetStdHandle,stdcall,Int32,(Uint32,Ptr{Void}),
                    $(-10-unix_fd),_get_osfhandle(_fd(stream)).handle) :
                dup(_fd(stream),  RawFD($unix_fd)) )
            $x = stream
        end
        function ($f)(handle::AsyncStream)
            $(_f)(handle)
            unsafe_store!(cglobal($(Expr(:quote,c_symbol)),Ptr{Void}),
                handle.handle)
            handle
        end
        ($f)(handle::IOStream) = ($_f)(handle)
        function ($f)()
            read,write = (Pipe(C_NULL), Pipe(C_NULL))
            link_pipe(read,$(writable),write,$(!writable))
            ($f)($(writable? :write : :read))
            (read,write)
        end
    end
end

max_stars_repo_name: aplund/julia

max_issues_repo_name: aplund/julia

max_forks_repo_name: aplund/julia


text:
#Module for D4M
module D4M

    #=
        Loaded Example
    =#
    module_dir = Pkg.dir("D4M")
    ParallelDatabase_dir = module_dir*"/examples/3Scaling/2ParallelDatabase"
    example_3Scaling_2ParallelDatabase_pDB01_DataTEST = () -> include( ParallelDatabase_dir *"/pDB01_DataTEST.jl")
    example_3Scaling_2ParallelDatabase_pDB02_FileTEST = () -> include( ParallelDatabase_dir *"/pDB02_FileTEST.jl")
    example_3Scaling_2ParallelDatabase_pDB03_AssocTEST = () -> include( ParallelDatabase_dir *"/pDB03_AssocTEST.jl")
    example_3Scaling_2ParallelDatabase__pDB04_DegreeTEST = () -> include( ParallelDatabase_dir *"/pDB04_DegreeTEST.jl")

    AssocIntro_dir = module_dir*"/examples/1Intro/1AssocIntro"
    example_1Intro_1AssocIntro_AI1_SetupTEST = () -> include( AssocIntro_dir * "/AI1_SetupTEST.jl")
    example_1Intro_1AssocIntro_AI2_SubsrefTEST = () -> include( AssocIntro_dir * "/AI2_SubsrefTEST.jl")
    example_1Intro_1AssocIntro_AI3_MathTEST = () -> include( AssocIntro_dir * "/AI3_MathTEST.jl")
    example_1Intro_1AssocIntro_AI4_AdvConstructTEST = () -> include( AssocIntro_dir * "/AI4_AdvConstructTEST.jl")

    EdgeArt_dir = module_dir*"/examples/1Intro/2EdgeArt"
    example_1Intro_2EdgeArt_EA1_GraphTEST = () -> include( EdgeArt_dir * "/EA1_GraphTEST.jl")
    example_1Intro_2EdgeArt_EA2_SubsrefTEST = () -> include( EdgeArt_dir * "/EA2_SubsrefTEST.jl")
    example_1Intro_2EdgeArt_EA3_SubGraphTEST = () -> include( EdgeArt_dir * "/EA3_SubGraphTEST.jl")

    EntityAnalysis_dir = module_dir*"/examples/2Apps/1EntityAnalysis"
    example_2Apps_1EntityAnalysis_EA1_GraphTEST = () -> include( EntityAnalysis_dir * "/EA1_ReadTEST.jl")
    example_2Apps_1EntityAnalysis_EA2_SubsrefTEST = () -> include( EntityAnalysis_dir * "/EA2_StatTEST.jl")
    example_2Apps_1EntityAnalysis_EA3_FacetTEST = () -> include( EntityAnalysis_dir * "/EA3_FacetTEST.jl")
    example_2Apps_1EntityAnalysis_EA4_GraphTEST = () -> include( EntityAnalysis_dir * "/EA4_GraphTEST.jl")
    example_2Apps_1EntityAnalysis_EA5_GraphQueryTEST = () -> include( EntityAnalysis_dir * "/EA5_GraphQueryTEST.jl")

    MatrixPerformance_dir = module_dir*"/examples/3Scaling/3MatrixPerformance"
    example_3Scaling_3MatrixPerformance_MP1_DenseTEST = () -> include( MatrixPerformance_dir * "/MP1_DenseTEST.jl")
    example_3Scaling_3MatrixPerformance_MP2_SparseTEST = () -> include( MatrixPerformance_dir * "/MP2_SparseTEST.jl")
    example_3Scaling_3MatrixPerformance_MP3_AssocTEST = () -> include( MatrixPerformance_dir * "/MP3_AssocTEST.jl")
    example_3Scaling_3MatrixPerformance_MP4_AssocCatKeyTEST = () -> include( MatrixPerformance_dir * "/MP4_AssocCatKeyTEST.jl")
    example_3Scaling_3MatrixPerformance_MP5_AssocCatValKeyTEST = () -> include( MatrixPerformance_dir * "/MP5_AssocCatValKeyTEST.jl")
    example_3Scaling_3MatrixPerformance_MP6_AssocPlusTEST = () -> include( MatrixPerformance_dir * "/MP6_AssocPlusTEST.jl")

end

include("Assoc_orig.jl") # Associative Array
include("WriteCSV.jl") #load Assoc from CSV file
include("ReadCSV.jl") #Dump Assoc into a CSV file
include("OutDegree.jl") #Calculate OutDegree
include("NumStr.jl") #Quickly calculate the length of string sequence separated by single-character
include("CatKeyMul.jl") #CatKeyMultiply
include("CatValMul.jl") #CatValMultiply
include("CatStr.jl")    #Cat String
include("SplitStr.jl")
include("col2type.jl")

#Operation that gain performance benefit on sorted inputs
include("sortedintersect.jl")
include("sortedintersectmapping.jl")
include("sortedunion.jl")
include("searchsortedmapping.jl")

########################################################
# D4M: Dynamic Distributed Dimensional Data Model
# Architect: Dr. Jeremy Kepner (kepner@ll.mit.edu)
# Software Engineer: Alexander Chen (alexc89@mit.edu)
########################################################

ext: jl


text:
@testset "Outer product" begin
    t1 = Tensor([1, 2, 4])
    t2 = Tensor([1.0, 0.5, 0.25])

    # the outer product of two 1-dimensional tensors yields a 2-dimensional tensor
    t3 = t1 ⊗ t2
    @test eltype(t3) === Float64
    @test size(t3) == (size(t1)..., size(t2)...)
    @test t3 == Tensor([
        1.0 0.5 0.25
        2.0 1.0 0.5
        4.0 2.0 1.0
    ])

    t4 = Tensor(rand(3, 3))
    @test t4 ⊗ Tensor(3.14) == t4 * 3.14
    @test t4 ⊗ Tensor([1, 1]) == Tensor(cat(t4.data, t4.data; dims=3))
    @test Tensor(rand(4, 5, 6)) ⊗ Tensor(zeros(3, 2)) == Tensor(zeros(4, 5, 6, 3, 2))

    # you can multiply multiple tensors at once
    t5 = ⊗(Tensor(rand(3, 6)), Tensor(rand(2, 8, 2)), Tensor(rand(9, 1, 7)))
    @test size(t5) == (3, 6, 2, 8, 2, 9, 1, 7)
end

@testset "Inner product" begin
    # for 1-dimensional tensors the inner product is equivalent to the vector dot product
    t1 = Tensor([1, 2, 3])
    t2 = Tensor([5, 2, -3])
    @test t1 ⋅ t2 == Tensor(0)

    t3 = Tensor([
        1 -1 1
        -1 1 -1
        1 -1 1
    ])
    @test t1 ⋅ t3 == Tensor([2, -2, 2])
    @test t2 ⋅ t3 == Tensor([0, 0, 0])

    # for 2-dimensional tensors the inner product is equivalent to the matrix product
    t4 = Tensor(rand(3, 6))
    t5 = Tensor(rand(6, 5))
    @test t4 ⋅ t5 ≈ Tensor(t4.data * t5.data)

    # analogously for a matrix-vector product
    t6 = Tensor(rand(6))
    @test t4 ⋅ t6 ≈ Tensor(t4.data * t6.data)
end

@testset "Contraction" begin
    a1 = reshape(Vector(1:12), 2, 2, 3)
    t1 = Tensor(a1)
    @test contract(t1, 1, 2) == Tensor([5, 13, 21])
    @test contract(t1, 1, 2) == contract(t1, 2, 1)

    a2 = rand(6, 3, 6, 2)
    t2 = Tensor(a2)
    t3 = contract(t2, 1, 3)
    @test size(t3) == (3, 2)
    @test t3 == Tensor(
        [
            sum(a2[i, 1, i, 1] for i in axes(a2, 1)) sum(a2[i, 1, i, 2] for i in axes(a2, 1))
            sum(a2[i, 2, i, 1] for i in axes(a2, 1)) sum(a2[i, 2, i, 2] for i in axes(a2, 1))
            sum(a2[i, 3, i, 1] for i in axes(a2, 1)) sum(a2[i, 3, i, 2] for i in axes(a2, 1))
        ],
    )

    # contracting is only possible along distinct, valid dimensions of the same size
    @test_throws ArgumentError contract(t1, 2, 2) # not distinct
    @test_throws BoundsError contract(t1, 2, 4) # 4 not a valid dimension for rank 3 tensor
    @test_throws DimensionMismatch contract(t1, 2, 3) # dimensions 2 and 3 are not equal

    # for rank 2 tensors we can use the trace method instead of contract
    @test trace(t1[:, :, 3]) == Tensor(21)
    @test trace(t1[2, :, 1:2]) == Tensor(10)
    @test_throws DimensionMismatch trace(t1[2, :, :])
    @test trace(t2[:, 2, :, 1]) == contract(t2[:, 2, :, 1], 1, 2)
    @test trace(t2[:, 2, :, 1]) == contract(t2[:, 2, :, 1], 2, 1)
end

@testset "Überschiebung" begin
    t1 = Tensor(reshape(Vector(1:12), 2, 3, 2))
    t2 = Tensor([
        1.1 2.2 3.3
        -3.3 -2.2 -1.1
    ])

    t3 = pushover(t1, t2, 2, 2)
    @test size(t3) == (2, 2, 2)
    @test t3 ≈ Tensor(reshape([24.2, 30.8, 63.8, 70.4, -15.4, -22, -55, -61.6], 2, 2, 2))

    t4 = pushover(t1, t2, 1, 1)
    @test size(t4) == (3, 2, 3)
    @test t4 ≈ Tensor(
        reshape(
            [
                -5.5,
                -9.9,
                -14.3,
                -18.7,
                -23.1,
                -27.5,
                -2.2,
                -2.2,
                -2.2,
                -2.2,
                -2.2,
                -2.2,
                1.1,
                5.5,
                9.9,
                14.3,
                18.7,
                23.1,
            ],
            3,
            2,
            3,
        ),
    )

    t5 = pushover(t1, t2, 3, 1)
    @test size(t5) == (2, 3, 3)
    @test t5 ≈ Tensor(
        reshape(
            [
                -22,
                -24.2,
                -26.4,
                -28.6,
                -30.8,
                -33,
                -13.2,
                -13.2,
                -13.2,
                -13.2,
                -13.2,
                -13.2,
                -4.4,
                -2.2,
                0,
                2.2,
                4.4,
                6.6,
            ],
            2,
            3,
            3,
        ),
    )

    # push-over is only possible along valid dimensions of the same size
    @test_throws BoundsError pushover(t1, t2, 3, 3)
    @test_throws DimensionMismatch pushover(t1, t2, 1, 2)

    # push-over satifies many interesting tensor-algebraic identities
    t6 = Tensor(rand(4, 7, 5))
    t7 = Tensor(rand(6, 3, 4))
    @test pushover(t6, t7, 1, 3) == permutedims(pushover(t7, t6, 3, 1), (3, 4, 1, 2))
    @test pushover(t6, t7, 1, 3) == contract(t6 ⊗ t7, 1, 6)
    @test pushover(t6, t7, 1, 3) == permutedims(t6, (2, 3, 1)) ⋅ permutedims(t7, (3, 1, 2))
end

# TODO: Test cross product (a.k.a. vector product)


text:
"""
`randomSDR()` creates a dataset of N SDRs with sparsity between 2-20%
"""
function randomSDR(n; sparsity=nothing)
  len= 32*32
  data= falses(len,n)
  for i= 1:n
    if isnothing(sparsity)
      k= floor(Int, (rand()*0.18 + 0.02) * len)  # number of 1s
    else
      k= floor(Int, sparsity * len)  # number of 1s
    end
    @views data[ randperm(len)[1:k] ,i].= true
  end
  data
end

max_stars_repo_name: Oblynx/HTMexperiments

max_issues_repo_name: Oblynx/HTMexperiments

max_forks_repo_name: Oblynx/HTMexperiments


text:
module Interact

using Reexport

@reexport using InteractBase
import InteractBase: notifications
import Widgets: Widget, @layout, @nodeps
import Observables: @on, @map!, @map

@reexport using OrderedCollections
@reexport using Observables
@reexport using Knockout
@reexport using CSSUtil
@reexport using WebIO
@reexport using Widgets

struct Bulma<:InteractBase.WidgetTheme; end

const notebookdir = joinpath(@__DIR__, "..", "doc", "notebooks")

const bulma_css = joinpath(@__DIR__, "..", "assets", "bulma.min.css")
const bulma_confined_css = joinpath(@__DIR__, "..", "assets", "bulma_confined.min.css")

function InteractBase.libraries(::Bulma)
    bulmalib = InteractBase.isijulia() ? bulma_confined_css : bulma_css
    vcat(InteractBase.font_awesome, InteractBase.style_css, bulmalib)
end

const themes = Dict(
    :nativehtml => InteractBase.NativeHTML(),
    :bulma => Bulma()
)

function InteractBase.settheme!(s::Symbol)
    (s in keys(themes)) || error("Theme $s is not supported")
    settheme!(themes[s])
end

function __init__()
    settheme!(:bulma)
    nothing
end

end


text:
using Documenter
using GeometricIntegrators

makedocs(;
    authors="Michael Kraus",
    sitename="SPARK Methods for Degenerate Lagrangian Systems",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://michakraus.github.io/papers-spark-methods-for-degenerate-lagrangian-systems",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Integrators" => "integrators.md",
        "Tableaus" => "tableaus.md",

        "Lotka-Volterra 2d" => [
            "Gauss-Legendre VPRK Methods" => "lotka-volterra-2d/lotka-volterra-2d-spark-glvprk.md",
            "Gauss-Legendre SPARK Methods" => "lotka-volterra-2d/lotka-volterra-2d-spark-glspark.md",
            "Gauss-Lobatto SPARK Methods" => "lotka-volterra-2d/lotka-volterra-2d-spark-lobspark.md",
            "Internal Stage Projection"  => "lotka-volterra-2d/lotka-volterra-2d-spark-pinternal.md",
            "Lobatto-IIIA-IIIB Projection"  => "lotka-volterra-2d/lotka-volterra-2d-spark-plobatto.md",
            "Midpoint Projection"  => "lotka-volterra-2d/lotka-volterra-2d-spark-pmidpoint.md",
            "Symmetric Projection" => "lotka-volterra-2d/lotka-volterra-2d-spark-psymmetric.md",
        ],

        "Massless Charged Particle" => [
            "Gauss-Legendre VPRK Methods" => "massless-charged-particle/massless-charged-particle-spark-glvprk.md",
            "Gauss-Legendre SPARK Methods" => "massless-charged-particle/massless-charged-particle-spark-glspark.md",
            "Gauss-Lobatto SPARK Methods" => "massless-charged-particle/massless-charged-particle-spark-lobspark.md",
            "Internal Stage Projection"  => "massless-charged-particle/massless-charged-particle-spark-pinternal.md",
            "Lobatto-IIIA-IIIB Projection"  => "massless-charged-particle/massless-charged-particle-spark-plobatto.md",
            "Midpoint Projection"  => "massless-charged-particle/massless-charged-particle-spark-pmidpoint.md",
            "Symmetric Projection" => "massless-charged-particle/massless-charged-particle-spark-psymmetric.md",
        ],
    ],
)

deploydocs(;
    repo="github.com/michakraus/paper-spark-methods-for-degenerate-lagrangian-systems",
    devbranch="main"
)


text:
using SurveyAnalysis
using Test

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


text:
const VERSION_DIR_NAME = ".versions"
const TEMP_CALC_DIR = "outputs"

@enum JobState BootFail Pending Running Completed Cancelled Deadline Failed NodeFail OutOfMemory Preempted Requeued Resizing Revoked Suspended Timeout Submitted Unknown

"""
    Job(name::String, structure::Structure;
          calculations      ::Vector{Calculation} = Calculation[],
          dir               ::String = pwd(),
          header            ::Vector{String} = getdefault_jobheader(),
          metadata          ::Dict = Dict(),
          version           ::Int = last_job_version(dir),
          copy_temp_folders ::Bool = false,
          server            ::String = getdefault_server(),
          environment::String ="")

A [`Job`](@ref) embodies a set of [`Calculations`](@ref Calculation) to be ran in directory `dir`, with the [`Structure`](@ref) as the subject.
## Keywords/further attributes
- `calculations`: calculations to calculations that will be run sequentially.
- `dir`: the directory where the calculations will be run.
- `header`: lines that will be pasted at the head of the job script, e.g. exports `export OMP_NUM_THREADS=1`, slurm settings`#SBATCH`, etc.
- `metadata`: various additional information, will be saved in `.metadata.jld2` in the `dir`.
- `version`: the current version of the job.
- `copy_temp_folders`: whether or not the temporary directory associated with intermediate calculation results should be copied when storing a job version. *CAUTION* These can be quite large.
- `server`: [`Server`](@ref) where to run the [`Job`](@ref).
- `environment`: [`Environment`](@ref) to be used for running the [`Job`](@ref).

    Job(job_name::String, structure::Structure, calculations::Vector{<:Calculation}, common_flags::Pair{Symbol, <:Any}...; kwargs...)

Creates a new job. The common flags will be attempted to be set in each of the `calculations`. The `kwargs...` are passed to the [`Job`](@ref) constructor.

    Job(job_dir::String, job_script="job.tt"; version=nothing, kwargs...)

Loads the job in the `dir`.
If `job_dir` is not a valid job path, the previously saved jobs will be scanned for a job with a `dir` that
partly includes `job_dir`. If `version` is specified the corresponding job version will be returned if it exists.
The `kwargs...` will be passed to the [`Job`](@ref) constructor.
"""
@with_kw_noshow mutable struct Job
    name::String = ""
    structure::Structure = Structure()
    calculations::Vector{Calculation} = Calculation[]
    dir::String = pwd()
    header::Vector{String} = String[]
    metadata::Dict{Symbol,Any} = Dict{Symbol,Any}()
    version::Int = -1
    copy_temp_folders::Bool = false
    server::String = "localhost"
    environment::String = ""
    function Job(name, structure, calculations, dir, header, metadata, version,
                 copy_temp_folders, server, environment)
        if dir[end] == '/'
            dir = dir[1:end-1]
        end
        out = new(name, structure, calculations, dir, header, metadata, version,
                  copy_temp_folders, server, environment)
        return out
    end
end

#TODO implement abinit
function Job(job_name::String, structure::Structure, calculations::Vector{<:Calculation},
             common_flags::Pair{Symbol,<:Any}...; kwargs...)
    out = Job(; name = job_name, structure = structure, calculations = calculations,
              kwargs...)
    for (f, v) in common_flags
        out[f] = v
    end
    return out
end

StructTypes.StructType(::Type{Job}) = StructTypes.Mutable()

#-------------------BEGINNING GENERAL SECTION-------------#
scriptpath(job::Job) = joinpath(job, "job.tt")
starttime(job::Job)  = mtime(scriptpath(job))

runslocal(job::Job) = job.server == "localhost"
isarchived(job::Job) = occursin(".archived", job.dir)

Base.abspath(job::Job) =
    isabspath(job.dir) ? job.dir : joinpath(Server(job.server), job.dir)

Base.ispath(job::Job, p...) =
    runslocal(job) ? joinpath(job, p...) : ispath(DFC.Servers.maybe_start_server(DFC.Server(job.server)), joinpath(job, p...))

"""
    joinpath(job::Job, args...)

`joinpath(job.dir, args...)`.
"""
Base.joinpath(job::Job, args...) = joinpath(abspath(job), args...)
Base.readdir(job::Job) = runslocal(job) ? readdir(abspath(job)) : readdir(Server(job.server), job.dir)

function Base.pop!(job::Job, name::String)
    i = findfirst(x -> x.name == name, job.calculations)
    if i === nothing
        error("No calculation with name $name found.")
    end
    out = job.calculations[i]
    deleteat!(job.calculations, i)
    return out
end

function Utils.searchdir(job::Job, str::AbstractString)
    if runslocal(job)
        return searchdir(abspath(job), str)
    else
        return searchdir(Server(job.server), abspath(job), str)
    end
end

function Base.setindex!(job::Job, value, key::Symbol)
    for c in job.calculations
        c[key] = value
    end
end

function Base.getindex(job::Job, flg::Symbol)
    outdict = Dict()
    for i in job.calculations
        tfl = get(i, flg, nothing)
        if tfl !== nothing
            outdict[i.name] = tfl
        end
    end
    return outdict
end

function Base.pop!(job::Job, args...)
    out = Dict{String, Any}()
    for c in job.calculations
        out[c.name] = pop!(c, args...)
    end
    return out
end

"""
    insert!(job::Job, i::Int, calculation::Calculation) = insert!(job.calculations, i, calculation)
"""
function Base.insert!(job::Job, index::Int, calculation::Calculation)
    return insert!(job.calculations, index, calculation)
end

"""
    push!(job::Job, calculation::Calculation) = push!(job.calculations, calculation)
"""
Base.push!(job::Job, calculation::Calculation) = push!(job.calculations, calculation)

"""
    pop!(job::Job) = pop!(job.calculations)
"""
Base.pop!(job::Job) = pop!(job.calculations)

"""
    append!(job::Job, args...) = append!(job.calculations, args...)
"""
Base.append!(job::Job, args...) = append!(job.calculations, args...)
Base.length(job::Job) = length(job.calculations)
Base.lastindex(job::Job) = length(job)
Base.getindex(job::Job, i::Integer) = job.calculations[i]

"""
    getindex(job::Job, name::String)

Returns the `Calculation` with the specified `name`.

    getindex(job::Job, i::Integer)

Returns the i'th `Calculation` in the job.
"""
function Base.getindex(job::Job, id::String)
    tmp = getfirst(x -> x.name == id, job.calculations)
    if tmp != nothing
        return tmp
    else
        error("Calculation $id not found.")
    end
end

Base.getindex(job::Job, el::Element) = job.structure[el]

"""
    main_job_dir(dir::AbstractString)
    main_job_dir(job::Job)

Returns the main directory of the job, also when the job's version is not the one
in the main directory.
"""
function main_job_dir(dir::AbstractString)
    d = split(dir, Jobs.VERSION_DIR_NAME)[1]
    return d[end] == '/' ? d[1:end-1] : d
end
main_job_dir(job::Job) = isabspath(job.dir) ? main_job_dir(job.dir) : main_job_dir(joinpath(Server(job), job.dir))

"""
    gencalc_wan(job::Job, min_window_determinator::Real, extra_wan_flags...; kwargs...)

Automates the generation of wannier calculations based on the `job`.
When a projwfc calculation is present in the `job`, `min_window_determinator` will be used to
determine the threshold value for including a band in the window based on the projections, otherwise
it will be used as the `Emin` value from which to start counting the number of bands needed for all
projections.
`extra_wan_flags` can be any extra flags for the Wannier90 calculation such as `write_hr` etc.
"""
function Calculations.gencalc_wan(job::Job, min_window_determinator::Real,
                                  extra_wan_flags...; kwargs...)
    nscf_calc = getfirst(x -> Calculations.isnscf(x), job.calculations)
    nscf_calc === nothing && "Please first run an nscf calculation."
    if get(nscf_calc, :nosym, false) != true
        @info "'nosym' flag was not set in the nscf calculation.
                If this was not intended please set it and rerun the nscf calculation.
                This generally gives errors because of omitted kpoints, needed for pw2wannier90.x"
    end
    projwfc_calc = getfirst(x -> Calculations.isprojwfc(x), job.calculations)
    outdat = outputdata(job)
    if projwfc_calc === nothing || !haskey(outdat, projwfc_calc.name)
        @info "No projwfc calculation found with valid output, using $min_window_determinator as Emin"
        Emin = min_window_determinator
    else
        @info "Valid projwfc output found, using $min_window_determinator as the dos threshold."
        Emin = Calculations.Emin_from_projwfc(job.structure, outdat[projwfc_calc.name][:states], outdat[projwfc_calc.name][:bands], min_window_determinator)
    end
    return Calculations.gencalc_wan(nscf_calc, job.structure, outdat[nscf_calc.name][:bands], Emin, extra_wan_flags...; kwargs...)
end

"""
    set_flow!(job::Job, should_runs::Pair{String, Bool}...)

Sets whether or not calculations should be scheduled to run.
The `name` of each calculation in the job will be checked against the string in each pair of `should_runs`, and the
`calculation.run` will be set accordingly.

Example:
```julia
set_flow!(job, "" => false, "scf" => true)
```
would un-schedule all calculations in the job, and schedule the "scf" and "nscf" calculations to run.
"""
function set_flow!(job::Job, should_runs...)
    for (name, run) in should_runs
        for calculation in filter(x -> occursin(name, x.name), job.calculations)
            calculation.run = run
        end
    end
    return job
end

"""
    set_headerword!(job::Job, old_new::Pair{String, String})

Replaces the specified word in the header with the new word.
"""
function set_headerword!(job::Job, old_new::Pair{String,String}; print = true)
    for (i, line) in enumerate(job.header)
        if occursin(first(old_new), line)
            job.header[i] = replace(line, old_new)
            s = """Old line:
            $line
            New line:
            $(job.header[i])
            """
            print && (@info s)
        end
    end
    return job
end

#TODO
"""
    set_wanenergies!(job::Job, nscf::Calculation{QE}, projwfc::Calculation{QE}, threshold::Real; Epad=5.0)

Will set the energy window limits correctly according to the projections specified in the
structure of the job. The output of `projwfc` and the `threshold` will be used to determine
the minimum limit of the frozen energy window such that the interesting DOS of inside it exceeds
the threshold. `nscf` will be used to determine the DOS, and what the upper limit of the frozen window
needs to be to fit enough bands inside it, depending on the projections.
"""
function set_wanenergies!(job::Job, nscf::Calculation, projwfc::Calculation,
                          threshold::Real; Epad = 5.0)
    hasoutput_assert(projwfc)
    @assert Calculations.isprojwfc(projwfc) "Please specify a valid projwfc calculation."
    @assert Calculations.isnscf(nscf) "Please specify a valid nscf calculation."
    Emin = Emin_from_projwfc(job.structure, projwfc, threshold)
    return set_wanenergies!(job, nscf, Emin; Epad = Epad)
end

"""
    set_wanenergies!(job::Job, min_window_determinator::Real; kwargs...)

Sets the energy windows of wannier calculations based on the `job`.
When a projwfc calculation is present in the `job`, `min_window_determinator` will be used to
determine the threshold value for including a band in the window based on the projections, otherwise
it will be used as the `Emin` value from which to start counting the number of bands needed for all
projections.
"""
function set_wanenergies!(job::Job, min_window_determinator::Real; kwargs...)
    nscf_calc = getfirst(Calculations.isnscf, job.calculations)
    projwfc_calc = getfirst(Calculations.isprojwfc, job.calculations)
    if projwfc_calc === nothing || !hasoutput(projwfc_calculation)
        @info "No projwfc calculation found with valid output, using $min_window_determinator as Emin"
        return set_wanenergies!(job, nscf_calc, min_window_determinator; kwargs...)
    else
        @info "Valid projwfc output found, using $min_window_determinator as the dos threshold."
        return set_wanenergies!(job, nscf_calc, projwfc_calculation,
                                min_window_determinator; kwargs...)
    end
end

"""
    bandgap(job::Job, fermi=nothing)

Calculates the bandgap (possibly indirect) around the fermi level.
Uses the first found bands calculation, if there is none it uses the first found nscf calculation.
"""
function bandgap(job::Job, fermi = nothing)
    band_calcs = filter(!isnothing, getfirst.([Calculations.isbands, Calculations.isnscf, Calculations.isscf], (job.calculations,)))
    if isempty(band_calcs)
        error("No valid calculation found to calculate the bandgap.\nMake sure the job has either a valid bands or nscf calculation.")
    end

    outdat = outputdata(job)
    if fermi === nothing
        fermi_calcs = filter(!isnothing, getfirst.([Calculations.isnscf, Calculations.isscf], (job.calculations,)))
        if isempty(fermi_calcs)
            error("No valid calculation found to extract the fermi level.\nPlease supply the fermi level manually.")
        end
        fermi = maximum(x->get(get(outdat, x.name, Dict()), :fermi, -Inf), fermi_calcs)
    end

    bands = map(x->outdata[x.name][:bands], filter(x -> haskey(outdat, x.name), band_calcs))
    return minimum(bandgap.(bands, fermi))
end

function readfermi(job::Job, outdat)
    ins = filter(x -> (Calculations.isvcrelax(x) || Calculations.isscf(x) || Calculations.isnscf(x)), job.calculations)
    @assert isempty(ins) !== nothing "Job does not have a valid scf or nscf output."
    for i in ins
        if haskey(outdat, i.name)
            o = outdat[i.name]
            if haskey(o, :fermi)
                return o[:fermi]
            end
        end
    end
    @warn "No output files with fermi level found."
    return 0.0
end

function readbands(job::Job, outdat)
    calc = getfirst(x -> Calculations.isbands(x), job.calculations)
    outdat = outputdata(job)
    if calc === nothing || !haskey(outdat, calc.name)
        calc = getfirst(x -> Calculations.isnscf(x), job.calculations)
        if calc === nothing || !haskey(outdat, calc.name)
            @warn "Job does not have a valid bands output."
            return nothing
        end
        @warn "No bands calculation found, return bands from nscf calculation."
        return outdat[calc.name][:bands]
    end
    return outdat[calc.name][:bands]
end

#TODO: only for QE
"""
    set_present!(job::Job, func::Function)
    set_present!(job::Job, func::String)
    set_present!(job::Job, func::Expr)

Sets a function with the call signature `func(job)` which can be later called using the [`@present`](@ref) macro.
"""
function set_present!(job::Job, func::Function)
    try
        str = loaded_modules_string() * @code_string func(job)
        set_present!(job, str)
    catch
        error("Could not generate the source string for the supplied function.\nIf you are running in a Jupyter notebook, please supply the source code as a string, or expression to set_present!.\nDon't forget to include the right using statements.")
    end
end
function set_present!(job::Job, func::AbstractString)
    open(joinpath(job, ".present.jl"), "w") do f
        return write(f, func)
    end
end
function set_present!(job::Job, func::Expr)
    funcstr = string(func)
    if funcstr[1:5] == "begin"
        funcstr = funcstr[findfirst(isequal('\n'), funcstr)+1:findlast(isequal('\n'), funcstr)-1]
    end
    return set_present!(job, funcstr)
end

"""
    present(job)

Calls a present function if it was previously saved using [`set_present!`](@ref) or [`archive`](@ref).
"""
macro present(job)
    return esc(quote
                   if ispath(joinpath($job, ".present.jl"))
                       t = include(joinpath($job, ".present.jl"))
                       DFControl.with_logger(DFControl.MinLevelLogger(DFControl.current_logger(),
                                                                      DFControl.Logging.Error)) do
                           return t($job)
                       end
                   else
                       @error "No presentation function defined.\n Please set it with `set_present!`."
                   end
               end)
end

"""
    archive(job::Job, archive_directory::AbstractString, description::String=""; present = nothing, version=job.version)

Archives `job` by copying it's contents to `archive_directory` alongside a `results.jld2` file with all the parseable results as a Dict. `description` will be saved in a `description.txt` file in the `archive_directory`. A different job version can be copied using the `version` kwarg, and with the `present` kwarg a function can be specified that can be later called with the [`@present`](@ref) macro.
"""
function archive(job::Job, archive_directory::AbstractString, description::String = "";
                 present = nothing, version = job.version)
    @assert !isarchived(job) "Job was already archived"
    final_dir = config_path(".archived", archive_directory)
    @assert !ispath(final_dir) "A archived job already exists in $archive_directory"
    mkpath(final_dir)

    present !== nothing && set_present!(job, present)
    out = outputdata(job)
    tj = deepcopy(job)
    switch_version!(tj, version)
    cp(tj, final_dir)
    tj.dir = final_dir

    JLD2.jldsave(joinpath(final_dir, "results.jld2"); outputdata=out)

    !isempty(description) && write(joinpath(final_dir, "description.txt"), description)
    push!(JOB_REGISTRY.archived, tj.dir)
    @info "Archived job at $(tj.dir). If you're done with this one, it is safe to delete the directory at $(job.dir)."
    write_job_registry()
    return nothing
end

for (f, strs) in zip((:cp, :mv), (("copy", "Copies"), ("move", "Moves")))
    @eval begin
        """
            $($f)(job::Job, dest::AbstractString; all=false, temp=false, kwargs...)

        $($(strs[2])) the contents of `job.dir` to `dest`. If `all=true`, it will also $($(strs[1])) the
        `.version` directory with all previous versions. If `temp=true` it will override
        `job.copy_temp_folders` and $($(strs[1])) also the temporary calculation directories.
        The `kwargs...` are passed to `Base.$($f)`.
        """
        function Base.$f(job::Job, dest::AbstractString; all = false, temp = false,
                         kwargs...)
            if !ispath(dest)
                mkpath(dest)
            end
            for file in readdir(job)
                if !all
                    if file == VERSION_DIR_NAME
                        continue
                    elseif file == TEMP_CALC_DIR && !(temp || job.copy_temp_folders)
                        continue
                    end
                end
                if joinpath(job, file) == abspath(dest)
                    continue
                end
                $f(joinpath(job, file), joinpath(dest, file); kwargs...)
            end
        end
    end
end

function pseudo_cutoffs(str::String)
    lines = split(str, "\n")
    id = findfirst(x->occursin("Suggested minimum cutoff for wavefunctions:", x), lines)
    if id !== nothing
        return parse(Float64, split(lines[id])[end-1]), parse(Float64, split(lines[id+1])[end-1])
    else
        return 0.0, 0.0
    end
end

function sanitize_cutoffs!(job::Job)
    ψcut, ρcut = 0.0, 0.0
    # the assumption is that the most important cutoff calculation is the scf/vcrelax that is ran first
    ψ_cut_calc = getfirst(x -> haskey(x, Calculations.ψ_cutoff_flag(x)), job.calculations)
    if ψ_cut_calc !== nothing
        ψcut = ψ_cut_calc[Calculations.ψ_cutoff_flag(ψ_cut_calc)]
    else
        pseudos = unique(map(x -> x.pseudo, job.structure.atoms))
        all_cuts = map(x -> pseudo_cutoffs(x), pseudos)
        ψcut = maximum(x -> x[1], all_cuts)
        ρcut = maximum(x -> x[2], all_cuts)
        @assert ψcut != 0.0 "No energy cutoff was specified in any calculation, and the calculated cutoff from the pseudopotentials was 0.0.\nPlease manually set one."
        @info "No energy cutoff was specified in the scf calculation.\nCalculated ψcut=$ψcut."
    end
    ρ_cut_calc = getfirst(x -> Calculations.hasflag(x, Calculations.ρ_cutoff_flag(x)),
                          job.calculations)
    if ρ_cut_calc !== nothing
        ρcut = ρ_cut_calc[Calculations.ρ_cutoff_flag(ρ_cut_calc)]
    end
    for i in job.calculations
        ψflag = Calculations.ψ_cutoff_flag(i)
        ψflag !== nothing &&
            !haskey(i, ψflag) &&
            Calculations.set_flags!(i, ψflag => ψcut; print = false)
        ρflag = Calculations.ρ_cutoff_flag(i)
        ρflag !== nothing && !haskey(i, ρflag) && ρcut != 0.0 &&
            Calculations.set_flags!(i, ρflag => ρcut; print = false)
    end
end


text:
#
# Correctness Tests
#

module TestCategoricalArrays
    fatalerrors = length(ARGS) > 0 && ARGS[1] == "-f"
    quiet = length(ARGS) > 0 && ARGS[1] == "-q"
    anyerrors = false

    using Base.Test
    using CategoricalArrays

    tests = [
        "01_typedef.jl",
        "02_buildorder.jl",
        "03_buildfields.jl",
        "03_constructors.jl",
        "05_convert.jl",
        "06_show.jl",
        "06_length.jl",
        "07_levels.jl",
        "08_equality.jl",
        "09_hash.jl",
        "10_isless.jl",
        "11_array.jl",
        "12_nullablearray.jl",
        "13_arraycommon.jl"
    ]

    println("Running tests:")

    for test in tests
        try
            include(test)
            println("\t\033[1m\033[32mPASSED\033[0m: $(test)")
        catch e
            anyerrors = true
            println("\t\033[1m\033[31mFAILED\033[0m: $(test)")
            if fatalerrors
                rethrow(e)
            elseif !quiet
                showerror(STDOUT, e, backtrace())
                println()
            end
        end
    end

    if anyerrors
        throw("Tests failed")
    end
end


text:
# Support for physical units: I/O and conversion

module Units

import Base: *, /, convert, show

export  SIPrefix,
        Quantity,
        Unit,
# Prefixes
        Yocto,
        Zepto,
        Atto,
        Femto,
        Pico,
        Nano,
        Micro,
        Milli,
        Centi,
        Deci,
        SINone,
        Deca,
        Hecto,
        Kilo,
        Mega,
        Giga,
        Tera,
        Peta,
        Exa,
        Zetta,
        Yotta,
# Units
        UnitBase,
        Unknown,
# Length units
        Meter,
        Angstrom,
        Inch,
        Foot,
        Yard,
        Mile,
        LightYear,
        Parsec,
        AstronomicalUnit,
# Time units
        Second,
        Minute,
        Hour,
        Day,
        Week,
        YearJulian,
        PlanckTime,
# Rate units
        Hertz,
# Mass units
        Gram,
        AtomicMassUnit,
        Dalton,
        PlanckMass,
# Electric current
        Ampere,
# Temperature
        Kelvin,
        Celsius,
        Fahrenheit,
# Luminosity
        Candela,
# Amount
        Mole,
        Entities,
#
# Functions
        base,
        fshow,
        lshow,
        prefix,
        pshow,
        parse_quantity,
        sshow

abstract SIPrefix
immutable Yocto <: SIPrefix end
immutable Zepto <: SIPrefix end
immutable Atto <: SIPrefix end
immutable Femto <: SIPrefix end
immutable Pico <: SIPrefix end
immutable Nano <: SIPrefix end
immutable Micro <: SIPrefix end
immutable Milli <: SIPrefix end
immutable Centi <: SIPrefix end
immutable Deci <: SIPrefix end
immutable SINone <: SIPrefix end
immutable Deca <: SIPrefix end
immutable Hecto <: SIPrefix end
immutable Kilo <: SIPrefix end
immutable Mega <: SIPrefix end
immutable Giga <: SIPrefix end
immutable Tera <: SIPrefix end
immutable Peta <: SIPrefix end
immutable Exa <: SIPrefix end
immutable Zetta <: SIPrefix end
immutable Yotta <: SIPrefix end

# PrettyShow and PrettyString
pshow(x) = pshow(STDOUT::IO, x)
function pstring(x)
    s = IOBuffer()
    pshow(s, x)
    takebuf_string(s)
end
# LatexShow and LatexString
lshow(x) = lshow(STDOUT::IO, x)
function lstring(x)
    s = IOBuffer()
    lshow(s, x)
    takebuf_string(s)
end
# FullShow and FullString
fshow(x) = fshow(STDOUT::IO, x)
function fstring(x)
    s = IOBuffer()
    fshow(s, x)
    takebuf_string(s)
end

let
#  Prefix        ToSINone      SimpleShow    PrettyShow   LatexShow  Full
const prefix_table = {
  (Yocto,        1e-24,        "y",          "y",         "y",       "yocto")
  (Zepto,        1e-21,        "z",          "z",         "z",       "zepto")
  (Atto,         1e-18,        "a",          "a",         "a",       "atto")
  (Femto,        1e-15,        "f",          "f",         "f",       "femto")
  (Pico,         1e-12,        "p",          "p",         "p",       "pico")
  (Nano,         1e-9,         "n",          "n",         "n",       "nano")
  (Micro,        1e-6,         "u",          "\u03bc",    "\$\\mu\$","micro")
  (Milli,        1e-3,         "m",          "m",         "m",       "milli")
  (Centi,        1e-2,         "c",          "c",         "c",       "centi")
  (Deci,         1e-1,         "d",          "d",         "d",       "deci")
  (Deca,         1e1,          "da",         "da",        "da",      "deca")
  (Hecto,        1e2,          "h",          "h",         "h",       "hecto")
  (Kilo,         1e3,          "k",          "k",         "k",       "kilo")
  (Mega,         1e6,          "M",          "M",         "M",       "mega")
  (Giga,         1e9,          "G",          "G",         "G",       "giga")
  (Tera,         1e12,         "T",          "T",         "T",       "tera")
  (Peta,         1e15,         "P",          "P",         "P",       "peta")
  (Exa,          1e18,         "E",          "E",         "E",       "exa")
  (Zetta,        1e21,         "Z",          "Z",         "Z",       "zetta")
  (Yotta,        1e24,         "Y",          "Y",         "Y",       "yotta")
}
global to_reference
global pshow
global sshow
global lshow
global fshow
global _unit_si_prefixes
for (t, f, ss, ps, ls, fs) in prefix_table
    @eval to_reference(::Type{$t}) = $f
    @eval sshow(io::IO, ::Type{$t}) = print(io, $ss)
    @eval pshow(io::IO, ::Type{$t}) = print(io, $ps)
    @eval lshow(io::IO, ::Type{$t}) = print(io, $ls)
    @eval fshow(io::IO, ::Type{$t}) = print(io, $fs)
    # Default to prettyshow
    @eval show(io::IO, ::Type{$t}) = print(io, $ps)
end
function return_prefix(i::Int)
    if i <= length(prefix_table)
        return prefix_table[i][1]
    elseif i == length(prefix_table)+1
        return SINone
    else
        error("Something is wrong")
    end
end
_unit_si_prefixes = ntuple(length(prefix_table)+1, i->return_prefix(i))
end  # let
to_reference(::Type{SINone}) = 1
sshow(io::IO, ::Type{SINone}) = nothing
pshow(io::IO, ::Type{SINone}) = nothing
lshow(io::IO, ::Type{SINone}) = nothing
fshow(io::IO, ::Type{SINone}) = nothing

# Units
abstract UnitBase
immutable Unit{TP<:SIPrefix, TU<:UnitBase} end
Unit{TP<:SIPrefix, TU<:UnitBase}(tp::Type{TP}, tu::Type{TU}) = Unit{tp, tu}
Unit{TU<:UnitBase}(tu::Type{TU}) = Unit{SINone, tu}
function show{TP<:SIPrefix, TU<:UnitBase}(io::IO, tu::Unit{TP, TU})
    print(io, TP)
    print(io, TU)
end

###################################
###         Quantity            ###
###################################
# A quantity is a value with a unit; it preserves the memory
# of the unit specified by the user. However, quantities cannot
# be combined to make compound units such as Joule
immutable Quantity{TP<:SIPrefix, TU<:UnitBase, Tdata}
    value::Tdata
end
Quantity{TP<:SIPrefix, TU<:UnitBase, Tdata}(tp::Type{TP}, tu::Type{TU}, val::Tdata) = Quantity{tp, tu, Tdata}(val)
Quantity{TU<:UnitBase, Tdata}(tu::Type{TU}, val::Tdata) = Quantity{SINone, tu, Tdata}(val)
prefix{TP<:SIPrefix, TU<:UnitBase, Tdata}(q::Quantity{TP, TU, Tdata}) = TP
base{TP<:SIPrefix, TU<:UnitBase, Tdata}(q::Quantity{TP, TU, Tdata}) = TU

function sshow{TP<:SIPrefix, TU<:UnitBase}(io::IO, q::Quantity{TP, TU})
    print(io, q.value, " ")
    sshow(io, TP)
    sshow(io, TU)
end
function pshow{TP<:SIPrefix, TU<:UnitBase}(io::IO, q::Quantity{TP, TU})
    print(io, q.value, " ")
    pshow(io, TP)
    pshow(io, TU)
end
function lshow{TP<:SIPrefix, TU<:UnitBase}(io::IO, q::Quantity{TP, TU})
    print(io, q.value, " ")
    lshow(io, TP)
    lshow(io, TU)
end
function fshow{TP<:SIPrefix, TU<:UnitBase}(io::IO, q::Quantity{TP, TU})
    print(io, q.value, " ")
    fshow(io, TP)
    fshow(io, TU)
end
show(io::IO, q::Quantity) = pshow(io, q)

# Functions and dictionaries for parsing
_unit_string_dict = (String=>Tuple)[]

function _unit_gen_func_multiplicative(table)
    for (t, r, to_r, s, ps, ls, fs) in table
        @eval reference(::Type{$t}) = $r
        @eval to_reference(::Type{$t}) = x->x*$to_r
        @eval from_reference(::Type{$t}) = x->x/$to_r
        @eval show(io::IO, ::Type{$t}) = print(io, $s)
        @eval pshow(io::IO, ::Type{$t}) = print(io, $ps)
        @eval lshow(io::IO, ::Type{$t}) = print(io, $ls)
        @eval fshow(io::IO, ::Type{$t}) = print(io, $fs)
    end
end

function _unit_gen_func(table)
    for (t, r, to_func, from_func, s, ps, ls, fs) in table
        @eval reference(::Type{$t}) = $r
        @eval to_reference(::Type{$t}) = $to_func
        @eval from_reference(::Type{$t}) = $from_func
        @eval show(io::IO, ::Type{$t}) = print(io, $s)
        @eval pshow(io::IO, ::Type{$t}) = print(io, $ps)
        @eval lshow(io::IO, ::Type{$t}) = print(io, $ls)
        @eval fshow(io::IO, ::Type{$t}) = print(io, $fs)
    end
end

function _unit_gen_dict_with_prefix(table)
    for (u, rest) in table
        ustr = string(u)
        upstr = pstring(u)
        for p in _unit_si_prefixes
            pstr = string(p)
            key = pstr*ustr
#            println(key, ": ", p, ", ", u)
            _unit_string_dict[key] = (p, u)
            # Add "pretty" variants, too
            ppstr = pstring(p)
            _unit_string_dict[pstr*upstr] = (p, u)
            _unit_string_dict[ppstr*ustr] = (p, u)
            _unit_string_dict[ppstr*upstr] = (p, u)
        end
    end
end

function _unit_gen_dict(table)
    for (u, rest) in table
        key = string(u)
#        println(key, ": ", u)
        _unit_string_dict[key] = (SINone, u)
        upstr = pstring(u)
        if key != upstr
            _unit_string_dict[upstr] = (SINone, u)
        end
    end
end

# Unknown unit
immutable Unknown <: UnitBase end

# Length units
abstract UnitLength <: UnitBase
immutable Meter <: UnitLength end
immutable Angstrom <: UnitLength end
immutable Inch <: UnitLength end
immutable Foot <: UnitLength end
immutable Yard <: UnitLength end
immutable Mile <: UnitLength end
immutable LightYear <: UnitLength end
immutable Parsec <: UnitLength end
immutable AstronomicalUnit <: UnitLength end
let
  # Unit     RefUnit    ToRef            show      pshow     lshow   fshow
const utable = {
  (Meter,    Meter,    1,                "m",      "m",      "m",    "meter")
  (Angstrom, Meter,    1e-10,            "A",      "\u212b", "\$\\AA\$","angstrom")
  (Inch,     Meter,    25.4/1000,        "in",     "in",     "in",   "inch")
  (Foot,     Meter,    12*25.4/1000,     "ft",     "ft",     "ft",   "foot")
  (Yard,     Meter,    36*25.4/1000,     "yd",     "yd",     "yd",   "yard")
  (Mile,     Meter,    5280*12*25.4/1000,"mi",     "mi",     "mi",   "mile")
  (LightYear,Meter,    9.4605284e15,     "ly",     "ly",     "ly",   "lightyear")
  (Parsec,   Meter,    3.08568025e16,    "pc",     "pc",     "pc",   "parsec")
  (AstronomicalUnit, Meter, 149_597_870_700, "AU", "AU",     "AU",   "astronomical unit")
}
_unit_gen_func_multiplicative(utable)
_unit_gen_dict_with_prefix({utable[1]})  # parse prefixes only for Meter
_unit_gen_dict(utable[2:end])
end

# Time units
abstract UnitTime <: UnitBase
immutable Second <: UnitTime end
immutable Minute <: UnitTime end
immutable Hour <: UnitTime end
immutable Day <: UnitTime end
immutable Week <: UnitTime end
immutable YearJulian <: UnitTime end
immutable PlanckTime <: UnitTime end
let
  # Unit       RefUnit     ToRef            show      pshow     lshow  fshow
const utable = {
  (Second,     Second,      1,             "s",      "s",      "s",    "second")
  (Minute,     Second,      60,            "min",    "min",    "min",  "minute")
  (Hour,       Second,      3600,          "hr",     "hr",     "hr",   "hour")
  (Day,        Second,      86400,         "d",      "d",      "d",    "day")
  (YearJulian, Second,      365.25*86400,  "yr",     "yr",     "yr",   "year")
  (PlanckTime, Second,      5.3910632e-44, "tP",     "tP",     "\$t_P\$","Planck time")
}
_unit_gen_func_multiplicative(utable)
_unit_gen_dict_with_prefix({utable[1]})   # prefixes are only common for seconds
_unit_gen_dict(utable[2:end])
end

# Rate units
immutable Hertz <: UnitBase end
let
  # Unit       RefUnit     ToRef        show      pshow     lshow   fshow
const utable = {
  (Hertz,      Hertz,      1,           "Hz",     "Hz",     "Hz",   "Hertz")
}
_unit_gen_func_multiplicative(utable)
_unit_gen_dict_with_prefix(utable)
end

# Mass units
# (note English units like pounds are technically weight, not mass)
abstract UnitMass <: UnitBase
immutable Gram <: UnitMass end
immutable AtomicMassUnit <: UnitMass end
immutable Dalton <: UnitMass end
immutable PlanckMass <: UnitMass end
let
  # Unit       RefUnit     ToRef           show      pshow     lshow  fshow
const utable = {
  (Gram,       Gram,        1,             "g",      "g",      "g",   "gram")
  (AtomicMassUnit, Gram,    1.66053892173e-24, "amu","amu",    "amu", "atomic mass unit")
  (Dalton,     Gram,        1.66053892173e-24, "Da", "Da",     "Da",  "Dalton")
  (PlanckMass, Gram,        2.1765113e-5,  "mP",     "mP",     "\$m_P\$","Planck mass")
}
_unit_gen_func_multiplicative(utable)
_unit_gen_dict_with_prefix({utable[1]})
_unit_gen_dict(utable[2:end])
end

# Electric current units
immutable Ampere <: UnitBase end
let
  # Unit       RefUnit     ToRef        show      pshow     lshow   fshow
const utable = {
  (Ampere,     Ampere,     1,           "A",      "A",      "A",   "ampere")
}
_unit_gen_func_multiplicative(utable)
_unit_gen_dict_with_prefix(utable)
end

# Temperature units
# The conversions for these are not just a product
immutable Kelvin <: UnitBase end
immutable Celsius <: UnitBase end
immutable Fahrenheit <: UnitBase end
let
  # Unit       RefUnit  ToRef              FromRef     show    pshow  lshow  fshow
const utable = {
  (Kelvin,     Kelvin,  x->x,              x->x,       "K",    "K",   "K",   "Kelvin")
  (Celsius,    Kelvin,  x->x+273.15,       x->x-273.15,"C",    "C",   "C",   "Celsius")
  (Fahrenheit, Kelvin,  x->(x+459.67)*5/9, x->x*9/5-459.67,"F","F",   "F",   "Fahrenheit")
}
_unit_gen_func(utable)
_unit_gen_dict_with_prefix({utable[1]})   # prefixes are only common for Kelvin
_unit_gen_dict(utable[2:end])
end

# Luminosity units
immutable Candela <: UnitBase end
let
  # Unit       RefUnit     ToRef        show      pshow     lshow   fshow
const utable = {
  (Candela,    Candela,    1,           "cd",     "cd",     "cd",   "candela")
}
_unit_gen_func_multiplicative(utable)
_unit_gen_dict_with_prefix(utable)
end

# Amount units
immutable Mole <: UnitBase end
immutable Entities <: UnitBase end
let
  # Unit       RefUnit     ToRef        show      pshow     lshow   fshow
const utable = {
  (Mole,       Mole,       1,           "mol",    "mol",    "mol",  "mole")
  (Entities,   Mole,       1/6.0221417930e23, "", "",       "",     "entities")
}
_unit_gen_func_multiplicative(utable)
_unit_gen_dict_with_prefix({utable[1]})
end

# Parsing string to extract Quantity
function parse_quantity(s::String, strict::Bool)
    # Find the last character of the numeric component
    m = match(r"[0-9\.\+-](?![0-9\.\+-])", s)
    if m == nothing
        error("String does not have a 'value unit' structure")
    end
    val = float64(s[1:m.offset])
    ustr = strip(s[m.offset+1:end])
    if isempty(ustr)
        if strict
            error("String does not have a 'value unit' structure")
        else
            return Quantity(SINone, Unknown, val)
        end
    end
    (prefix, unit) = _unit_string_dict[ustr]
    return Quantity(prefix, unit, val)
end
parse_quantity(s::String) = parse_quantity(s, true)

# Conversion between units
function convert{Uin<:UnitBase, Uout<:UnitBase, Pin<:SIPrefix, Pout<:SIPrefix, T}(::Type{Unit{Pout, Uout}}, q::Quantity{Pin, Uin, T})
    if reference(Uin) != reference(Uout)
        error("Not convertable")
    end
    return Quantity(Pout, Uout, from_reference(Uout)(to_reference(Uin)(q.value*to_reference(Pin))) / to_reference(Pout))
end

# Arithmetic with Quantities
(*){TP,TU}(n::Number, q::Quantity{TP,TU}) = Quantity(TP, TU, n*q.value)
(*){TP,TU}(q::Quantity{TP,TU}, n::Number) = Quantity(TP, TU, n*q.value)
(/){TP,TU}(q::Quantity{TP,TU}, n::Number) = Quantity(TP, TU, q.value/n)

(/){TP,TU<:UnitLength}(q1::Quantity{TP,TU}, q2::Quantity{TP,TU}) = q1.value/q2.value
(/){TP1,TP2,TU<:UnitLength}(q1::Quantity{TP1,TU}, q2::Quantity{TP2,TU}) = q1.value/q2.value*to_reference(TP1)/to_reference(TP2)
(/){TP1,TP2,TU1<:UnitLength,TU2<:UnitLength}(q1::Quantity{TP1,TU1}, q2::Quantity{TP2,TU2}) = from_reference(TU2)(to_reference(TU1)(q1.value/q2.value*to_reference(TP1)/to_reference(TP2)))
(/){TP,TU<:UnitTime}(q1::Quantity{TP,TU}, q2::Quantity{TP,TU}) = q1.value/q2.value
(/){TP1,TP2,TU<:UnitTime}(q1::Quantity{TP1,TU}, q2::Quantity{TP2,TU}) = q1.value/q2.value*to_reference(TP1)/to_reference(TP2)
(/){TP1,TP2,TU1<:UnitTime,TU2<:UnitTime}(q1::Quantity{TP1,TU1}, q2::Quantity{TP2,TU2}) = from_reference(TU2)(to_reference(TU1)(q1.value/q2.value*to_reference(TP1)/to_reference(TP2)))
(/){TP,TU<:UnitMass}(q1::Quantity{TP,TU}, q2::Quantity{TP,TU}) = q1.value/q2.value
(/){TP1,TP2,TU<:UnitMass}(q1::Quantity{TP1,TU}, q2::Quantity{TP2,TU}) = q1.value/q2.value*to_reference(TP1)/to_reference(TP2)
(/){TP1,TP2,TU1<:UnitMass,TU2<:UnitMass}(q1::Quantity{TP1,TU1}, q2::Quantity{TP2,TU2}) = from_reference(TU2)(to_reference(TU1)(q1.value/q2.value*to_reference(TP1)/to_reference(TP2)))

(/){TP,TU}(q1::Quantity{TP,TU}, q2::Quantity{TP,TU}) = q1.value/q2.value
(/){TP1,TP2,TU}(q1::Quantity{TP1,TU}, q2::Quantity{TP2,TU}) = q1.value/q2.value*to_reference(TP1)/to_reference(TP2)

end


text:
function do_unsupported()
    a = 1
    # dict comprehension ((key + 1), (value + 1)) unimplemented on line 9:4;
    b = Bool(a)
    println(join([b ? ("True") : ("False")], " "))
end

function main()
    do_unsupported()
end

main()


text:
@testset "fit_adc.jl" begin

end


text:
# ==============================================================================
# Modele JuMP pour calculer la relaxation linéaire du 2SPA sur un objectif donne, avec eventuellement une ϵ-contrainte
function computeLinearRelax2SPA(nbvar::Int, nbcontraintes::Int, L::Array{Int,2}, c1::Array{Int,1}, c2::Array{Int,1}, epsilon, obj::Int)

    model = Model(with_optimizer(GLPK.Optimizer))
#    model = Model(with_optimizer(Gurobi.Optimizer, GUROBI_ENV))
#    set_optimizer_attributes(model, "OutputFlag" => 0)

    @variable(model, 0.0 <= x[1:nbvar] <= 1.0 )
    @constraint(model, [i=1:nbcontraintes],(sum((x[j]*L[i,j]) for j in 1:nbvar)) == 1)
    if obj == 1
        @objective(model, Min, sum((c1[i])*x[i] for i in 1:nbvar))
        @constraint(model, sum((c2[i])*x[i] for i in 1:nbvar) <= epsilon)
    else
        @objective(model, Min, sum((c2[i])*x[i] for i in 1:nbvar))
        @constraint(model, sum((c1[i])*x[i] for i in 1:nbvar) <= epsilon)
    end
    optimize!(model)
    return objective_value(model), value.(x)
end

max_stars_repo_name: ducanh-le/GravityMachine

max_issues_repo_name: ducanh-le/GravityMachine

max_forks_repo_name: ducanh-le/GravityMachine


text:
@testset "Elastic Cholesky" begin
    T = Float64
    max_size = 4
    A = rand(MersenneTwister(0), max_size, max_size)
    A = A'A

    # Test incrementally adding columns of `A` to `ElasticCholesky`
    for colperms in [[1,2,3,4], [4,3,2,1], [1,3,2,4], [4,1,3,2]]
        C_copy_into_A = ElasticCholesky{T}(max_size)
        C_wrap_A = ElasticCholesky(A)
        for j in 1:max_size
            # Add (permuted) column `colperms[j]` of `A` to the `colperms[j]`th column of `C.A`
            cholesky!(C_copy_into_A, colperms[j], A[colperms[1:j], colperms[j]], Val(true))
            cholesky!(C_wrap_A, colperms[j])

            for C in [C_copy_into_A, C_wrap_A]
                # Check the fields of `C`
                J = colperms[1:C.ncols[]]
                @test Hermitian(C.A[J, J], :U) ≈ A[J, J]
                @test UpperTriangular(C.U[J, J]) ≈ cholesky(A[J, J]).U
                @test C.colperms[1:C.ncols[]] == J
                @test C.ncols[] == j

                # Check `ldiv!` gives same result as `LinearAlgebra.cholesky`
                b = rand(j)
                @test ldiv!(similar(b), C, b) ≈ cholesky(A[J, J]) \ b
            end
        end
    end

    # Test incremental factorization of wrapped array
    C = ElasticCholesky(copy(A))
    cholesky!(C)

    @test Hermitian(C.A, :U) ≈ A
    @test UpperTriangular(C.U) ≈ cholesky(A).U
    @test C.colperms == 1:max_size
    @test C.ncols[] == max_size

    b = rand(max_size)
    @test ldiv!(similar(b), C, b) ≈ cholesky(A) \ b
end

function random_nodes(n, T, max_size)
    min_bound   = -rand(SVector{n,T})
    max_bound   = rand(SVector{n,T})
    rand_node() = min_bound .+ rand(SVector{n,T}) .* (max_bound .- min_bound)
    rand_dir()  = (x = rand_node(); x/norm(x))

    nodes     = [min_bound, max_bound, (rand_node() for i in 3:max_size)...]
    values    = rand(T, max_size)
    d_nodes   = [rand_node() for _ in 1:n*max_size]
    d_dirs    = [rand_dir() for _ in 1:n*max_size]
    d_values  = rand(T, n*max_size)
    return (; min_bound, max_bound, nodes, values, d_nodes, d_dirs, d_values)
end

@testset "Elastic Gram Matrix" begin
    nhs = NormalHermiteSplines
    max_size = 4
    T = Float64

    @testset "RK_H0" begin
        kernel = RK_H0(0.5 + rand())
        for n in 1:3
            # Test incremental building of Gram matrix, inserting nodes in random order
            @unpack nodes = random_nodes(n, T, max_size)
            n₁ = length(nodes)
            A  = nhs._gram(nodes, kernel)
            A′ = zeros(T, n₁, n₁)
            J  = Int[]
            for j in randperm(n₁)
                push!(J, j)
                nhs._gram!(view(A′, J, J), nodes[J[end]], nodes[J[1:end-1]], kernel)
                @test Hermitian(A′[J,J], :U) ≈ A[J,J] # `A′[J,J]` is upper triangular, `A` is Hermitian
            end
        end
    end

    @testset "RK_H1" begin
        kernel = RK_H1(0.5 + rand())
        for n in 1:3
            # Test incremental building of Gram matrix, inserting (derivative-)nodes in random order
            @unpack nodes, d_nodes, d_dirs = random_nodes(n, T, max_size)
            n₁, n₂ = length(nodes), length(d_nodes)
            A  = nhs._gram(nodes, d_nodes, d_dirs, kernel)
            A′ = zeros(T, n₁+n₂, n₁+n₂)
            J, J₁, J₂ = Int[], Int[], Int[]
            for j in randperm(n₁+n₂)
                if j <= n₁
                    push!(J₁, j)
                    insert!(J, length(J₁), j)
                    nhs._gram!(view(A′, J, J), nodes[J₁[end]], nodes[J₁[1:end-1]], d_nodes[J₂], d_dirs[J₂], kernel)
                else
                    push!(J₂, j-n₁)
                    push!(J, j)
                    nhs._gram!(view(A′, J, J), d_nodes[J₂[end]], d_dirs[J₂[end]], nodes[J₁], d_nodes[J₂[1:end-1]], d_dirs[J₂[1:end-1]], kernel)
                end
                @test Hermitian(A′[J,J], :U) ≈ A[J,J] # `A′[J,J]` is upper triangular, `A` is Hermitian
            end
        end
    end
end

@testset "Elastic Normal Spline" begin
    nhs = NormalHermiteSplines
    max_size = 3
    T = Float64

    for n in 1:3
        @unpack min_bound, max_bound, nodes, values, d_nodes, d_dirs, d_values = random_nodes(n, T, max_size)
        rk_H0 = RK_H0(0.5 + rand())
        rk_H1 = RK_H1(0.5 + rand())
        espl_H0 = ElasticNormalSpline(min_bound, max_bound, max_size, rk_H0)
        espl_H1_0 = ElasticNormalSpline(min_bound, max_bound, max_size, rk_H1)
        espl_H1_1 = ElasticNormalSpline(min_bound, max_bound, max_size, rk_H1)

        for i in 1:max_size
            # Update `ElasticNormalSpline`
            nodes′, values′ = nodes[1:i], values[1:i]
            d_nodes′, d_dirs′, d_values′ = d_nodes[1:n*i], d_dirs[1:n*i], d_values[1:n*i]

            # Insert regular node
            insert!(espl_H0, nodes′[i], values′[i])
            insert!(espl_H1_0, nodes′[i], values′[i])

            # Insert `n` derivative nodes
            insert!(espl_H1_1, nodes′[i], values′[i])
            for k in n*(i-1).+(1:n)
                insert!(espl_H1_1, d_nodes′[k], d_dirs′[k], d_values′[k])
            end

            # Compute `NormalSpline`
            i == 1 && continue # `NormalSpline` requires at least two nodes′
            spl_H0  = interpolate(nodes′, values′, rk_H0)
            spl_H1_0 = interpolate(nodes′, values′, rk_H1)
            spl_H1_1 = interpolate(nodes′, values′, d_nodes′, d_dirs′, d_values′, rk_H1)

            for (espl, spl) in [
                    (espl_H0, spl_H0)
                    (espl_H1_0, spl_H1_0)
                    (espl_H1_1, spl_H1_1)
                ]
                C = espl._chol
                n₁ = espl._num_nodes[]
                J = C.colperms[1:C.ncols[]]

                # Cholesky factorization of `ElasticNormalSpline` is built incrementally in arbitrary column order of the underlying Gram matrix;
                # compare with the Cholesky factorization of the gram matrix from `NormalSpline`, permuted appropriately
                J′ = (j -> ifelse(j > max_size,  j - max_size + n₁, j)).(J)
                C′ = cholesky(nhs._get_gram(spl)[J′, J′])
                b = randn(C.ncols[])

                @test nhs._get_kernel(espl)      == nhs._get_kernel(spl)
                @test nhs._get_nodes(espl)       ≈ nhs._get_nodes(spl)
                @test nhs._get_values(espl)      ≈ nhs._get_values(spl)
                @test nhs._get_d_nodes(espl)     ≈ nhs._get_d_nodes(spl)
                @test nhs._get_d_dirs(espl)      ≈ nhs._get_d_dirs(spl)
                @test nhs._get_d_values(espl)    ≈ nhs._get_d_values(spl)
                @test nhs._get_mu(espl)          ≈ nhs._get_mu(spl)
                @test nhs._get_gram(espl)        ≈ nhs._get_gram(spl)
                @test UpperTriangular(C.U[J, J]) ≈ C′.U
                @test ldiv!(similar(b), C, b)    ≈ C′ \ b
                # @test nhs._get_cond(espl)      ≈ nhs._get_cond(spl)
                @test nhs._get_min_bound(espl)   ≈ nhs._get_min_bound(spl)
                @test nhs._get_max_bound(espl)   ≈ nhs._get_max_bound(spl)
                @test nhs._get_scale(espl)       ≈ nhs._get_scale(spl)
            end
        end
    end
end


text:
include("learning_rate/alrsvi.jl")
include("learning_rate/inversedecay.jl")
include("analytic.jl")
include("analyticVI.jl")
include("gibbssampling.jl")
include("numericalVI.jl")

function post_process!(model::AbstractGP{<:Likelihood,<:Inference})
    nothing
end


text:
export RedPitayaClusterView, master, readDataPeriods, numChan, readDataSlow, readFrames, readPeriods, readPipelinedSamples, startPipelinedData, collectSamples!, readData, readDataPeriods, modeDAC, translateViewToClusterChannel

import Base: length, iterate, getindex, firstindex, lastindex

struct RedPitayaClusterView
    rpc::RedPitayaCluster
    view::Vector{Integer}
end

function RedPitayaClusterView(rpc::RedPitayaCluster, selection::Vector{Bool})
    length(selection) == length(rpc) || throw(DimensionMismatch("Length of boolean vector must match length of cluster"))
    view = findall(selection)
    return RedPitayaClusterView(rpc, view)
end

function RedPitayaClusterView(hosts::Vector{String}, selection, port=5025)
    rpc = RedPitayaCluster(hosts, port)
    return RedPitayaClusterView(rpc, selection)
end

length(rpcv::RedPitayaClusterView) = length(rpcv.view)
numChan(rpcv::RedPitayaClusterView) = 2 * length(rpcv.view)
master(rpcv::RedPitayaClusterView) = master(rpcv.rpc)

# Index
function getindex(rpcv::RedPitayaClusterView, index::Integer)
    1 <= index <= length(rpcv) || throw(BoundsError(rpcv.view, index))
    viewIndex = rpcv.view[index]
    1 <= viewIndex <= length(rpcv.rpc) || throw(BoundsError(rpcv.rpc, viewIndex))
    return rpcv.rpc[viewIndex]
end
firstindex(rpcv::RedPitayaClusterView) = start_(rpcv)
lastindex(rpcv::RedPitayaClusterView) = length(rpcv.view)

# Iterate
start_(rpcv::RedPitayaClusterView) = 1
next_(rpcv::RedPitayaClusterView, state) = (rpcv[state], state + 1)
done_(rpcv::RedPitayaClusterView, state) = state > length(rpcv)
iterate(rpcv::RedPitayaClusterView, s=start_(rpcv)) = done_(rpcv, s) ? nothing : next_(rpcv, s)

# Pass function calls on to RedPitayaCluster
for op in [:currentFrame, :currentPeriod, :currentWP, :connectADC, :startADC, :stopADC, :bufferSize]
    @eval $op(rpcv::RedPitayaClusterView) = $op(rpcv.rpc)
end

for op in [:periodsPerFrame, :samplesPerPeriod, :decimation, :keepAliveReset,
    :triggerMode, :slowDACStepsPerSequence, :samplesPerSlowDACStep,
    :slowDACStepsPerFrame, :ramWriterMode, :numSlowADCChan, :numSlowDACChan,
    :passPDMToFastDAC, :modeDAC, :masterTrigger]
    @eval $op(rpcv::RedPitayaClusterView) = $op(master(rpcv))
    @eval $op(rpcv::RedPitayaClusterView, value) = $op(rpcv.rpc, value)
end

for op in [:amplitudeDAC, :amplitudeDACNext, :frequencyDAC, :phaseDAC, :modulusFactorDAC, :modulusDAC]
    @eval $op(rpcv::RedPitayaClusterView, chan::Integer, component::Integer) = $op(rpcv.rpc, chan, component)
    @eval $op(rpcv::RedPitayaClusterView, chan::Integer, component::Integer, value) = $op(rpcv.rpc, chan, component, value)
end

for op in [:signalTypeDAC,  :DCSignDAC, :jumpSharpnessDAC, :setSlowDAC, :getSlowDAC, :offsetDAC]
    @eval $op(rpcv::RedPitayaClusterView, chan::Integer) = $op(rpcv.rpc, chan)
    @eval $op(rpcv::RedPitayaClusterView, chan::Integer, value) = $op(rpcv.rpc, chan, value)
end

function enableSlowDAC(rpcv::RedPitayaClusterView, enable::Bool, numFrames::Int64=0,
    ffRampUpTime::Float64=0.4, ffRampUpFraction::Float64=0.8)
    enableSlowDAC(rpcv.rpc, enable, numFrames, ffRampUpTime, ffRampUpFraction)
end

slowDACInterpolation(rpcv::RedPitayaClusterView, enable::Bool) = slowDACInterpolation(rpcv.rpc, enable)

# Returns the channel number in the cluster for a given channel number in the view
function translateViewToClusterChannel(rpcv::RedPitayaClusterView, chan)
    view = rpcv.view
    idxInView = div(chan -1, 2) + 1
    chanRP = mod1(chan, 2)
    return 2*(view[idxInView] - 1) + chanRP
end

# Helper functions for readData functions in Cluster.jl
function startPipelinedData(rpcv::RedPitayaClusterView, reqWP::Int64, numSamples::Int64, chunkSize::Int64)
    @sync for rp in rpcv
        @async startPipelinedData(rp, reqWP, numSamples, chunkSize)
    end
end

max_stars_repo_name: tknopp/RedPitayaDAQServer

max_issues_repo_name: tknopp/RedPitayaDAQServer

max_forks_repo_name: tknopp/RedPitayaDAQServer


text:
"""
COSINEKNN
---------
    compute the k-nearest neighbors similarity metric between the
    vertices of A or the upper half of a bipartite graph A.

Functions
---------
- S = cosineknn{T}(A::SparseMatrixCSC{T,Int64},K::Int64)
- S = cosineknn(A::MatrixNetwork,K::Int64)

Example
-------
~~~
A = load_matrix_network("bfs_example")
S = cosineknn(A,2)
~~~
"""

## CSC support
function cosineknn{T}(A::SparseMatrixCSC{T,Int64},K::Int64)
    (rp,ci,ai) = sparse_to_csr(A)
    (rpt,cit,ait) = sparse_to_csr(A')
    (m,n) = size(A)
    return cosineknn_internal(rp,ci,ai,rpt,cit,ait,m,K)
end

##  MatrixNetwork support
function cosineknn(A::MatrixNetwork,K::Int64)
    # for the original matrix
    (rp,ci,ai) = (A.rp,A.ci,A.vals)
    # for the transposed matrix:
    M = sparse_transpose(A)
    At = MatrixNetwork(M)
    (rpt,cit,ait) = (At.rp,At.ci,At.vals)
    m = A.n
    return cosineknn_internal(rp,ci,ai,rpt,cit,ait,m,K)

end

function cosineknn_internal{T}(rp::Vector{Int64},ci::Vector{Int64},ai::Vector{T},
                rpt::Vector{Int64},cit::Vector{Int64},ait::Vector{T},m::Int64,K::Int64)

    # accumarray
    rn = zeros(Float64,maximum(cit))
    for i = 1:length(ait)
        rn[cit[i]] += ait[i]^2
    end
    rn = sqrt(rn)

    rn[rn.>eps()] = 1./rn[rn.>eps()] # do division once

    si = zeros(Int64,m*K)
    sj = zeros(Int64,m*K)
    sv = zeros(Float64,m*K)
    nused = 0

    dwork = zeros(Float64,m)
    iwork = zeros(Int64,m)
    iused = zeros(Int64,m)

    for i = 1:m

        # for each row of A, compute an inner product against all rows of A.
        # to do this efficiently, we know that the inner-product will
        # only be non-zero if two rows of A overlap in a column.  Thus, in
        # row i of A, we look at all non-zero columns, and then look at all
        # rows touched by those columns in the A' matrix.  We  compute
        # the similarity metric, then sort and only store the top k.

        curused = 0 # track how many non-zeros we compute during this operation
        for ri = rp[i]:rp[i+1]-1
            # find all columns j in row i
            j = ci[ri]
            aij = ai[ri]*rn[i]
            for rit = rpt[j]:rpt[j+1]-1
                # find all rows k in column j
                k = cit[rit]
                if k == i
                    continue
                end # skip over the standard entries
                akj = ait[rit]*rn[k]
                if iwork[k]>0
                    # we already have a non-zero between row i and row k
                    dwork[iwork[k]] += aij*akj
                else
                    # we need to add a non-zero betwen row i and row k
                    curused = curused + 1
                    iwork[k] = curused
                    dwork[curused] = aij*akj
                    iused[curused] = k
                end
            end
        end

        # don't sort if fewer than K elements
        if curused < K
            sperm = 1:curused
        else
            sperm = sortperm(dwork[1:curused],rev=true)
        end

        for k = 1:min(K,length(sperm))
            nused = nused + 1
            si[nused] = i
            sj[nused] = iused[sperm[k]]
            sv[nused] = dwork[sperm[k]]
        end

        # reset the iwork array, no need to reset dwork, as it's just
        # overwritten
        for k = 1:curused
            iwork[iused[k]] = 0
        end
    end

    ei = si[1:nused]
    ej = sj[1:nused]
    ev = sv[1:nused]

    S = sparse(ei,ej,ev,m,m)

    return S
end


text:
choose_mesh = 5
#=
* * *
* * * ->    1
* * *

* * * * *
* * * * * ->  2
* * * * *

  ^
  |
* * *
* * *
* * *     3
* * *
* * *

* * *
* * *
* * *->     4
* * *
* * *
=#
if choose_mesh ==1
nodes =
[
0 0 0
1 0 0
1 1 0
0 1 0
.5 0 0
1 .5 0
.5 1 0
0 .5 0
.5 .5 0
]
elem =
[
1 2 3 4 5 6 7 8 9
]
edge =
[
1 5 2
2 6 3
3 7 4
4 8 1
]
bc =
[
3
-1
1
-1
]
elseif choose_mesh==2
nodes =
[
0 0 0
2 0 0
2 1 0
0 1 0
1 0 0
1 1 0
0 .5 0
.5 0 0
1 .5 0
.5 1 0
.5 .5 0
2 .5 0
1.5 1 0
1.5 0 0
1.5 .5 0
]
elem =
[
1 5 6 4 8 9 10 7 11
5 2 3 6 14 12 13 9 15
]
edge =
[
1 8 5
5 14 2
2 12 3
1 7 4
4 10 6
6 13 3
]
bc =
[
3
-1
4
-1
]
elseif choose_mesh==3
nodes =
[
0 0 0
1 0 0
1 1 0
0 1 0
1 .5 0
.5 .5 0
0 .5 0
.5 0 0
1 .25 0
.5 .25 0
0 .25 0
.5 1 0
1 .75 0
.5 .75 0
0 .75 0
]
elem =
[
1 2 5 7 8 9 6 11 10
7 5 3 4 6 13 12 15 14
]
edge =
[
1 8 2
2 9 5
5 13 3
3 12 4
4 15 7
7 11 1
]
bc =
[
4
-1
1
-1
]
elseif choose_mesh==4
nodes =
[
0 0 0
1 0 0
1 1 0
0 1 0
1 .5 0
.5 .5 0
0 .5 0
.5 0 0
1 .25 0
.5 .25 0
0 .25 0
.5 1 0
1 .75 0
.5 .75 0
0 .75 0
]
elem =
[
1 2 5 7 8 9 6 11 10
7 5 3 4 6 13 12 15 14
]
edge =
[
1 8 2
2 9 5
5 13 3
3 12 4
4 15 7
7 11 1
]
bc =
[
2
3
-1
5
6
-1
]
elseif choose_mesh==5
nodes=
[
0.0   0.0   0.0
0.0   1.0   0.0
1.0   0.0   0.0
1.0   1.0   0.0
0.0   0.5   0.0
0.5   0.0   0.0
1.0   0.5   0.0
0.5   1.0   0.0
0.5   0.5   0.0
0.0   0.25  0.0
0.0   0.75  0.0
0.25  0.0   0.0
0.75  0.0   0.0
1.0   0.25  0.0
1.0   0.75  0.0
0.25  1.0   0.0
0.75  1.0   0.0
0.25  0.5   0.0
0.5   0.75  0.0
0.25  0.75  0.0
0.75  0.5   0.0
0.75  0.75  0.0
0.5   0.25  0.0
0.25  0.25  0.0
0.75  0.25  0.0
]
elem =
[
5  9  8  2  18  19  16  11  20
9  7  4  8  21  15  17  19  22
1  6  9  5  12  23  18  10  24
6  3  7  9  13  14  21  23  25
]
edge =
[
1  5  10
5  2  11
1  6  12
6  3  13
3  7  14
7  4  15
2  8  16
8  4  17
]
bc =
[
5
6
-1
1
2
-1
]
end

max_stars_repo_name: billyon/first

max_issues_repo_name: billyon/first

max_forks_repo_name: billyon/first