text:
import Base: ==

export InlineString, InlineStringType

abstract type InlineString <: AbstractString end

for sz in (1, 4, 8, 16, 32, 64, 128, 256)
    nm = Symbol(:InlineString, max(1, sz - 1))
    @eval begin
        primitive type $nm <: InlineString $(sz * 8) end
        export $nm
    end
end

_bswap(x::T) where {T <: InlineString} = T === InlineString1 ? x : Base.bswap_int(x)

const InlineStrings = Union{InlineString1,
                            InlineString3,
                            InlineString7,
                            InlineString15,
                            InlineString31,
                            InlineString63,
                            InlineString127,
                            InlineString255}

function Base.promote_type(::Type{T}, ::Type{S}) where {T <: InlineString, S <: InlineString}
    T === InlineString1 && return S
    S === InlineString1 && return T
    T === InlineString3 && return S
    S === InlineString3 && return T
    T === InlineString7 && return S
    S === InlineString7 && return T
    T === InlineString15 && return S
    S === InlineString15 && return T
    T === InlineString31 && return S
    S === InlineString31 && return T
    T === InlineString63 && return S
    S === InlineString63 && return T
    T === InlineString127 && return S
    S === InlineString127 && return T
    return InlineString255
end

Base.promote_type(::Type{T}, ::Type{String}) where {T <: InlineString} = String

Base.widen(::Type{InlineString1}) = InlineString3
Base.widen(::Type{InlineString3}) = InlineString7
Base.widen(::Type{InlineString7}) = InlineString15
Base.widen(::Type{InlineString15}) = InlineString31
Base.widen(::Type{InlineString31}) = InlineString63
Base.widen(::Type{InlineString63}) = InlineString127
Base.widen(::Type{InlineString127}) = InlineString255
Base.widen(::Type{InlineString255}) = String

Base.ncodeunits(::InlineString1) = 1
Base.ncodeunits(x::InlineString) = Int(Base.trunc_int(UInt8, x))
Base.codeunit(::InlineString) = UInt8

Base.@propagate_inbounds function Base.codeunit(x::T, i::Int) where {T <: InlineString}
    @boundscheck checkbounds(Bool, x, i) || throw(BoundsError(x, i))
    if T === InlineString1
        return Base.bitcast(UInt8, x)
    else
        return Base.trunc_int(UInt8, Base.lshr_int(x, 8 * (sizeof(T) - i)))
    end
end

function Base.String(x::T) where {T <: InlineString}
    len = ncodeunits(x)
    out = Base._string_n(len)
    if T === InlineString1
        GC.@preserve out unsafe_store!(pointer(out), codeunit(x, 1))
        return out
    end
    ref = Ref{T}(_bswap(x))
    GC.@preserve ref out begin
        ptr = convert(Ptr{UInt8}, Base.unsafe_convert(Ptr{T}, ref))
        unsafe_copyto!(pointer(out), ptr, len)
    end
    return out
end

function Base.Symbol(x::T) where {T <: InlineString}
    ref = Ref{T}(_bswap(x))
    return ccall(:jl_symbol_n, Ref{Symbol},
        (Ref{T}, Int), ref, sizeof(x))
end

# add a codeunit to end of string method
function addcodeunit(x::T, b::UInt8) where {T <: InlineString}
    if T === InlineString1
        return x, false
    end
    len = ncodeunits(x)
    (len + 1) == sizeof(T) && return x, true
    x = Base.or_int(x, Base.shl_int(Base.zext_int(T, b), 8 * (sizeof(T) - len - 1)))
    return Base.add_int(x, Base.zext_int(T, 0x01)), false
end

# from String
InlineString1(byte::UInt8=0x00) = Base.bitcast(InlineString1, byte)
(::Type{T})() where {T <: InlineString} = Base.zext_int(T, 0x00)

function (::Type{T})(x::AbstractString) where {T <: InlineString}
    if T === InlineString1
        sizeof(x) == 1 || stringtoolong(T, sizeof(x))
        return Base.bitcast(InlineString1, codeunit(x, 1))
    elseif typeof(x) === String && sizeof(T) <= sizeof(UInt)
        len = sizeof(x)
        len < sizeof(T) || stringtoolong(T, len)
        y = GC.@preserve x unsafe_load(convert(Ptr{T}, pointer(x)))
        sz = 8 * (sizeof(T) - len)
        return Base.or_int(Base.shl_int(Base.lshr_int(_bswap(y), sz), sz), Base.zext_int(T, UInt8(len)))
    else
        len = ncodeunits(x)
        len < sizeof(T) || stringtoolong(T, len)
        y = T()
        for i = 1:len
            @inbounds y, _ = addcodeunit(y, codeunit(x, i))
        end
        return y
    end
end

@noinline stringtoolong(T, n) = throw(ArgumentError("string too large ($n) to convert to $T"))

function InlineStringType(n::Integer)
    n > 255 && stringtoolong(InlineString, n)
    return n == 1  ? InlineString1   : n < 4  ? InlineString3  :
           n < 8   ? InlineString7   : n < 16 ? InlineString15 :
           n < 32  ? InlineString31  : n < 64 ? InlineString63 :
           n < 128 ? InlineString127 : InlineString255
end

InlineString(x::AbstractString)::InlineStrings = (InlineStringType(ncodeunits(x)))(x)

# between InlineStrings
function (::Type{T})(x::S) where {T <: InlineString, S <: InlineString}
    if T === S
        return x
    elseif T === InlineString1
        sizeof(x) == 1 || stringtoolong(T, sizeof(x))
        return Base.bitcast(InlineString1, codeunit(x, 1))
    elseif sizeof(T) < sizeof(S)
        # trying to compress
        len = sizeof(x)
        len > (sizeof(T) - 1) && stringtoolong(T, len)
        y = Base.trunc_int(T, Base.lshr_int(x, 8 * (sizeof(S) - sizeof(T))))
        return Base.add_int(y, Base.zext_int(T, UInt8(len)))
    else
        # promoting smaller InlineString to larger
        if S === InlineString1
            y = Base.shl_int(Base.zext_int(T, x), 8 * (sizeof(T) - sizeof(S)))
        else
            y = Base.shl_int(Base.zext_int(T, Base.lshr_int(x, 8)), 8 * (sizeof(T) - sizeof(S) + 1))
        end
        return Base.add_int(y, Base.zext_int(T, UInt8(sizeof(x))))
    end
end

(==)(x::T, y::T) where {T <: InlineString} = Base.eq_int(x, y)
function ==(x::String, y::T) where {T <: InlineString}
    sizeof(x) == sizeof(y) || return false
    ref = Ref{T}(_bswap(y))
    return ccall(:memcmp, Cint, (Ptr{UInt8}, Ref{T}, Csize_t),
            pointer(x), ref, sizeof(x)) == 0
end
==(y::InlineString, x::String) = x == y

function Base.hash(x::T, h::UInt) where {T <: InlineString}
    h += Base.memhash_seed
    ref = Ref{T}(_bswap(x))
    return ccall(Base.memhash, UInt,
        (Ref{T}, Csize_t, UInt32),
        ref, sizeof(x), h % UInt32) + h
end

function Base.write(io::IO, x::T) where {T <: InlineString}
    ref = Ref{T}(_bswap(x))
    return GC.@preserve ref begin
        ptr = convert(Ptr{UInt8}, Base.unsafe_convert(Ptr{T}, ref))
        Int(unsafe_write(io, ptr, reinterpret(UInt, sizeof(x))))::Int
    end
end

Base.print(io::IO, s::InlineString) = (write(io, s); nothing)

function Base.isascii(x::T) where {T <: InlineString}
    if T === InlineString1
        return codeunit(x, 1) < 0x80
    end
    len = ncodeunits(x)
    x = Base.lshr_int(x, 8 * (sizeof(T) - len))
    for _ = 1:len
        Base.trunc_int(UInt8, x) >= 0x80 && return false
        x = Base.lshr_int(x, 8)
    end
    return true
end

# copy/pasted from substring.jl
function Base.reverse(s::InlineString)
    # Read characters forwards from `s` and write backwards to `out`
    out = Base._string_n(sizeof(s))
    offs = sizeof(s) + 1
    for c in s
        offs -= ncodeunits(c)
        Base.__unsafe_string!(out, c, offs)
    end
    return out
end

@inline function Base.__unsafe_string!(out, x::T, offs::Integer) where {T <: InlineString}
    n = sizeof(x)
    ref = Ref{T}(_bswap(x))
    GC.@preserve ref out begin
        ptr = convert(Ptr{UInt8}, Base.unsafe_convert(Ptr{T}, ref))
        unsafe_copyto!(pointer(out, offs), ptr, n)
    end
    return n
end

Base.string(a::InlineString) = a
Base.string(a::InlineString...) = _string(a...)
Base.string(a::BaseStrs, b::InlineString) = _string(a, b)
Base.string(a::BaseStrs, b::BaseStrs, c::InlineString) = _string(a, b, c)
@inline function _string(a::Union{BaseStrs, InlineString}...)
    n = 0
    for v in a
        if v isa Char
            n += ncodeunits(v)
        else
            n += sizeof(v)
        end
    end
    out = Base._string_n(n)
    offs = 1
    for v in a
        offs += Base.__unsafe_string!(out, v, offs)
    end
    return out
end

function Base.repeat(x::T, r::Integer) where {T <: InlineString}
    r < 0 && throw(ArgumentError("can't repeat a string $r times"))
    r == 0 && return ""
    r == 1 && return s
    n = sizeof(x)
    out = Base._string_n(n * r)
    if n == 1 # common case: repeating a single-byte string
        @inbounds b = codeunit(x, 1)
        ccall(:memset, Ptr{Cvoid}, (Ptr{UInt8}, Cint, Csize_t), out, b, r)
    else
        for i = 0:r-1
            ref = Ref{T}(_bswap(x))
            GC.@preserve ref out begin
                ptr = convert(Ptr{UInt8}, Base.unsafe_convert(Ptr{T}, ref))
                unsafe_copyto!(pointer(out, i * n + 1), ptr, n)
            end
        end
    end
    return out
end

# copy/pasted from strings/util.jl
function Base.startswith(a::T, b::Union{String, SubString{String}, InlineString}) where {T <: InlineString}
    cub = ncodeunits(b)
    ncodeunits(a) < cub && return false
    ref = Ref{T}(_bswap(a))
    return GC.@preserve ref begin
        ptr = convert(Ptr{UInt8}, Base.unsafe_convert(Ptr{T}, ref))
        if Base._memcmp(ptr, b, sizeof(b)) == 0
            nextind(a, cub) == cub + 1
        else
            false
        end
    end
end

function Base.endswith(a::T, b::Union{String, SubString{String}, InlineString}) where {T <: InlineString}
    cub = ncodeunits(b)
    astart = ncodeunits(a) - ncodeunits(b) + 1
    astart < 1 && return false
    ref = Ref{T}(_bswap(a))
    return GC.@preserve ref begin
        ptr = convert(Ptr{UInt8}, Base.unsafe_convert(Ptr{T}, ref))
        if Base._memcmp(ptr + (astart - 1), b, sizeof(b)) == 0
            thisind(a, astart) == astart
        else
            false
        end
    end
end

Base.match(r::Regex, s::InlineString, i::Integer) = match(r, String(s), i)

# the rest of these methods are copy/pasted from Base strings/string.jl file
# for efficiency
Base.@propagate_inbounds function Base.isvalid(x::InlineString, i::Int)
    @boundscheck checkbounds(Bool, x, i) || throw(BoundsError(x, i))
    return @inbounds thisind(x, i) == i
end

Base.@propagate_inbounds function Base.thisind(s::InlineString, i::Int)
    i == 0 && return 0
    n = ncodeunits(s)
    i == n + 1 && return i
    @boundscheck Base.between(i, 1, n) || throw(BoundsError(s, i))
    @inbounds b = codeunit(s, i)
    (b & 0xc0 == 0x80) & (i-1 > 0) || return i
    @inbounds b = codeunit(s, i-1)
    Base.between(b, 0b11000000, 0b11110111) && return i-1
    (b & 0xc0 == 0x80) & (i-2 > 0) || return i
    @inbounds b = codeunit(s, i-2)
    Base.between(b, 0b11100000, 0b11110111) && return i-2
    (b & 0xc0 == 0x80) & (i-3 > 0) || return i
    @inbounds b = codeunit(s, i-3)
    Base.between(b, 0b11110000, 0b11110111) && return i-3
    return i
end

Base.@propagate_inbounds function Base.nextind(s::InlineString, i::Int)
    i == 0 && return 1
    n = ncodeunits(s)
    @boundscheck Base.between(i, 1, n) || throw(BoundsError(s, i))
    @inbounds l = codeunit(s, i)
    (l < 0x80) | (0xf8 ≤ l) && return i+1
    if l < 0xc0
        i′ = @inbounds thisind(s, i)
        return i′ < i ? @inbounds(nextind(s, i′)) : i+1
    end
    # first continuation byte
    (i += 1) > n && return i
    @inbounds b = codeunit(s, i)
    b & 0xc0 ≠ 0x80 && return i
    ((i += 1) > n) | (l < 0xe0) && return i
    # second continuation byte
    @inbounds b = codeunit(s, i)
    b & 0xc0 ≠ 0x80 && return i
    ((i += 1) > n) | (l < 0xf0) && return i
    # third continuation byte
    @inbounds b = codeunit(s, i)
    ifelse(b & 0xc0 ≠ 0x80, i, i+1)
end

Base.@propagate_inbounds function Base.iterate(s::InlineString, i::Int=firstindex(s))
    (i % UInt) - 1 < ncodeunits(s) || return nothing
    b = @inbounds codeunit(s, i)
    u = UInt32(b) << 24
    Base.between(b, 0x80, 0xf7) || return reinterpret(Char, u), i+1
    return iterate_continued(s, i, u)
end

function iterate_continued(s::InlineString, i::Int, u::UInt32)
    u < 0xc0000000 && (i += 1; @goto ret)
    n = ncodeunits(s)
    # first continuation byte
    (i += 1) > n && @goto ret
    @inbounds b = codeunit(s, i)
    b & 0xc0 == 0x80 || @goto ret
    u |= UInt32(b) << 16
    # second continuation byte
    ((i += 1) > n) | (u < 0xe0000000) && @goto ret
    @inbounds b = codeunit(s, i)
    b & 0xc0 == 0x80 || @goto ret
    u |= UInt32(b) << 8
    # third continuation byte
    ((i += 1) > n) | (u < 0xf0000000) && @goto ret
    @inbounds b = codeunit(s, i)
    b & 0xc0 == 0x80 || @goto ret
    u |= UInt32(b); i += 1
@label ret
    return reinterpret(Char, u), i
end

Base.@propagate_inbounds function Base.getindex(s::InlineString, i::Int)
    b = codeunit(s, i)
    u = UInt32(b) << 24
    Base.between(b, 0x80, 0xf7) || return reinterpret(Char, u)
    return getindex_continued(s, i, u)
end

function getindex_continued(s::InlineString, i::Int, u::UInt32)
    if u < 0xc0000000
        # called from `getindex` which checks bounds
        @inbounds isvalid(s, i) && @goto ret
        Base.string_index_err(s, i)
    end
    n = ncodeunits(s)

    (i += 1) > n && @goto ret
    @inbounds b = codeunit(s, i) # cont byte 1
    b & 0xc0 == 0x80 || @goto ret
    u |= UInt32(b) << 16

    ((i += 1) > n) | (u < 0xe0000000) && @goto ret
    @inbounds b = codeunit(s, i) # cont byte 2
    b & 0xc0 == 0x80 || @goto ret
    u |= UInt32(b) << 8

    ((i += 1) > n) | (u < 0xf0000000) && @goto ret
    @inbounds b = codeunit(s, i) # cont byte 3
    b & 0xc0 == 0x80 || @goto ret
    u |= UInt32(b)
@label ret
    return reinterpret(Char, u)
end

Base.length(s::InlineString) = length_continued(s, 1, ncodeunits(s), ncodeunits(s))

Base.@propagate_inbounds function Base.length(s::InlineString, i::Int, j::Int)
    @boundscheck begin
        0 < i ≤ ncodeunits(s)+1 || throw(BoundsError(s, i))
        0 ≤ j < ncodeunits(s)+1 || throw(BoundsError(s, j))
    end
    j < i && return 0
    @inbounds i, k = thisind(s, i), i
    c = j - i + (i == k)
    length_continued(s, i, j, c)
end

@inline function length_continued(s::InlineString, i::Int, n::Int, c::Int)
    i < n || return c
    @inbounds b = codeunit(s, i)
    @inbounds while true
        while true
            (i += 1) ≤ n || return c
            0xc0 ≤ b ≤ 0xf7 && break
            b = codeunit(s, i)
        end
        l = b
        b = codeunit(s, i) # cont byte 1
        c -= (x = b & 0xc0 == 0x80)
        x & (l ≥ 0xe0) || continue

        (i += 1) ≤ n || return c
        b = codeunit(s, i) # cont byte 2
        c -= (x = b & 0xc0 == 0x80)
        x & (l ≥ 0xf0) || continue

        (i += 1) ≤ n || return c
        b = codeunit(s, i) # cont byte 3
        c -= (b & 0xc0 == 0x80)
    end
end

# parsers
# this is mostly copy-pasta from Parsers.jl main xparse function
import Parsers: SENTINEL, OK, EOF, OVERFLOW, QUOTED, DELIMITED, INVALID_QUOTED_FIELD, ESCAPED_STRING, NEWLINE, SUCCESS, peekbyte, incr!, checksentinel, checkdelim, checkcmtemptylines

@inline function Parsers.xparse(::Type{T}, source::Union{AbstractVector{UInt8}, IO}, pos, len, options::Parsers.Options{ignorerepeated, ignoreemptylines, Q, debug, S, D, DF}) where {T <: InlineString, ignorerepeated, ignoreemptylines, Q, debug, S, D, DF}
    startpos = vstartpos = vpos = pos
    sentstart = sentinelpos = 0
    code = SUCCESS
    sentinel = options.sentinel
    quoted = overflowed = false
    x = T()
    if debug
        println("parsing $T, pos=$pos, len=$len")
    end
    if Parsers.eof(source, pos, len)
        code = (sentinel === missing ? SENTINEL : OK) | EOF
        if T === InlineString1
            # InlineString1 must be exactly 1 byte, so for empty string
            # this is an "underflow" situation
            code |= OVERFLOW
        end
        @goto donedone
    end
    b = peekbyte(source, pos)
    if debug
        println("string 1) parsed: '$(escape_string(string(Char(b))))'")
    end
    # strip leading whitespace
    while b == options.wh1 || b == options.wh2
        if debug
            println("stripping leading whitespace")
        end
        x, overflowed = addcodeunit(x, b)
        pos += 1
        incr!(source)
        if Parsers.eof(source, pos, len)
            code |= EOF
            @goto donedone
        end
        b = peekbyte(source, pos)
        if debug
            println("string 2) parsed: '$(escape_string(string(Char(b))))'")
        end
    end
    # check for start of quoted field
    if Q
        quoted = b == options.oq
        if quoted
            if debug
                println("detected open quote character")
            end
            code = QUOTED
            x = T() # start our parsed value back over
            pos += 1
            vstartpos = pos
            incr!(source)
            if Parsers.eof(source, pos, len)
                code |= INVALID_QUOTED_FIELD
                @goto donedone
            end
            b = peekbyte(source, pos)
            if debug
                println("string 3) parsed: '$(escape_string(string(Char(b))))'")
            end
            # ignore whitespace within quoted field
            while b == options.wh1 || b == options.wh2
                if debug
                    println("stripping whitespace within quoted field")
                end
                x, overflowed = addcodeunit(x, b)
                pos += 1
                incr!(source)
                if Parsers.eof(source, pos, len)
                    code |= INVALID_QUOTED_FIELD | EOF
                    @goto donedone
                end
                b = peekbyte(source, pos)
                if debug
                    println("string 4) parsed: '$(escape_string(string(Char(b))))'")
                end
            end
        end
    end
    # check for sentinel values if applicable
    if sentinel !== nothing && sentinel !== missing
        if debug
            println("checking for sentinel value")
        end
        sentstart = pos
        sentinelpos = checksentinel(source, pos, len, sentinel, debug)
    end
    vpos = pos
    if Q
        # for quoted fields, find the closing quote character
        if quoted
            if debug
                println("looking for close quote character")
            end
            same = options.cq == options.e
            while true
                vpos = pos
                pos += 1
                incr!(source)
                if same && b == options.e
                    if Parsers.eof(source, pos, len)
                        code |= EOF
                        @goto donedone
                    elseif peekbyte(source, pos) != options.cq
                        break
                    end
                    code |= ESCAPED_STRING
                    b = peekbyte(source, pos)
                    pos += 1
                    incr!(source)
                elseif b == options.e
                    if Parsers.eof(source, pos, len)
                        code |= INVALID_QUOTED_FIELD | EOF
                        @goto donedone
                    end
                    code |= ESCAPED_STRING
                    b = peekbyte(source, pos)
                    pos += 1
                    incr!(source)
                elseif b == options.cq
                    if Parsers.eof(source, pos, len)
                        code |= EOF
                        @goto donedone
                    end
                    break
                end
                if Parsers.eof(source, pos, len)
                    code |= INVALID_QUOTED_FIELD | EOF
                    @goto donedone
                end
                x, overflowed = addcodeunit(x, b)
                b = peekbyte(source, pos)
                if debug
                    println("string 9) parsed: '$(escape_string(string(Char(b))))'")
                end
            end
            b = peekbyte(source, pos)
            if debug
                println("string 10) parsed: '$(escape_string(string(Char(b))))'")
            end
            # ignore whitespace after quoted field
            while b == options.wh1 || b == options.wh2
                if debug
                    println("stripping trailing whitespace after close quote character")
                end
                pos += 1
                incr!(source)
                if Parsers.eof(source, pos, len)
                    code |= EOF
                    @goto donedone
                end
                b = peekbyte(source, pos)
                if debug
                    println("string 11) parsed: '$(escape_string(string(Char(b))))'")
                end
            end
        end
    end
    if options.delim !== nothing
        delim = options.delim
        quo = Int(!quoted)
        # now we check for a delimiter; if we don't find it, keep parsing until we do
        if debug
            println("checking for delimiter: pos=$pos")
        end
        while true
            if !ignorerepeated
                if delim isa UInt8
                    if b == delim
                        pos += 1
                        incr!(source)
                        code |= DELIMITED
                        @goto donedone
                    end
                else
                    predelimpos = pos
                    pos = checkdelim(source, pos, len, delim)
                    if pos > predelimpos
                        # found the delimiter we were looking for
                        code |= DELIMITED
                        @goto donedone
                    end
                end
            else
                if delim isa UInt8
                    matched = false
                    matchednewline = false
                    while true
                        if b == delim
                            matched = true
                            code |= DELIMITED
                            pos += 1
                            incr!(source)
                        elseif !matchednewline && b == UInt8('\n')
                            matchednewline = matched = true
                            pos += 1
                            incr!(source)
                            pos = checkcmtemptylines(source, pos, len, options)
                            code |= NEWLINE | ifelse(Parsers.eof(source, pos, len), EOF, SUCCESS)
                        elseif !matchednewline && b == UInt8('\r')
                            matchednewline = matched = true
                            pos += 1
                            incr!(source)
                            if !Parsers.eof(source, pos, len) && peekbyte(source, pos) == UInt8('\n')
                                pos += 1
                                incr!(source)
                            end
                            pos = checkcmtemptylines(source, pos, len, options)
                            code |= NEWLINE | ifelse(Parsers.eof(source, pos, len), EOF, SUCCESS)
                        else
                            break
                        end
                        if Parsers.eof(source, pos, len)
                            @goto donedone
                        end
                        b = peekbyte(source, pos)
                        if debug
                            println("14) parsed: '$(escape_string(string(Char(b))))'")
                        end
                    end
                    if matched
                        @goto donedone
                    end
                else
                    matched = false
                    matchednewline = false
                    while true
                        predelimpos = pos
                        pos = checkdelim(source, pos, len, delim)
                        if pos > predelimpos
                            matched = true
                            code |= DELIMITED
                        elseif !matchednewline && b == UInt8('\n')
                            matchednewline = matched = true
                            pos += 1
                            incr!(source)
                            pos = checkcmtemptylines(source, pos, len, options)
                            code |= NEWLINE | ifelse(Parsers.eof(source, pos, len), EOF, SUCCESS)
                        elseif !matchednewline && b == UInt8('\r')
                            matchednewline = matched = true
                            pos += 1
                            incr!(source)
                            if !Parsers.eof(source, pos, len) && peekbyte(source, pos) == UInt8('\n')
                                pos += 1
                                incr!(source)
                            end
                            pos = checkcmtemptylines(source, pos, len, options)
                            code |= NEWLINE | ifelse(Parsers.eof(source, pos, len), EOF, SUCCESS)
                        else
                            break
                        end
                        if Parsers.eof(source, pos, len)
                            @goto donedone
                        end
                        b = peekbyte(source, pos)
                        if debug
                            println("14) parsed: '$(escape_string(string(Char(b))))'")
                        end
                    end
                    if matched
                        @goto donedone
                    end
                end
            end
            # didn't find delimiter, but let's check for a newline character
            if b == UInt8('\n')
                pos += 1
                incr!(source)
                pos = checkcmtemptylines(source, pos, len, options)
                code |= NEWLINE | ifelse(Parsers.eof(source, pos, len), EOF, SUCCESS)
                @goto donedone
            elseif b == UInt8('\r')
                pos += 1
                incr!(source)
                if !Parsers.eof(source, pos, len) && peekbyte(source, pos) == UInt8('\n')
                    pos += 1
                    incr!(source)
                end
                pos = checkcmtemptylines(source, pos, len, options)
                code |= NEWLINE | ifelse(Parsers.eof(source, pos, len), EOF, SUCCESS)
                @goto donedone
            end
            # didn't find delimiter nor newline, so increment and check the next byte
            x, overflowed = addcodeunit(x, b)
            pos += 1
            vpos += quo
            incr!(source)
            if Parsers.eof(source, pos, len)
                code |= EOF
                @goto donedone
            end
            b = peekbyte(source, pos)
        end
    else
        # no delimiter, so read until EOF
        while !Parsers.eof(source, pos, len)
            b = peekbyte(source, pos)
            x, overflowed = addcodeunit(x, b)
            pos += 1
            vpos += 1
            incr!(source)
        end
        code |= EOF
    end

@label donedone
    if sentinel !== nothing && sentinel !== missing && sentstart == vstartpos && sentinelpos == vpos
        # if we matched a sentinel value that was as long or longer than our type value
        code |= SENTINEL
    elseif sentinel === missing && vstartpos == vpos
        code |= SENTINEL
    else
        code |= OK
        if T === InlineString1
            if (vpos - vstartpos) != 1
                overflowed = true
            else
                Parsers.fastseek!(source, vstartpos)
                x = InlineString1(peekbyte(source, vstartpos))
                Parsers.fastseek!(source, pos - 1)
            end
        end
    end
    if overflowed
        code |= OVERFLOW
    end
    if debug
        println("finished parsing: $(codes(code))")
    end
    return x, code, Int64(vstartpos), Int64(vpos - vstartpos), Int64(pos - startpos)
end


text:
include("AbstractInterpolationGrid.jl")
include("RegularGrid.jl")
include("IrregularGrid.jl")

ext: jl


text:
# This file was generated, do not modify it. # hide
using Gnuplot
@gp 1:10
#save("plt2_ex1.gp")
#save(term="cairolatex standalone", output="plt2_ex1.tex")
#save(term="pdf", output="plt2_ex1.pdf")
save(term="png", output="plt2_ex1.png")


text:
using DMUStudent
using Test
using DataFrames
import JSON

@testset "email check" begin
    @test DMUStudent.is_identikey_colorado_email("zasu3213@colorado.edu")
    @test !DMUStudent.is_identikey_colorado_email("zachary.sunberg@colorado.edu")
    @test !DMUStudent.is_identikey_colorado_email("sunbergzach@gmail.com")
end

# include("hw1.jl")
# include("hw2.jl")
# include("hw3.jl")
# include("hw4.jl")
# include("hw5.jl")
include("hw6.jl")


text:
export track, topN
using Combinatorics
using DataStructures

abstract type Tracking end
@with_kw struct SimpleTracking <: Tracking
  tc::typeof(1.0s) = 1s
end
Tracking(::Val{:simple};params...) = SimpleTracking(;params...)

# TODO: this makes me realize I can make a more general metadata array that
# isn't an axis array at all
const WithAxes{Ax} = AxisArray{<:Any,<:Any,<:Any,Ax}
const SourceTracking =
  MetaArray{<:WithAxes{<:Tuple{Axis{:scale},Axis{:freq},
                               Axis{:component},Axis{:time}}}}
component_means(S::SourceTracking) = vec(mean(S,dims=[1,2,4]))
function component_means(S::SubArray{<:Any,<:Any,<:SourceTracking})
  vec(mean(S,dims=[1,2,4]))
end

function track(C::Coherence;method=:simple,progressbar=true,params...)
  method = Tracking(C,Val{method}();params...)
  track(C,method,progressbar)
end

function topN(by,N,xs)
  state = start(xs)
  if done(xs,state)
    error("Empty iterable.")
  end
  x, state = next(xs,state)
  by_x = by(x)
  queue = PriorityQueue{typeof(by_x),typeof(x)}()
  enqueue!(queue,x,by_x)

  while !done(xs,state)
    x, state = next(xs,state)
    by_x = by(x)
    enqueue!(queue,x,by_x)
    if length(queue) > N
      dequeue!(queue)
    end
  end

  result = Array{eltype(xs)}(length(queue))
  i = 0
  while !isempty(queue)
    result[i += 1] = dequeue!(queue)
  end
  result
end

function maximumby(by,xs)
  itr = iterate(xs)
  if itr == nothing
    error("Empty iterator.")
  end

  x,state = itr
  result = x
  maxval = by(x)
  itr = iterate(xs,state)
  while itr !== nothing
    x, state = itr
    val = by(x)
    if isless(maxval,val)
      result = x
      maxval = val
    end
    itr = iterate(xs,state)
  end
  result, maxval
end

function bestordering(x,y)
  @assert length(x) == length(y)
  @assert length(x) <= 6 "No more than 6 components supported"*
    " (due to combinatorial explosion of possible solutions.)"
  # minimize cosine difference
  diff(x,y) = norm(vec(x) .- vec(y),1)

  maximumby(permutations(1:length(y))) do order
    -sum(diff(xi,yi) for (xi,yi) in zip(x,y[order]))
  end
end

function track_progress(progressbar,n,name)
  if progressbar
    Progress(n,desc="Source Tracking ($name): ")
  end
end

max_stars_repo_name: JuliaTagBot/AuditoryBistabilityLE

max_issues_repo_name: JuliaTagBot/AuditoryBistabilityLE

max_forks_repo_name: JuliaTagBot/AuditoryBistabilityLE


text:
"""
    Exo2SplitBinion{T <: Real} = Parabolic{Exo1SplitBinion{T}}

A 2-exo-split-binion is a parabolic Cayley-Dickson construct with 1-exo-split-binions.
"""
const Exo2SplitBinion{T <: Real} = Parabolic{Exo1SplitBinion{T}}

function Exo2SplitBinion(z::Exo1SplitBinion{T}) where T <: Real
    Exo2SplitBinion{T}(z, zero(Exo1SplitBinion{T}))
end

function Exo2SplitBinion(z::SplitBinion{T}) where T <: Real
    Exo2SplitBinion{T}(
        Exo1SplitBinion{T}(z, zero(SplitBinion{T})),
        zero(Exo1SplitBinion{T})
    )
end

function Exo2SplitBinion(a::T, b::T, c::T, d::T, f::T, g::T, h::T, j::T) where T <: Real
    Exo2SplitBinion{T}(
        Exo1SplitBinion{T}(
            SplitBinion{T}(a, b),
            SplitBinion{T}(c, d)
        ),
        Exo1SplitBinion{T}(
            SplitBinion{T}(f, g),
            SplitBinion{T}(h, j)
        )
    )
end

function Exo2SplitBinion(a::T, b::T, c::T, d::T, f::T, g::T, h::T) where T <: Real
    Exo2SplitBinion{T}(
        Exo1SplitBinion{T}(
            SplitBinion{T}(a, b),
            SplitBinion{T}(c, d)
        ),
        Exo1SplitBinion{T}(
            SplitBinion{T}(f, g),
            SplitBinion{T}(h, zero(T))
        )
    )
end

function Exo2SplitBinion(a::T, b::T, c::T, d::T, f::T, g::T) where T <: Real
    Exo2SplitBinion{T}(
        Exo1SplitBinion{T}(
            SplitBinion{T}(a, b),
            SplitBinion{T}(c, d)
        ),
        Exo1SplitBinion{T}(
            SplitBinion{T}(f, g),
            zero(SplitBinion{T})
        )
    )
end

function Exo2SplitBinion(a::T, b::T, c::T, d::T, f::T) where T <: Real
    Exo2SplitBinion{T}(
        Exo1SplitBinion{T}(
            SplitBinion{T}(a, b),
            SplitBinion{T}(c, d)
        ),
        Exo1SplitBinion{T}(
            SplitBinion{T}(f, zero(T)),
            zero(SplitBinion{T})
        )
    )
end

function Exo2SplitBinion(a::T, b::T, c::T, d::T) where T <: Real
    Exo2SplitBinion{T}(
        Exo1SplitBinion{T}(
            SplitBinion{T}(a, b),
            SplitBinion{T}(c, d)
        ),
        zero(Exo1SplitBinion{T})
    )
end

function Exo2SplitBinion(a::T, b::T, c::T) where T <: Real
    Exo2SplitBinion{T}(
        Exo1SplitBinion{T}(
            SplitBinion{T}(a, b),
            SplitBinion{T}(c, zero(T))
        ),
        zero(Exo1SplitBinion{T})
    )
end

function Exo2SplitBinion(a::T, b::T) where T <: Real
    Exo2SplitBinion{T}(
        Exo1SplitBinion{T}(
            SplitBinion{T}(a, b),
            zero(SplitBinion{T})
        ),
        zero(Exo1SplitBinion{T})
    )
end

function Exo2SplitBinion(a::T) where T <: Real
    Exo2SplitBinion{T}(
        Exo1SplitBinion{T}(
            SplitBinion{T}(a, zero(T)),
            zero(SplitBinion{T})
        ),
        zero(Exo1SplitBinion{T})
    )
end

function Exo2SplitBinion(a::Real, b::Real)
    Exo2SplitBinion(promote(a, b)...)
end

function Exo2SplitBinion(a::Real, b::Real, c::Real)
    Exo2SplitBinion(promote(a, b, c)...)
end

function Exo2SplitBinion(a::Real, b::Real, c::Real, d::Real)
    Exo2SplitBinion(promote(a, b, c, d)...)
end

function Exo2SplitBinion(a::Real, b::Real, c::Real, d::Real, f::Real)
    Exo2SplitBinion(promote(a, b, c, d, f)...)
end

function Exo2SplitBinion(a::Real, b::Real, c::Real, d::Real, f::Real, g::Real)
    Exo2SplitBinion(promote(a, b, c, d, f, g)...)
end

function Exo2SplitBinion(a::Real, b::Real, c::Real, d::Real, f::Real, g::Real, h::Real)
    Exo2SplitBinion(promote(a, b, c, d, f, g, h)...)
end

function Exo2SplitBinion(a::Real, b::Real, c::Real, d::Real, f::Real, g::Real, h::Real, j::Real)
    Exo2SplitBinion(promote(a, b, c, d, f, g, h, j)...)
end

function Exo2SplitBinion(z::Exo2Real{T}) where T <: Real
    Exo2SplitBinion(z.l.l, zero(T), z.l.r, zero(T), z.r.l, zero(T), z.r.r)
end

"""
    unreal(z::Exo2SplitBinion)

The unreal part of a 2-exo-split-binion is a 7-vector.
"""
function unreal(z::Exo2SplitBinion)
    vcat(unreal(z.l), array(z.r))
end

function show(io::IO, z::Exo2SplitBinion)
    print(io, "[1: ")
    print(io, z.l.l.l)
    print(io, ", s: ")
    print(io, z.l.l.r)
    print(io, ", a: ")
    print(io, z.l.r.l)
    print(io, ", sa: ")
    print(io, z.l.r.r)
    print(io, ", b: ")
    print(io, z.r.l.l)
    print(io, ", sb: ")
    print(io, z.r.l.r)
    print(io, ", ab: ")
    print(io, z.r.r.l)
    print(io, ", (sa)b: ")
    print(io, z.r.r.r)
    print(io, "]")
end

function random(::Type{Exo2SplitBinion{T}}) where T <: Real
    Exo2SplitBinion{T}(random(Exo1SplitBinion{T}), random(Exo1SplitBinion{T}))
end


text:
function compute()
	"""
	s=Array{Char,1}()
	for i=1:1000000
		for c in string(i)
			push!(s,c)
		end
	end
	"""
	s=String[];
	for i=1:1000000
		push!(s,string(i))
	end
	s=string(s...)

	ans=1
	for i=1:6
		aa=parse(Int,s[10^i])
		ans *= aa
	end
	return ans

end

println(compute())


text:
using VersatileHDPMixtureModels

nips_path = "docword.nips.txt"

function text_file_to_dict(path)
	data_dict = Dict()
	cur_doc = 0
	open(path) do file
		cur_vec =[]
		for ln in eachline(file)
			items = [parse(Int64,x) for x in split(ln)]
			if length(items) < 3
				continue
			end
			if cur_doc != items[1]
				data_dict[cur_doc] = Float32.(cur_vec)[:,:]'
				cur_vec = []
				cur_doc = items[1]
			end
			for i=1:items[3]
				push!(cur_vec,items[2])
			end
		end
		data_dict[cur_doc] = Float32.(cur_vec)[:,:]'
	end

	delete!(data_dict,0)
	return data_dict
end

nips_data = text_file_to_dict(nips_path)
nips_gprior = topic_modeling_hyper(Float64.(ones(12419))*0.1)
### Unlike the HDP, here to promote splits you need to increase both γ and α.
model = hdp_fit(nips_data,100.0,100.0,nips_gprior,80,1,10)
avg_word = VersatileHDPMixtureModels.calc_avg_word(model[1])


text:
module JacobianUtils

using SparseArrays
using LinearAlgebra
using Statistics
using Rotations
using StaticArrays
using HDF5
using Arpack
using DelimitedFiles

using LinAlgUtils

export BAProblemInstance
export readproblem, splitjacobian, estimate_parameter_scales
export estimate_camera_center_scale
export compute_measurement_covariance
export form_nullspace, build_schur_complement_matrix
export compute_covariance_evecs
export compute_evecs_deflation, compute_evecs_projection, compute_evecs_dropping

mutable struct BAProblemInstance
    num_cameras::Int
    num_points::Int
    num_observations::Int
    camera_parameters::Array{Float64}
    point_parameters::Array{Float64}
    residuals::Array{Float64}
    jacobian::SparseMatrixCSC{Float64,Int}
    scale_matrix::Diagonal{Float64}
end

"""
    readproblem(filename::String)

Unpack a HDF5 file containing a bundle adjustment problem.
"""
function readproblem(filename::String) :: BAProblemInstance
    num_cameras::Int                  = h5read(filename, "num_cameras")
    num_points::Int                   = h5read(filename, "num_points")
    num_observations::Int             = h5read(filename, "num_observations")

    camera_parameters::Array{Float64} = h5read(filename, "camera_parameters")
    point_parameters::Array{Float64}  = h5read(filename, "point_parameters")
    residuals::Array{Float64}         = h5read(filename, "residuals")

    rows::Array{Int}                  = h5read(filename, "jacobian_sparse_row")
    cols::Array{Int}                  = h5read(filename, "jacobian_sparse_column")
    vals::Array{Float64}              = h5read(filename, "jacobian_sparse_value")
    jacobian = sparse(rows .+ 1, cols .+ 1, vals)

    scale_matrix = Diagonal(ones(Float64, 6*num_cameras))

    return BAProblemInstance(
        num_cameras, num_points, num_observations,
        camera_parameters, point_parameters,
        residuals, jacobian,
        scale_matrix
    )
end

"""
Estimate the size of a SfM problem. Compute the distances from camera centers
to the mean camera location, and return the median of those distances.
"""
function estimate_camera_center_scale(problem::BAProblemInstance) :: Float64
    camera_centers = zeros(Float64, problem.num_cameras, 3)
    for k = 0:problem.num_cameras-1
        rot = problem.camera_parameters[9*k+1:9*k+3]
        θ = norm(rot)
        rot_aa = AngleAxis(θ, rot[1], rot[2], rot[3])
        t = problem.camera_parameters[9*k+4:9*k+6]
        camera_centers[k+1,:] = -rot_aa * t
    end

    average_cam_center = mean(camera_centers, dims=1)
    return median(sqrt.(sum((camera_centers .- average_cam_center).^2, dims=2)))
end

"""
Estimate the size of a SfM problem. Compute mean camera orientation, and
compute distances from that. Return the median orientation distance
"""
function estimate_camera_rot_scale(problem::BAProblemInstance) :: Float64
    # unpack the camera rotations from the problem params
    camera_rots = Vector{RotMatrix{3, Float64}}(undef, problem.num_cameras)
    for k = 0:problem.num_cameras-1
        rot = problem.camera_parameters[9*k+1:9*k+3]
        θ = norm(rot)
        camera_rots[k+1] = RotMatrix(AngleAxis(θ, rot[1], rot[2], rot[3]))
    end

    # compute the mean orientation
    Ravg = average_rotmat_l2(camera_rots)

    # compute distances to the mean rotation
    distances = Vector{Float64}(undef, problem.num_cameras)
    for k = 1:problem.num_cameras
        distances[k] = rotation_angle(Ravg' * camera_rots[k])
    end
    return median(distances)
end

"""
Split the jacobian into two. Jc are the columns representing camera params
and Jp are the columns representing point params. Discard all columns for
camera intrinsics.
"""
function splitjacobian(problem::BAProblemInstance) :: Tuple{SparseMatrixCSC{Float64,Int},SparseMatrixCSC{Float64,Int}}
    mask_cameras = falses(problem.jacobian.n)
    mask_cameras[1:9*problem.num_cameras] .= true

    mask_intrinsics = falses(problem.jacobian.n)
    for i in range(7, stop=9*problem.num_cameras, step=9)
        mask_intrinsics[i]   = true
        mask_intrinsics[i+1] = true
        mask_intrinsics[i+2] = true
    end

    Jc = problem.jacobian[:, mask_cameras .& .~mask_intrinsics]
    Jp = problem.jacobian[:, .~mask_cameras]
    return Jc, Jp
end

"""
Estimate the relative scales of rotations and translations, and save it as
a Diagonal matrix.
"""
function estimate_parameter_scales(problem::BAProblemInstance) :: Diagonal{Float64}
    scale_r = estimate_camera_rot_scale(problem)
    scale_x = estimate_camera_center_scale(problem)
    scale_factors = Vector{Float64}(undef, 6*problem.num_cameras)
    for i = 1:problem.num_cameras
        scale_factors[6*i-5:6*i-3] .= scale_r
        scale_factors[6*i-2:6*i-0] .= scale_x
    end
    camera_scale_matrix = Diagonal(scale_factors)
    problem.scale_matrix = camera_scale_matrix
    return camera_scale_matrix
end

"""
Estimate the covariance in the residual noise from the values of the
residuals. (This is very coarse.) Return a sparse diagonal matrix.

See Michal Polic's thesis (pg 10) for an overview of these methods.
method=
  :bishop: use diag(abs(residuals)) as an estimate of measurement cov.
  :lhuillier: Use a dof-adjusted norm of the residuals. Assume iid.
"""
function compute_measurement_covariance(
        problem::BAProblemInstance;
        method=:bishop) :: Diagonal{Float64}

    if method == :bishop
        return Diagonal(abs.(problem.residuals))

    elseif method == :lhuillier
        num_camera_params = 9
        gauge_rank = 7
        k = problem.num_observations
        n = problem.num_cameras
        σ = norm(problem.residuals) / (2*k - num_camera_params * n - gauge_rank)

        return Diagonal(σ .* ones(2*k))
        # return σ .* sparse(I, 2*k, 2*k)
    else
        throw(ArgumentError("method=$method is not supported"))
    end
end

"""
Create a num_params x 7 matrix that represents the kernel of a similarity
transformation. Vectors orthogonal to the columns of this matrix are
normal to the gauge subspace.
"""
function form_nullspace(problem::BAProblemInstance)
    num_params = 6 * problem.num_cameras
    nullspace = Array{Float64, 2}(undef, 6*problem.num_cameras, 7)
    for k in range(0, stop=problem.num_cameras-1)
        nullspace[6*k+1:6*k+6, 1:6] = Matrix{Float64}(I, 6, 6)
        # old code: this is a bug!
        # nullspace[6*k+1:6*k+3, 7] = problem.camera_parameters[9*k+1:9*k+3]
        # nullspace[6*k+4:6*k+6, 7] = zeros(Float64, 3, 1)
        rot = problem.camera_parameters[9*k+1:9*k+3]
        θ = norm(rot)
        rot_aa = AngleAxis(θ, rot[1], rot[2], rot[3])
        t = problem.camera_parameters[9*k+4:9*k+6]

        nullspace[6*k+1:6*k+3, 7] = zeros(Float64, 3, 1)
        nullspace[6*k+4:6*k+6, 7] = -rot_aa * t
    end

    # orthonormalize it
    fact = qr(problem.scale_matrix * nullspace)
    return Matrix(fact.Q)
end

"""
The Schur complement matrix is defined as D - (B' * A^-1 * B) where
    [A  B ] = J' Σ J
    [B' D ]
and the Jacobian J is partitioned: J = [Jp, Jc].

This function explicitly builds the Schur complement matrix as a dense matrix.
"""
function build_schur_complement_matrix(
        Jc :: SparseMatrixCSC{Float64,Int64},
        Jp :: SparseMatrixCSC{Float64,Int64},
        Sc :: Diagonal{Float64},
        Σ_input :: Diagonal{Float64}) :: Matrix{Float64}

    Σ_input_inv = inv(Σ_input)

    N  :: Int64 = size(Jc, 1)    # num rows
    Mp :: Int64 = size(Jp, 2)  # num point cols
    Mc :: Int64 = size(Jc, 2)  # num cam cols

    size(Jp,          1) == N || throw(DimensionMismatch())
    size(Σ_input_inv, 1) == N || throw(DimensionMismatch())
    size(Σ_input_inv, 2) == N || throw(DimensionMismatch())

    JpT = transpose(Jp)
    JcT = transpose(Jc)
    A = JpT * Σ_input_inv * Jp             # Mp x Mp
    B = JpT * Σ_input_inv * Jc             # Mp x Mc
    D = JcT * Σ_input_inv * Jc             # Mc x Mc

    # The A matrix is block-diagonal with 3x3 blocks
    # Create Ainv using a specialized 3x3 inversion function.
    # No nice way to reform a sparse matrix...
    num_pts = convert(Int64, Mp / 3)
    row = Vector{Int64}(undef, 9 * num_pts)
    col = Vector{Int64}(undef, 9 * num_pts)
    val = Vector{Float64}(undef, 9 * num_pts)
    for i = 1:num_pts
        # get the block and invert it
        Aidx = 3 * (i-1) + 1
        block = invertSym3x3(Array(A[Aidx:Aidx+2, Aidx:Aidx+2]))

        # stuff into the sparse storage (will this compile better?)
        Lidx = 9 * (i-1) + 1
        row[Lidx + 0] = Aidx + 0
        row[Lidx + 1] = Aidx + 0
        row[Lidx + 2] = Aidx + 0
        row[Lidx + 3] = Aidx + 1
        row[Lidx + 4] = Aidx + 1
        row[Lidx + 5] = Aidx + 1
        row[Lidx + 6] = Aidx + 2
        row[Lidx + 7] = Aidx + 2
        row[Lidx + 8] = Aidx + 2

        col[Lidx + 0] = Aidx + 0
        col[Lidx + 1] = Aidx + 1
        col[Lidx + 2] = Aidx + 2
        col[Lidx + 3] = Aidx + 0
        col[Lidx + 4] = Aidx + 1
        col[Lidx + 5] = Aidx + 2
        col[Lidx + 6] = Aidx + 0
        col[Lidx + 7] = Aidx + 1
        col[Lidx + 8] = Aidx + 2

        val[Lidx + 0] = block[1,1]
        val[Lidx + 1] = block[1,2]
        val[Lidx + 2] = block[1,3]
        val[Lidx + 3] = block[2,1]
        val[Lidx + 4] = block[2,2]
        val[Lidx + 5] = block[2,3]
        val[Lidx + 6] = block[3,1]
        val[Lidx + 7] = block[3,2]
        val[Lidx + 8] = block[3,3]
    end
    Ainv = sparse(row, col, val, Mp, Mp)

    # Compose the final operator: D - B'A^{-1}B
    Z = D - B' * Ainv * B

    return Matrix(Sc * Z * Sc)
end

"""
    compute_covariance_evecs

Compute the eigenvectors corresponding to the largest eigenvalues of the
covariance matrix.

Params:
  problem::BAProblemInstance - the problem instance
  num_vecs::Int - the number of eigenvectors to compute
Keywords:
  method:
    :projection - project the gauge space to null. Run ARPACK and discard
        the first few vectors.
    :deflation (default) - use Golub deflation to collapse the gauge space.
        Run ARPACK and then re-inflate the results.
    :dropping - just drop the lowest 7 eigenvectors
  debug::bool - display debugging output or not
"""
function compute_covariance_evecs(
    problem::BAProblemInstance,
    num_vecs::Int;
    method=:deflation::Symbol,
    debug=false::Bool
    )

    if method == :projection
        return compute_evecs_projection(problem, num_vecs, debug=debug)
    elseif method == :deflation
        return compute_evecs_deflation(problem, num_vecs, debug=debug)
    else
        return compute_evecs_dropping(problem, num_vecs, debug=debug)
    end
end

function compute_evecs_projection(
    problem::BAProblemInstance,
    num_vecs::Int;
    debug=false::Bool
    ) :: Tuple{Vector{Float64}, Matrix{Float64}, Vector{Float64}}

    Jc, Jp = splitjacobian(problem)
    Sc = estimate_parameter_scales(problem)
    Σ = compute_measurement_covariance(problem, method=:bishop)
    Z = build_schur_complement_matrix(Jc, Jp, Sc, Σ)
    N = form_nullspace(problem)

    # ensure symmetry
    Z = (Z + Z') / 2

    # project Z off of N
    P = Matrix{Float64}(I, size(Z, 1), size(Z, 1)) - N * N'
    PZP = Hermitian(P * Z * P)

    # now run ARPACK to get the eigensystem
    λ, V, nconv, niter, nmult, resid = eigs(PZP, nev=num_vecs+size(N,2), which=:SM, ritzvec=true, maxiter=1000)
    i = sortperm(λ)
    λ = λ[i]
    V = V[:,i]

    if debug
        println("Report: arpack + projection")
        println("  null eigs: $(λ[1:size(N, 2)])")
        println("  real eigs: $(λ[size(N,2)+1:end])")
        println("  nconv: $nconv")
        println("  niter: $niter")
        println("  nmult: $nmult")
        evec_accuracy_report(PZP, V, λ, constraints=N)
    end
    residuals = eigensystemresiduals(P*Z*P, V[:,size(N,2)+1:end], λ[size(N,2)+1:end])

    # Reverse the parameter scaling transformation
    V = problem.scale_matrix * V[:,size(N,2)+1:end]
    V = V ./ sqrt.(sum(abs2.(V), dims=1))

    return λ[size(N,2)+1:end].^-1, V, residuals
end

function compute_evecs_dropping(
    problem::BAProblemInstance,
    num_vecs::Int;
    debug=false::Bool
    ) :: Tuple{Vector{Float64}, Matrix{Float64}, Vector{Float64}}

    Jc, Jp = splitjacobian(problem)
    Sc = estimate_parameter_scales(problem)
    Σ = compute_measurement_covariance(problem, method=:bishop)
    Z = build_schur_complement_matrix(Jc, Jp, Sc, Σ)
    N = form_nullspace(problem) # only really using this for diagnostics

    # ensure symmetry
    Z = (Z + Z') / 2

    # now run ARPACK to get the eigensystem
    λ, V, nconv, niter, nmult, resid = eigs(Z, nev=num_vecs+size(N,2), which=:SM, ritzvec=true, maxiter=1000)
    i = sortperm(λ)
    λ = λ[i]
    V = V[:,i]

    if debug
        println("Report: arpack + projection")
        println("  null eigs: $(λ[1:size(N, 2)])")
        println("  real eigs: $(λ[size(N,2)+1:end])")
        println("  nconv: $nconv")
        println("  niter: $niter")
        println("  nmult: $nmult")
        evec_accuracy_report(Z, V, λ, constraints=N)
        err = grassmanian_metric_chordal(V[:,1:7], N)
        println("null space grassmanian chordal err:", err)
    end
    residuals = eigensystemresiduals(Z, V[:,size(N,2)+1:end], λ[size(N,2)+1:end])

    # Reverse the parameter scaling transformation
    V = problem.scale_matrix * V[:,size(N,2)+1:end]
    V = V ./ sqrt.(sum(abs2.(V), dims=1))

    return λ[size(N,2)+1:end].^-1, V, residuals
end

"""
Use ARPACK to solve a constrained eigenvalue problem, this time using Golub
deflation to do hard constraint enforcement.
"""
function compute_evecs_deflation(
    problem::BAProblemInstance,
    num_vecs::Int;
    debug=false::Bool
    ) :: Tuple{Vector{Float64}, Matrix{Float64}, Vector{Float64}}

    Jc, Jp = splitjacobian(problem)
    Sc = estimate_parameter_scales(problem)
    Σ = compute_measurement_covariance(problem, method=:bishop)
    Z = build_schur_complement_matrix(Jc, Jp, Sc, Σ)
    N = form_nullspace(problem)

    # ensure symmetry
    Z = (Z + Z') / 2

    # deflate the schur complement wrt the constraint space
    m = size(N, 1)
    p = size(N, 2)
    fact = qr(N)
    Q = fact.Q * Matrix(I, m, m) # get the full Q from the factorization
    Q2 = Q[:,p+1:end]
    G = Symmetric(Q2' * Z * Q2)

    # now run ARPACK to get the eigensystem
    λ, V, nconv, niter, nmult, resid = eigs(G, nev=num_vecs, which=:SM, ritzvec=true, maxiter=1000)

    # reinflate the answer
    V = Q2 * V

    P = Matrix(I, m, m) - N * N'
    if debug
        println("Report: arpack + golub deflation")
        println("  eigs: $λ")
        println("  nconv: $nconv")
        println("  niter: $niter")
        println("  nmult: $nmult")
        evec_accuracy_report(P*Z*P, V, λ, constraints=N)
    end
    residuals = eigensystemresiduals(P*Z*P, V, λ)

    # check convergences?
    e = eps(real(eltype(G)))/2
    f = cholesky(G)
    v = V[:,1]
    x = f / v
    println("|Ainv*v - λinv*v|:" , norm(x - v./λ[1]))

    # Reverse the parameter scaling transformation
    V = problem.scale_matrix * V
    V = V ./ sqrt.(sum(abs2.(V), dims=1))

    return λ.^-1, V, residuals
end

end # module

max_stars_repo_name: wilsonkl/sfmflex-release

max_issues_repo_name: wilsonkl/sfmflex-release

max_forks_repo_name: wilsonkl/sfmflex-release


text:
export sosbuildsequence

function candidates(s::AbstractDiscreteSwitchedSystem, l, curstate)
    switchings(s, l, curstate, false)
end

function best_dynamic(s::AbstractSwitchedSystem, d, μs, p::GramMatrix, l, curstate, candidates)
    function rating(dyn, transformation)
        # The major part of the time spent in sosbuildsequence is done
        # in this line so the transformation is cached in `transformation`
        soslf = SetProg.apply_transformation(p, transformation)
        ν = measurefor(μs, dyn)
        # return dot(μ, soslf)
        return dot(getmat(ν), getmat(soslf))
    end
    best_dyn = nothing
    best = -Inf
    for (dyn, transformation) in candidates[curstate]
        cur = rating(dyn, transformation)
        if cur > best
            best = cur
            best_dyn = dyn
        end
    end
    if best_dyn === nothing
        error("$curstate does not have any incoming path of length $l.")
    end
    return best_dyn
end

function apply_map(sw::HybridSystems.DiscreteSwitchingSequence, p, new_vars, d)
    return SetProg.apply_matrix(p, sw.A, new_vars, d)
end

function sosbuilditeration(s::AbstractDiscreteSwitchedSystem, d, seq, μs, p_k::GramMatrix, l, Δt, curstate, iter, candidates)
    best_dyn = best_dynamic(s, d, μs, p_k, l, curstate, candidates)

    curstate = state(s, best_dyn.seq[1], false)
    prepend!(seq, best_dyn)
    x = variables(p_k)
    iter+l, curstate, apply_map(best_dyn, p_k, variables(s, source(s, best_dyn)), d)
end

#function sosbuilditeration(s::AbstractContinuousSwitchedSystem, seq, μs, p_prev, l, Δt, curstate, iter)
#    dyn = best_dynamic(s, μs, p_prev, l, curstate)
#    ub = Δt
#    x = variables(p_prev)
#    p_cur = p_prev(integratorfor(s, (dyn, ub)) * x, x)
#    best_dyn = best_dynamic(s, μs, p_cur, l, curstate)
#    if best_dyn != dyn
#        lb = 0
#        while ub-lb > 1e-5
#            mid = (lb+ub)/2
#            p_cur = p_prev(integratorfor(s, (dyn, mid)) * x, x)
#            best_dyn = best_dynamic(s, μs, p_cur, l, curstate)
#            if best_dyn == dyn
#                lb = mid
#            else
#                ub = mid
#            end
#        end
#    end
#
#    p_cur = p_prev(integratorfor(s, (dyn, ub)) * x, x)
#    if iter > 1 && dyn == seq.seq[iter-1][1] && duration(@view seq.seq[1:(iter-1)]) < 1000# && ub < 1e-3
#        seq.seq[iter-1] = (dyn, seq.seq[iter-1][2]+ub)
#    else
#        push!(seq, (dyn, ub))
#        curstate = state(s, dyn, false)
#        iter += 1
#    end
#    iter, curstate, p_cur
#end

function transformation(s::AbstractSwitchedSystem, path, d)
    xin = variables(s, source(s, path))
    xout = variables(s, target(s, path))
    return SetProg.transformation(monomials(xout, d), dynamicfort(s, path), xin, d)
end

# Extracting trajectory from Lyapunov

"""
    sosbuildsequence(s::AbstractSwitchedSystem, d::Integer;
                     v_0=:Random, p_0=:Random, l::Integer=1,
                     Δt::Float64=1., niter::Integer=42,
                     kws...)

Compute the truncation of length `l` of the high growth rate infinite sequence
produced by the algorithm introduced in [LJP17]. The trajectory ends at mode
`v_0` (or a random one if `v_0` is `:Random`) and is built backward as explained
in [LJP17]. The measures used to guide the construction are the infeasibility
certificates of highest growth rate computed by [`soslyap`](@ref) with
polynomials of degree `2d`. The starting polynomial is either `p_0`, or a random
strictly sum-of-squares polynomial if `p_0` is `:Random` or the primal solution
of [`soslyap`](@ref) certifying the best upper bound for mode `v_0`.

* [LJP17] B. Legat, R. M. Jungers, and P. A. Parrilo.
[Certifying unstability of Switched Systems using Sum of Squares Programming](https://arxiv.org/abs/1710.01814),
arXiv preprint arXiv:1710.01814, **2017**.
"""
function sosbuildsequence(s::AbstractSwitchedSystem, d::Integer;
                          v_0=:Random, p_0=:Random, l::Integer=1,
                          Δt::Float64=1., niter::Integer=42,
                          kws...)
    lyap = getlyap(s, d; kws...)

    if v_0 == :Random
        curstate = rand(states(s))
    else
        if !(v_0 in states(s))
            throw(ArgumentError("Invalid v_0=$v_0"))
        end
        curstate = v_0
    end

    if p_0 == :Primal
        p_k = lyap.primal[curstate].p
    elseif p_0 == :Random
        Z = monomials(variables(s, curstate), d)
        p_k = randsos(Z, monotype=:Gram, r=1)
    else
        # otherwise p_0 is assumed to be an sos polynomial given by the user
        p_k = p_0::GramMatrix
    end

    seq = HybridSystems.switchingsequence(s, niter, curstate)

    candidates = Dict(state => [(dyn, transformation(s, dyn, d)) for dyn in collect(switchings(s, l, state, false))]
                      for state in states(s))

    iter = 1
    while iter <= niter
        iter, curstate, p_k = sosbuilditeration(s, d, seq, lyap.dual, p_k, l, Δt, curstate, iter, candidates)
        # Avoid having it go to zero
        p_k = gram_operate(/, p_k, p_k(variables(p_k) => ones(Int, nvariables(p_k))))
    end
    @assert seq.len == length(seq.seq)
    seq
end


text:
"""
    euclidean(img)
    euclidean(img::BitArray)

Wrapper function for `ImageMorphology.feature_transform` and
`ImageMorphology.distance_transform`. Applies a true Euclidean
distance transform to the array elements and returns an array
with spatial information embedded in the elements.

# Arguments
- img: N-dimensional array to be transformed based on location
    to the nearest background (0) pixel

# Citation
'A Linear Time Algorithm for Computing Exact Euclidean Distance
Transforms of Binary Images in Arbitrary Dimensions' [Maurer et al.,
2003] (DOI: 10.1109/TPAMI.2003.1177156)
"""
euclidean(img) = distance_transform(feature_transform(.!(Bool.(img))))
euclidean(img::BitArray) = distance_transform(feature_transform(.!(img)))


text:
using NamedDims
using NamedDims: names
using SparseArrays
using Test

@testset "get the parent array that was wrapped" begin
    for orig in ([1 2; 3 4], spzeros(2, 2))
        @test parent(NamedDimsArray(orig, (:x, :y))) === orig

        @test unname(NamedDimsArray(orig, (:x, :y))) === orig
        @test unname(orig) === orig
    end
end

@testset "get the named array that was wrapped" begin
    @test names(NamedDimsArray([10 20; 30 40], (:x, :y))) === (:x, :y)
end

@testset "Name-asserting constructor" begin
    orig_full = NamedDimsArray(ones(3, 4, 5), (:a, :b, :c))
    @test names(NamedDimsArray(orig_full, (:a, :b, :c))) == (:a, :b, :c)
    @test names(NamedDimsArray(orig_full, (:a, :b, :_))) == (:a, :b, :c)
    @test_throws DimensionMismatch NamedDimsArray(orig_full, (:a, :b, :wrong))
    @test_throws DimensionMismatch NamedDimsArray(orig_full, (:c, :a, :b))

    orig_partial = NamedDimsArray(ones(3, 4, 5), (:a, :_, :c))
    @test names(NamedDimsArray(orig_partial, (:a, :b, :c))) == (:a, :b, :c)
    @test names(NamedDimsArray(orig_partial, (:a, :_, :c))) == (:a, :_, :c)
end

@testset "getindex" begin
    nda = NamedDimsArray([10 20; 30 40], (:x, :y))

    @test nda[x=1, y=1] == nda[y=1, x=1] == nda[1, 1] == nda[CartesianIndex(1, 1)] == 10
    @test nda[y=end, x=end] == nda[end, end] == 40

    # Unspecified dims become slices
    @test nda[y=1] == nda[y=1, x=:] == nda[:, 1] == [10; 30]

    @test_broken nda[CartesianIndex(1), 1] == nda[1, 1]

    @testset "name preservation" begin
        @test names(nda[y=1]) == (:x, )
        @test names(nda[y=1:1]) == (:x, :y)
    end

    # https://github.com/invenia/NamedDims.jl/issues/8
    @testset "with multiple-wildcards" begin
        nda_mw = NamedDimsArray{(:_, :_, :c)}(ones(10, 20, 30));
        @test nda_mw[c=2] == ones(10, 20)
        @test names(nda_mw[c=2]) == (:_, :_)
    end

    @testset "newaxis" begin
        newaxis = [CartesianIndex()];

        @test names(nda[:, newaxis, :]) == (:x, :_, :y)
        @test size(nda[:, newaxis, :]) == (2, 1, 2)

        @test names(nda[1, newaxis, 1]) == (:_, )
        @test size(nda[1, newaxis, 1]) == (1, )

        @test names(nda[CartesianIndex(1, 1), newaxis]) == (:_, )
        @test size(nda[CartesianIndex(1, 1), newaxis]) == (1, )
        @test nda[CartesianIndex(1, 1), newaxis][1] == nda[1, 1]
    end
end

@testset "views" begin
    nda = NamedDimsArray([10 20; 30 40], (:x, :y))

    @test @view(nda[y=1]) == @view(nda[y=1, x=:]) == @view(nda[:, 1]) == [10; 30]

    @testset "name preservation" begin
        @test names(nda[y=1]) == (:x, )
        @test names(nda[y=1:1]) == (:x, :y)
    end
end

@testset "setindex!" begin
    @testset "by name" begin
        nda = NamedDimsArray([10 20; 30 40], (:x, :y))

        nda[x=1, y=1] = 100
        @test nda == [100 20; 30 40]

        nda[x=1] .= 1000
        @test nda == [1000 1000; 30 40]
    end

    @testset "by position" begin
        nda = NamedDimsArray([10 20; 30 40], (:x, :y))

        nda[1, 1] = 100
        @test nda == [100 20; 30 40]

        nda[1, :] .= 1000
        @test nda == [1000 1000; 30 40]
    end
end

@testset "IndexStyle" begin
    nda = NamedDimsArray([10 20; 30 40], (:x, :y))
    @test IndexStyle(typeof(nda)) == IndexLinear()

    sparse_nda = NamedDimsArray(spzeros(4, 2), (:x, :y))
    @test IndexStyle(typeof(sparse_nda)) == IndexCartesian()
end

@testset "length/size/axes" begin
    nda = NamedDimsArray([10 20; 30 40; 50 60], (:x, :y))

    @test length(nda) == 6

    @test axes(nda) == (1:3, 1:2)
    @test axes(nda, :x) == (1:3) == axes(nda, 1)

    @test size(nda) == (3, 2)
    @test size(nda, :x) == 3 == size(nda, 1)
end

@testset "similar" begin
    nda = NamedDimsArray{(:a, :b, :c, :d)}(ones(10, 20, 30, 40))

    @testset "content" begin
        ndb = similar(nda)
        @test parent(ndb) !== parent(nda)
        @test eltype(ndb) == Float64
        @test size(ndb) == (10, 20 , 30, 40)
        @test names(ndb) == (:a, :b, :c, :d)
    end
    @testset "eltype" begin
        ndb = similar(nda, Char)
        @test parent(ndb) !== parent(nda)
        @test eltype(ndb) == Char
        @test size(ndb) == (10, 20, 30, 40)
        @test names(ndb) == (:a, :b, :c, :d)
    end
    @testset "size" begin
        ndb = similar(nda, Float64, (15, 25, 35, 45))
        @test parent(ndb) !== parent(nda)
        @test eltype(ndb) == Float64
        @test size(ndb) == (15, 25, 35, 45)
        @test names(ndb) == (:a, :b, :c, :d)
    end
    @testset "dim names" begin
        ndb = similar(nda, Float64, (:w, :x, :y, :z))
        @test parent(ndb) !== parent(nda)
        @test eltype(ndb) == Float64
        @test size(ndb) == (10, 20, 30, 40)
        @test names(ndb) == (:w, :x, :y, :z)
    end

    @testset "dimensions" begin
        ndb = similar(nda, Float64, (w=11, x=22))
        @test parent(ndb) !== parent(nda)
        @test eltype(ndb) == Float64
        @test size(ndb) == (11, 22)
        @test names(ndb) == (:w, :x)
    end
end

const cnda = NamedDimsArray([10 20; 30 40], (:x, :y))
@testset "allocations: wrapper" begin
    @test 0 == @allocated parent(cnda)
    @test 0 == @allocated NamedDims.names(cnda)

    @test 0 == @allocated NamedDimsArray(cnda, (:x, :y))
    if VERSION >= v"1.1"
        @test 0 == @allocated NamedDimsArray(cnda, (:x, :_))
    else
        @test_broken 0 == @allocated NamedDimsArray(cnda, (:x, :_))
    end

    # indexing
    @test 0 == @allocated cnda[1,1]
    @test 0 == @allocated cnda[1,1] = 55
    if v"1.1-" <= VERSION <= v"1.2-" # tests pass on 1.1, including 1.1.1
        @test 0 == @allocated cnda[x=1,y=1]
        @test @allocated(cnda[x=1]) == @allocated(cnda[1,:])
        @test 0 == @allocated cnda[x=1,y=1] = 66
    else
        @test_broken 0 == @allocated cnda[x=1,y=1]
        @test_broken @allocated(cnda[x=1]) == @allocated(cnda[1,:])
        @test_broken 0 == @allocated cnda[x=1,y=1] = 66
    end
end


text:
using CurrentPopulationSurvey
using DataDeps
using Tables
using Test

@testset "CurrentPopulationSurvey.jl" begin
    tbl = cpsdata(2019, 1, ["HRINTSTA"])
    @test isdir(@datadep_str "CPS 20191")
    @test Tables.istable(tbl) == true
    @test in(:hrintsta, getfield(tbl, :names))
    tbl_all = cpsdata(2020,1)
    @test in(:pworwgt, getfield(tbl_all, :names))
end


text:
# module MySQL_Octo.SQL

function arrange(::Type{DefaultArrangement}, io::IO, query::QueryBuilders.Queryable)
    write(io, QueryBuilders.colored_repr(query), "\n")
    for row in SQL.iterate(query)
        write(io, string(row, "\n"))
    end
end


text:
abstract type AbstractCovariate end

struct Linear <: AbstractCovariate end
(::Linear)(Y) = [ones(1, size(Y, 2)); Y]

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text: # include("test_actions.jl")


text:
function model_tofts(
    conv::Function = expconv;
    t::AbstractVector,
    params::NamedTuple,
    ca::AbstractVector,
)
    if haskey(params, :ve)
        @extract (kt, ve) params
        kep = kt / ve
    else
        @extract (kt, kep) params
    end
    ct = _model_tofts(conv, t, [kt, kep], ca)
    if haskey(params, :vp)
        ct .+= params.vp * ca
    end
    return ct
end

function model_tofts2(
    t::AbstractVector,
    params::NamedTuple,
    ca::AbstractVector,
    conv::Function = expconv,
)
    if haskey(params, :ve)
        @extract (kt, ve) params
        kep = kt / ve
    else
        @extract (kt, kep) params
    end
    ct = _model_tofts(conv, t, [kt, kep], ca)
    if haskey(params, :vp)
        ct .+= params.vp * ca
    end
    return ct
end

function _model_tofts(
    conv::Function,
    t::AbstractVector,
    params::AbstractVector,
    ca::AbstractVector,
)
    kt, kep = params
    ct = kt * conv(ca, kep, t)
    return ct
end

function _model_extended_tofts(
    conv::Function,
    t::AbstractVector,
    params::AbstractVector,
    ca::AbstractVector,
)
    kt, kep, vp = params
    ct = kt * conv(ca, kep, t)
    ct .+= vp * ca
    return ct
end

function compute_ve_or_kep_if_it_is_missing!(collection)
    for (key, values) in collection
        if !haskey(values, :kep)
            collection[key] = (values..., kep = values.kt ./ values.ve)
        end
        if !haskey(values, :ve)
            collection[key] = (values..., ve = values.kt ./ values.kep)
        end
    end
    return collection
end

function compute_ve_or_kep_if_it_is_missing(namedtuple::NamedTuple)
    if !haskey(namedtuple, :kep)
        return (namedtuple..., kep = namedtuple.kt ./ namedtuple.ve)
    elseif !haskey(namedtuple, :ve)
        return (namedtuple..., ve = namedtuple.kt ./ namedtuple.kep)
    else
        return namedtuple
    end
end

function fit_tofts_nls(
    conv::Function = expconv;
    t::AbstractVector,
    ca::AbstractVector,
    ct::AbstractArray,
    mask = true,
)
    (t, ct, mask, num_timepoints, volume_size) = resolve_fitting_inputs(; t, ca, ct, mask)
    kt, kep = (fill(NaN, volume_size...) for _ = 1:2)
    model(x, p) = _model_tofts(conv, x, p, ca)
    lls_est = fit_tofts_lls(; t, ca, ct, mask).est
    init_kt, init_kep = select(lls_est, (:kt, :kep))
    for idx in eachindex(IndexCartesian(), mask)
        if mask[idx] == false
            continue
        end
        initialvalues = [init_kt[idx], init_kep[idx]]
        (kt[idx], kep[idx]) = curve_fit(model, t, ct[idx, :], initialvalues).param
    end
    est = compute_ve_or_kep_if_it_is_missing((; kt, kep))
    return (; est)
end

function fit_extendedtofts_nls(
    conv::Function = expconv;
    t::AbstractVector,
    ca::AbstractVector,
    ct::AbstractArray,
    mask = true,
)
    (t, ct, mask, num_timepoints, volume_size) = resolve_fitting_inputs(; t, ca, ct, mask)
    kt, kep, vp = (fill(NaN, volume_size...) for _ = 1:3)
    model(x, p) = _model_extended_tofts(conv, x, p, ca)
    lls_est = fit_extendedtofts_lls(; t, ca, ct, mask).est
    init_kt, init_kep, init_vp = select(lls_est, (:kt, :kep, :vp))
    for idx in eachindex(IndexCartesian(), mask)
        if mask[idx] == false
            continue
        end
        initialvalues = [init_kt[idx], init_kep[idx], init_vp[idx]]
        (kt[idx], kep[idx], vp[idx]) = curve_fit(model, t, ct[idx, :], initialvalues).param
    end
    est = compute_ve_or_kep_if_it_is_missing((; kt, kep, vp))
    return (; est)
end

function fit_tofts_lls(;
    t::AbstractVector,
    ca::AbstractVector,
    ct::AbstractArray,
    mask = true,
)
    (t, ct, mask, num_timepoints, volume_size) = resolve_fitting_inputs(; t, ca, ct, mask)
    kt, kep = (fill(NaN, volume_size...) for _ = 1:2)

    M = zeros(num_timepoints, 2)
    M[:, 1] = cumul_integrate(t, ca, TrapezoidalFast())
    for idx in eachindex(IndexCartesian(), mask)
        if mask[idx] == false
            continue
        end
        M[:, 2] = -cumul_integrate(t, ct[idx, :], TrapezoidalFast())
        (kt[idx], kep[idx]) = M \ ct[idx, :]
    end
    est = compute_ve_or_kep_if_it_is_missing((; kt, kep))
    return (; est)
end

function fit_extendedtofts_lls(;
    t::AbstractVector,
    ca::AbstractVector,
    ct::AbstractArray,
    mask = true,
)
    (t, ct, mask, num_timepoints, volume_size) = resolve_fitting_inputs(; t, ca, ct, mask)
    kt, kep, vp = (fill(NaN, volume_size...) for _ = 1:3)

    M = zeros(num_timepoints, 3)
    M[:, 1] = cumul_integrate(t, ca, TrapezoidalFast())
    M[:, 3] = ca
    for idx in eachindex(IndexCartesian(), mask)
        if mask[idx] == false
            continue
        end
        M[:, 2] = -cumul_integrate(t, ct[idx, :], TrapezoidalFast())
        (kt[idx], kep[idx], vp[idx]) = M \ ct[idx, :]
    end
    # Apply correction because fit actually returns: kt + kep * vp
    @. kt = kt - (kep * vp)
    est = compute_ve_or_kep_if_it_is_missing((; kt, kep, vp))
    return (; est)
end


text:
exists(filename::String) = (s = stat(filename); s.inode!=0)
if exists(Pkg.dir("METADATA"))
	pkgs = Pkg.installed();
	map(k -> println(k, ' ', pkgs[k]), sort(collect(keys(pkgs))))
end


text:
using RvSpectraKitLearn

# Likely need to update this to point to directory containing the spectra you wish to analyze
include(joinpath(Pkg.dir("RvSpectraKitLearn"),"path_to_spectra.jl"))
cd(path_to_spectra)

(datafiles, phases) = get_filenames();
(lambda, obs) = read_filelist(datafiles);

lambda_min = 4500
 lambda_max = 4520
 idx_min = findfirst(x->lambda_min<x<lambda_max,lambda)
 idx_max = findlast(x->lambda_min<x<lambda_max,lambda)
 lambda_range = idx_min:idx_max

Y = reshape(-10.0*cos.(phases*2pi),(length(phases),1))
 X = obs[lambda_range,:]'
 @assert size(X,1) == size(Y,1)
 n = size(X,1)
 p = size(X,2)
 q = size(Y,2)

@time (gamma_list, Z_list)  = fit_arp_slr_2d(X,Y,stepsize=1.01,gamma_init=0.001, max_itterations=1000, min_active=floor(Int64,p*q/10))
(L0,L1,L2,rmsdelta) = compute_alr_path_stats(Z_list)
#(chisq,chisq_refit) = compute_alr_path_chisq!(Z_list,X,Y)  # Refitting via least squares makes coefficients worse, since highly underdetermined
print_alr_path(gamma_list, Z_list,X,Y)

#=
# Plot which wavelengths are being used from the last itteration of the algorithm regularization path algorithm
=#
using PyPlot
Z_idx = size(Z_list,1) # floor(Int64,size(Z_list,1)/2)
  Z = reshape(Z_list[Z_idx,:,:],(size(Z_list,2),size(Z_list,3)))
  find(reshape(Z,(p,q)))
  zmax = maximum(Z)
  ymax = maximum(vec(X[1,:]))
  ymin = minimum(vec(X[1,:]))
  yrange = ymax-ymin
  plot(lambda[lambda_range],vec(X[1,:]),"r.")
  plot(lambda[lambda_range],ymin+yrange/2*(1+1/zmax*(reshape(Z,(p,q)))),"b+")
  plot(lambda[lambda_range],ymin+yrange/2*ones(length(lambda_range)),"-")
#=
=#

#=
# Check that the predictions are good
pred = X*reshape(Z,(p,q))
 plot(phases,Y,"r.")
 plot(phases,Y,"r+")
 plot(phases,pred,"b+")
=#

#=
# Random other diagnostic plots
plot(L0,chisq,"r-")
 plot(L0,chisq,"r.")
 plot(L0,chisq_refit,"b-")
 plot(L0,chisq_refit,"b.")
 #plot(L0,chisq_cv,"g-")
 #plot(L0,chisq_cv,"g.")

plt[:loglog](L0,gamma_list,"g.")
 plot(L0,L1./L0,"g-")
 plot(L0,L1./L0,"g.")

plot(L1,chisq,"r-")
 plot(L1,chisq,"r.")
 plot(L1,chisq_refit,"b-")
 plot(L1,chisq_refit,"b.")

plot(L0,chisq,"r-")
 plot(L0,chisq,"r.")

plt[:loglog](gamma_list,L0,"r-")
 plt[:loglog](gamma_list,L1/(p*q),"b-")
 plt[:loglog](gamma_list,sqrt(L2),"g-")

 plt[:loglog](gamma_list,sqrt(L2/(p*q)),"g+")
 plt[:loglog](gamma_list,sqrt(chisq/(n*q)),"r+")
 plt[:loglog](gamma_list,rmsdelta,"g+")

=#

max_stars_repo_name: eford/RVSpectraKitLearn

max_issues_repo_name: eford/RVSpectraKitLearn

max_forks_repo_name: eford/RVSpectraKitLearn


text:
####################################################################################
#  .dolokc          '                        :O                                    #
#  kl    c         :N..        ..  '.        .'         .. .'.           .'. ..    #
#  'xddc,          :X..        .kO..c.       ;N         .Xl..ko         dx..'Xc    #
#  .   .:N.        'X           dd           'N          X'  cO        .M,   0;    #
#  xl'..ck         .X.;,       .Ok'         .cW,        'Nl .xK.        ck'.,X;    #
#    ...    ......   .  ......        ......     ......          ...... .....K.    #
#          .......      .......      .......     ......         .......  ,'''      #
#                               .                   ',.              ..            #
#                              kMx       .    ...   ,Mc         ..   do   ...      #
#                             o.dMl  ,MO.xW' ,;,kN. ,Mc dMc o,cWl,c..XN .Wx;:'     #
#                            c,.'kM; .Mx cM: KX.oM: ,Mc  dWo' .:cOW: KN .::dMo ,,  #
#                           .,   .:;..:;.,:, 'c,':, ':,   k,  .,.,, .;:. ,.';. ,;  #
#                                                       dlc                        #
#                                                                                #
####################################################################################

# StringAnalysis.jl - Library for analyzing text written at 0x0α Research
#                     by Corneliu Cofaru, 2018

"""
A Julia library for working with text, hard-forked from TextAnalysis.jl.
"""
module StringAnalysis
    # Using
    using Unicode
    using SparseArrays
    using LinearAlgebra
    using Statistics
    using Dates
    using DelimitedFiles
    using DataStructures
    using TSVD
    using Languages
    using AutoHashEquals

    # Imports
    import Base: size, show, summary, names
    import Languages: name
    import WordTokenizers

    # Exports
    export AbstractDocument, Document, FileDocument, StringDocument,
           TokenDocument, NGramDocument, GenericDocument,
           AbstractMetadata, DocumentMetadata, metadata,
           text!, text, tokens!, tokens, ngrams!, ngrams,
           ngram_complexity
    export Corpus, DirectoryCorpus, documents, standardize!,
           lexicon, update_lexicon!, create_lexicon,
           lexical_frequency, lexicon_size,
           inverse_index, update_inverse_index!,
           create_inverse_index, index_size
    export DocumentTermMatrix, dtm, dtv, dtv_regex, each_dtv
    export hash_dtv, each_hash_dtv, hash_dtm
    export TextHashFunction, index_hash, cardinality,
           hash_function, hash_function!
    export CooMatrix, coom
    export Stemmer, stem!, stem, stemmer_types
    export tokenize, sentence_tokenize
    export tf!, tf, tf_idf!, tf_idf, bm_25!, bm_25
    export LSAModel, lsa, save_lsa_model, load_lsa_model
    export RPModel, rp, save_rp_model, load_rp_model
    export embed_document, cosine, similarity,
           vocabulary, in_vocabulary, index, get_vector
    export lda
    export frequent_terms, sparse_terms,
           prepare!, prepare,
           strip_patterns, strip_corrupt_utf8, strip_case,
           strip_accents, strip_punctuation, stem_words,
           strip_whitespace, strip_numbers, strip_non_ascii,
           strip_single_chars, strip_html_tags,
           strip_indefinite_articles, strip_definite_articles,
           strip_articles, strip_prepositions,
           strip_pronouns, strip_stopwords,
           strip_sparse_terms, strip_frequent_terms,
           strip_everything, strip_everything_stem

    # Include section
    include("defaults.jl")
    include("hash.jl")
    include("document.jl")
    include("corpus.jl")
    include("metadata.jl")
    include("stemmer.jl")
    include("tokenizer.jl")
    include("ngramizer.jl")
    include("dtm.jl")
    include("coom.jl")
    include("stats.jl")
    include("lsa.jl")
    include("rp.jl")
    include("lda.jl")
    include("preprocessing.jl")
    include("show.jl")

    # Load libstemmer from Snowball_jll
    using Snowball_jll

end


text:
using FactCheck

include("apitest.jl")
include("mathprogtest.jl")
include("testexamples.jl")
include("mathprogtestextra.jl")
include("testconvex.jl")
#include("test_consts.jl") # Something in Julia 0.5 make this test go nuts

FactCheck.exitstatus()


text:
# Random.seed!(1);
if "data" in keys(conf)
    idx_order = randperm(n_exp)
    pert = DataFrame(CSV.File(string(conf["data"],"/pert.csv"); header=false))
    μ_list = convert(Matrix,pert)[idx_order,:]
else
    μ_list = randomLHC(n_exp, ns) ./ n_exp;
    nμ = Int64(conf["n_mu"])
    for i = 1:n_exp
        nonzeros = findall(μ_list[i, :].>0)
        ind_zero = sample(nonzeros, max(0, length(nonzeros)-nμ), replace=false)
        μ_list[i, ind_zero] .= 0
    end
end

function gen_network(m; weight_params=(-1., 1.), sparsity=0., drop_range=(-1e-1, 1e-1))

    # uniform random for W matrix
    w = rand(Uniform(weight_params[1], weight_params[2]), (m, m))

    # Drop small values
    @inbounds for i in eachindex(w)
        w[i] = ifelse(drop_range[1]<=w[i]<=drop_range[2], 0, w[i])
    end

    # Add sparsity
    p = [sparsity, 1 - sparsity]
    w .*= sample([0, 1], weights(p), (m, m), replace=true)

    # Add α vector
    α = abs.(rand(Uniform(weight_params[1], weight_params[2]), (m))) .+ 0.5

    return hcat(α, w)
end

if "network" in keys(conf)
    df = DataFrame(CSV.File(conf["network"]))
    nodes = names(df)
    w = convert(Matrix, df[:,2:end])
    @assert size(w)[1] == size(w)[2]
    @assert size(w)[1] == ns
    if "randomize_network" in keys(conf)
        w_rand = rand(Normal(1, conf["randomize_network"]), (ns, ns))
        w = w .* w_rand
    end
    if "alpha" in keys(conf)
        α = ones(ns) .* conf["alpha"]
    else
        α = ones(ns) .* 0.2
    end
    p_gold = hcat(α, w)
elseif "data" in keys(conf)
    p_gold = hcat(ones(ns),zeros(ns,ns))
else
    p_gold = gen_network(ns, weight_params=(-1.0, 1.0),
                         sparsity=Float64(conf["sparsity"]),
                         drop_range=(Float64(conf["drop_range"]["lb"]), Float64(conf["drop_range"]["ub"])));
end
p = gen_network(ns; weight_params=(0.0, 0.01), sparsity=0);
# show_network(p)

function show_network(p)
    println("p_gold")
    show(stdout, "text/plain", round.(p_gold, digits=2))
    println("\np_learned")
    show(stdout, "text/plain", round.(p, digits=2))
end

function loss_network(p)
     # distalpha = cosine_dist(p_gold[:,1],p[:,1])
     # distw = cosine_dist(p_gold[:,2:end],p[:,2:end])
     @inbounds coralpha = cor(p_gold[:,1],p[:,1])
     @inbounds corw = cor([p_gold[:,2:end]...],[p[:,2:end]...])
     return coralpha, corw
end

function cellbox!(du, u, p, t)
    @inbounds du .= @view(p[:, 1]) .* tanh.(@view(p[:, 2:end]) * u - μ) .- u
end

tspan = (0, tfinal);
ts = 0:tspan[2]/ntotal:tspan[2];
ts = ts[2:end];
u0 = zeros(ns);
prob = ODEProblem(cellbox!, u0, tspan, saveat=ts);

function max_min(ode_data)
    return maximum(ode_data, dims=2) .- minimum(ode_data, dims=2)
end

if "data" in keys(conf)
    expr = DataFrame(CSV.File(string(conf["data"],"/expr.csv"); header=false))
    ode_data_list = zeros(Float64, (n_exp, ns, ntotal));
    ode_data_list[:,:,1] = convert(Matrix,expr)[idx_order,:]
else
    ode_data_list = zeros(Float64, (n_exp, ns, ntotal));
    yscale_list = [];
    for i = 1:n_exp
        global μ = μ_list[i, 1:ns]
        ode_data = Array(solve(prob, Tsit5(), u0=u0, p=p_gold))

        ode_data += randn(size(ode_data)) .* noise
        ode_data_list[i, :, :] = ode_data

        push!(yscale_list, max_min(ode_data))
    end
    yscale = maximum(hcat(yscale_list...), dims=2);
end

function predict_neuralode(u0, p, i_exp=1, batch=ntotal, saveat=true)
    global μ = μ_list[i_exp, 1:ns]
    if saveat
        @inbounds _prob = remake(prob, p=p, tspan=[0, ts[batch]])
        pred = Array(solve(_prob, Tsit5(), saveat=ts[1:batch],
                     sensealg=InterpolatingAdjoint()))
    else # for full trajectory plotting
        @inbounds _prob = remake(prob, p=p, tspan=[0, ts[end]])
        pred = Array(solve(_prob, Tsit5(), saveat=0:ts[end]/nplot:ts[end]))
    end
    return pred
end
predict_neuralode(u0, p, 1);

function loss_neuralode(p, i_exp=1, batch=ntotal)
    pred = predict_neuralode(u0, p, i_exp, batch)
    @inbounds loss = mae(@views(ode_data_list[i_exp, :, 1:batch]), pred)
    return loss
end
@show loss_neuralode(p, 1);


text:
"Propagate Abstract Information around Autumn Program"
module AbstractInterpretation

"Perform abstract interpretation on `aexpr`, propage until converged"
function ainterpret(aexpr, converged)
  while !converged()

  end
end

end


text:
type Flag1
end

type Flag2
end

flag1 = Flag1();
flag2 = Flag2();

immutable ExampleFlags
    flag::Uint8;
    function ExampleFlags(flag1::Flag1)
        return new(convert(Uint8,1));
    end
    function ExampleFlags(flag2::Flag2)
        return new(convert(Uint8,2));
    end
end

type FlaggedType
    data;
    flag::ExampleFlags;
end

function foo(input::FlaggedType)
    if input.flag == ExampleFlags(flag1)
        print("Flag is flag1\n");
    elseif infput.flag == ExampleFlags(flag2)
        print("Flag is flag2\n");
    end
end

max_stars_repo_name: ReidAtcheson/EnumsInJulia

max_issues_repo_name: ReidAtcheson/EnumsInJulia

max_forks_repo_name: ReidAtcheson/EnumsInJulia


text:
using PDESolver  # load main module
using EulerEquationMod  # load physics module

fname = "input_vals1.jl"  # name of input file in each directory

for i=1:9  # loop over directories

  println("\nrunning convergence study in directory m$i")

  cd("./m$i")

  run_solver(fname)  # invoke the solver on the input file in the current
                     # directory

  cd("../")
end


text:
@testset "SortPermutation2AllDifferent: $(fct_type), dimension $(array_dim), $(T)" for fct_type in ["vector of variables", "vector affine function"], array_dim in [2, 3], T in [Int, Float64]
    dim = 3 * array_dim
    mock = MOIU.MockOptimizer(AllDifferentIndexingModel{T}())
    model = COIB.SortPermutation2AllDifferent{T}(mock)

    if T == Int
        @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.Integer)
    end
    @test MOI.supports_constraint(
        model,
        MOI.VectorAffineFunction{T},
        CP.AllDifferent,
    )
    @test MOI.supports_constraint(
        model,
        MOI.VectorAffineFunction{T},
        CP.ElementVariableArray,
    )
    @test MOIB.supports_bridging_constraint(
        model,
        MOI.VectorAffineFunction{T},
        CP.SortPermutation,
    )

    if T == Int
        x, _ = MOI.add_constrained_variables(model, [MOI.Integer() for _ in 1:dim])
    elseif T == Float64
        x = MOI.add_variables(model, dim)
    end

    fct = if fct_type == "vector of variables"
        MOI.VectorOfVariables(x)
    elseif fct_type == "vector affine function"
        MOIU.vectorize(x)
    else
        @assert false
    end

    c = MOI.add_constraint(model, fct, CP.SortPermutation(array_dim))

    for i in 1:dim
        @test MOI.is_valid(model, x[i])
    end
    @test MOI.is_valid(model, c)

    bridge = MOIBC.bridges(model)[MOI.ConstraintIndex{MOI.VectorOfVariables, CP.SortPermutation}(-1)]

    @testset "Bridge properties" begin
        @test MOIBC.concrete_bridge_type(typeof(bridge), MOI.VectorOfVariables, CP.MinimumAmong) == typeof(bridge)
        @test MOIB.added_constrained_variable_types(typeof(bridge)) == Tuple{Type}[]
        @test MOIB.added_constraint_types(typeof(bridge)) == [
            (MOI.VectorAffineFunction{T}, CP.AllDifferent),
            (MOI.ScalarAffineFunction{T}, CP.ElementVariableArray),
            (MOI.ScalarAffineFunction{T}, MOI.LessThan{T}),
        ]

        @test MOI.get(bridge, MOI.NumberOfVariables()) == 0
        @test MOI.get(bridge, MOI.NumberOfConstraints{MOI.VectorAffineFunction{T}, CP.AllDifferent}()) == 1
        @test MOI.get(bridge, MOI.NumberOfConstraints{MOI.ScalarAffineFunction{T}, CP.ElementVariableArray}()) == array_dim
        @test MOI.get(bridge, MOI.NumberOfConstraints{MOI.ScalarAffineFunction{T}, MOI.LessThan{T}}()) == array_dim - 1

        @test MOI.get(bridge, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{T}, CP.AllDifferent}()) == [bridge.con_alldiff]
        @test MOI.get(bridge, MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{T}, CP.ElementVariableArray}()) == bridge.cons_value
        @test MOI.get(bridge, MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{T}, MOI.LessThan{T}}()) == bridge.cons_sort
    end

    @testset "Constraint: all difference" begin
        @test MOI.is_valid(model, bridge.con_alldiff)
        s = MOI.get(model, MOI.ConstraintSet(), bridge.con_alldiff)
        f = MOI.get(model, MOI.ConstraintFunction(), bridge.con_alldiff)

        @test typeof(s) == CP.AllDifferent
        @test length(f.terms) == array_dim
        @test f.constants == zeros(T, array_dim)

        for i in 1:array_dim
            t = f.terms[i]
            @test t.output_index === i
            @test t.scalar_term.coefficient === one(T)
            @test t.scalar_term.variable == x[2 * array_dim + i]
        end
    end

    @testset "Constraints: greater than" begin
        @test length(bridge.cons_sort) == array_dim - 1
        for i in 1:(array_dim - 1)
            @test MOI.is_valid(model, bridge.cons_sort[i])
            s = MOI.get(model, MOI.ConstraintSet(), bridge.cons_sort[i])
            f = MOI.get(model, MOI.ConstraintFunction(), bridge.cons_sort[i])

            @test typeof(s) == MOI.GreaterThan{T}
            @test s.lower == zero(T)
            @test length(f.terms) == 2
            @test f.constant == zero(T)

            t1 = f.terms[1]
            @test t1.coefficient === one(T)
            @test t1.variable == x[i]

            t2 = f.terms[2]
            @test t2.coefficient === -one(T)
            @test t2.variable == x[i + 1]
        end
    end

    @testset "Constraints: array index" begin
        @test length(bridge.cons_value) == array_dim
        for i in 1:array_dim
            @test MOI.is_valid(model, bridge.cons_value[i])
            s = MOI.get(model, MOI.ConstraintSet(), bridge.cons_value[i])
            f = MOI.get(model, MOI.ConstraintFunction(), bridge.cons_value[i])

            @test typeof(s) == CP.ElementVariableArray
            @test length(f.terms) == 2 + array_dim
            @test f.constants == zeros(T, 2 + array_dim)

            t1 = f.terms[1]
            @test t1.output_index === 1
            @test t1.scalar_term.coefficient === one(T)
            @test t1.scalar_term.variable == x[i]

            t2 = f.terms[2]
            @test t2.output_index === 2
            @test t2.scalar_term.coefficient === one(T)
            @test t2.scalar_term.variable == x[2 * array_dim + i]

            t3 = f.terms[3]
            @test t3.output_index === 3
            @test t3.scalar_term.coefficient === one(T)
            @test t3.scalar_term.variable == x[array_dim + 1]

            t4 = f.terms[4]
            @test t4.output_index === 4
            @test t4.scalar_term.coefficient === one(T)
            @test t4.scalar_term.variable == x[array_dim + 2]
        end
    end
end


text:
"""
    NoGradient <: NeutralLandscapeMaker

    NoGradient()

This type is used to generate a random landscape with no gradients
"""
struct NoGradient <: NeutralLandscapeMaker
end

_landscape!(mat, ::NoGradient) = rand!(mat)


text:
abstract type AbstractOutput end

show_output(::AbstractOutput, io, ::MIME"text/plain", t) = nothing

struct TrialClassicalOutputPre <: AbstractOutput end
struct TrialClassicalHistogram <: AbstractOutput end
struct TrialClassicalOutputPost <: AbstractOutput end

function show_output(::TrialClassicalOutputPre, io, ::MIME"text/plain", t::Trial)
    if length(t) > 0
        min = minimum(t)
        max = maximum(t)
        med = median(t)
        avg = mean(t)
        memorystr = string(prettymemory(memory(min)))
        allocsstr = string(iallocs(min))
        minstr = string(prettytime(time(min)), " (", prettypercent(gcratio(min)), " GC)")
        maxstr = string(prettytime(time(max)), " (", prettypercent(gcratio(max)), " GC)")
        medstr = string(prettytime(time(med)), " (", prettypercent(gcratio(med)), " GC)")
        meanstr = string(prettytime(time(avg)), " (", prettypercent(gcratio(avg)), " GC)")
    else
        memorystr = "N/A"
        allocsstr = "N/A"
        minstr = "N/A"
        maxstr = "N/A"
        medstr = "N/A"
        meanstr = "N/A"
    end
    println(io, "BenchmarkExt.Trial:")
    pad = get(io, :pad, "")
    println(io, pad, "  memory estimate:  ", memorystr)
    println(io, pad, "  allocs estimate:  ", allocsstr)
    println(io, pad, "  --------------")
    println(io, pad, "  minimum time:     ", minstr)
    println(io, pad, "  median time:      ", medstr)
    println(io, pad, "  mean time:        ", meanstr)
    println(io, pad, "  maximum time:     ", maxstr)
    println(io, pad, "  --------------")
    println(io, pad, "  samples:          ", length(t))
    print(io,   pad, "  evals/sample:     ", t.params.evals)
end


text:
using Distributions
import Distributions.logpdf

# TODO: iid is currently a bit weird. We'd like it to allow both a specific multiplicity and "on demand" (unspecified but flexible). It's clear how to do either, but not yet if this will allow a single representation.

export iid
struct iid
    size
    dist
end

# TODO: Clean up this hack
iid(n::Int) = dist -> iid(n,dist)

iid(dist) = iid(Nothing, dist)

function Base.iterate(d::iid)
    n = prod(d.size) * length(d.dist)
    return (d.dist, (n, 1))
end

function Base.iterate(d::iid, state)
    (n, count) = state
    count < n && return (d.dist, (n, count + 1))
    return nothing
end

# Distributions.cdf(iid, x) = cdf.(iid.dist, x)
# Distributions.quantile(iid, x) = quantile.(iid.dist, x)

import Base.length
Base.length(d::iid) = prod(d.size)

import Base.eltype
Base.eltype(d::iid) = typeof(d.dist)

rand(d::iid) = rand(d.dist,d.size)

function Distributions.logpdf(d::iid,x)
    s = Float64(0)
    Δs(xj) = logpdf(d.dist, xj)
    @inbounds @simd for j = 1:length(x)
        s += Δs(x[j])
    end
    s
end


text:
module TestConversions
    using Base.Test
    using DataFrames

    df = DataFrame()
    df[:A] = 1:5
    df[:B] = [:A, :B, :C, :D, :E]
    @test isa(array(df), Matrix{Any})
    # @test isa(array(df, Any), Matrix{Any})

    df = DataFrame()
    df[:A] = 1:5
    df[:B] = 1.0:5.0
    @test isa(array(df), Matrix{Real})
    # @test isa(array(df, Any), Matrix{Any})
    # @test isa(array(df, Float64), Matrix{Float64})

    df = DataFrame()
    df[:A] = 1.0:5.0
    df[:B] = 1.0:5.0
    @test isa(array(df), Matrix{Float64})
    # @test isa(matrix(df, Any), Matrix{Any})
    # @test isa(matrix(df, Int), Matrix{Int})
end


text:
##用于实现一些有用的函数

export  add_gaussian_noise,
        hilbert,
        find_fft_peak

#高斯白噪音注入,信噪比越大越好
function add_gaussian_noise(single::AbstractArray{T, 1}, SNR::Float64) where{T<:Number}
    ps = sum(single.^2) / length(single)
    pn_set = sqrt(ps / (10^(SNR / 10)))

    return single + pn_set * randn(Float64, size(single))
end

function add_gaussian_noise(single::AbstractArray{T, 2}, SNR::Float64; dims=1) where{T<:Number}
    buff = similar(single)
    if dims == 1
        @views Threads.@threads for i in 1:size(single, 2)
            buff[:, i] = add_gaussian_noise(single[:, i], SNR)
        end
    elseif dims == 2
        @views Threads.@threads for i in 1:size(single, 1)
            buff[i, :] = add_gaussian_noise(single[i, :], SNR)
        end
    else
        error("dims should be 1 or 2")
    end
    return @view buff[:,:]
end
#希尔伯特变换实现
function hilbert(x::AbstractArray{T, 1}) where{T<:Real}
    buffer = similar(x, Complex)
    buffer = fft(x, 1)
    len = length(x)
    if len %2 ==0
        buffer[2: div(len, 2)] .= 2 .*buffer[2: div(len, 2)]
        buffer[div(len, 2)+1: end] .= 0
    elseif len % 2 == 1
         buffer[2: div(len+1, 2)] .= 2 .*buffer[2: div(len+1, 2)]
         buffer[div(len+1, 2)+1:end] .= 0
    end
    buffer = ifft(buffer, 1)
    buffer1 = similar(x, Complex)
    map!((x, y) -> complex(x, imag(y)), buffer1, x, buffer)
    return view(buffer1,:)
end

function hilbert(x::AbstractArray{T, 2}; dims = 1) where{T <: Real}
    buff =  similar(x, Complex)
    if dims == 1
        @views Threads.@threads for i in 1:size(x, 2)
             buff[:, i] = hilbert(x[:, i])
        end
    elseif  dims == 2
        @views Threads.@threads for i in 1:size(x, 1)
             buff[i, :] = hilbert(x[i, :])
        end
    else
        error("dims should be 1 or 2")
    end

    return @view buff[:,:]
end

##进行傅里叶变换并寻找波峰值,找到前5个波峰
##如果调用上面那种阈值需要调整
function find_fft_peak(x::Vector{T}, fs, Threshold::Float64) where{T}
    xf = abs.(fft(x))
    index = Int[]
    L = length(xf)
    P = (xf ./ L)[1:Int(L/2)]
    @. P[2: end] = 2 * P[2: end]

    f = view(fs*(0:L-1) / L, 1:Int(L/2))

    ##直接for循环遍历
    Threads.@threads for i in 1:length(P)
        P[i] >= Threshold && push!(index, i)
    end

    return view(f, index)
end
function find_fft_peak(x::Vector,fs, n::Int)
    xf = abs.(fft(x))
    L = length(x)
    index = partialsortperm(view(xf, 1:Int(L/2)), 1:n, rev = true)
    f = (fs*(0:L-1)/L)[1:Int(L/2)]

    return view(f, index)
end


text:
include("../Strukture/klubovi.jl")

#Insertion Sort
function insertionSortKlubovi!(A)
    for j in 2:length(A)
        key = A[j];
        i = j-1;
        while(i > 0 && uporedi(key,A[i]))
            A[i+1] = A[i];
            i = i-1;
        end
        A[i+1] = key;
    end
    return A
end

#Merge Sort

function mergeSortKlubovi!(A)
    mergeSortStepKlubovi!(A,1,length(A))
end

function mergeSortStepKlubovi!(A,p,r)
    if p < r
        q = Int(floor((p+r)/2));
        mergeSortStepKlubovi!(A,p,q)
        mergeSortStepKlubovi!(A,q+1,r)
        mergeRezKlubovi!(A,p,q,r)
    end
end

function mergeRezKlubovi!(A,p,q,r)
    L = copy(A[p:q])
    R = copy(A[q+1:r])

    push!(L,Klub("-",typemin(Int64),typemin(Int64),typemin(Int64)))
    push!(R,Klub("-",typemin(Int64),typemin(Int64),typemin(Int64)))
    i = 1;
    j = 1;
    for k = p:r
        if uporedi(L[i],R[j])
            A[k] = L[i];
            i = i + 1;
        else
            A[k] = R[j];
            j = j + 1;
        end
    end
end
mergeSortKlubovi!(klubovi)
ispisKlubova(klubovi)


text:
using DiscreteEvents, Printf, Distributions, Random

const μ = 1/3       # service rate
const λ = 0.9       # arrival rate
count = [0]         # a job counter

# describe the server process
function serve(clk::Clock, id::Int, input::Channel, output::Channel, X::Distribution)
    job = take!(input)
    print(clk, @sprintf("%6.3f: server %d serving customer %d\n", tau(clk), id, job))
    @delay clk X
    print(clk, @sprintf("%6.3f: server %d finished serving %d\n", tau(clk), id, job))
    put!(output, job)
end

# model the arrivals
function arrive(c::Clock, input::Channel, cust::Vector{Int})
    cust[1] += 1
    @printf("%6.3f: customer %d arrived\n", tau(c), cust[1])
    put!(input, cust[1])
end

Random.seed!(123)  # set random number seed
clock = Clock()    # create a clock
input = Channel{Int}(Inf)
output = Channel{Int}(Inf)
for i in 1:3       # start three server processes
    @process serve(clock, i, input, output, Exponential(1/μ))
end
# create a repeating event for 10 arrivals
@event arrive(clock, input, count) every Exponential(1/λ) 10
@run! clock 20     # run the clock


text:
#requires flatten! and unflatten!
import ..flatten!
import ..unflatten!

function swarmoptimize{F}(mlmodel::MLAlgorithm{F}, input::Matrix, answers::Matrix, cost::Function)

  #return this as the stepping function.
  #generate a storage array.
  parameters = length(mlmodel)
  initial_array = Vector{F}(parameters)
  flatten!(initial_array, mlmodel)

  cf = (data_array) -> begin
    unflatten!(mlmodel, data_array)
    nn_results = mlmodel(input)
    cost(answers, nn_results)
  end

  init = (population, a, cf) -> begin
    for idx = 1:size(population,2)
      population[:, idx] = a .* randn(size(population, 1))
    end
    return true
  end

  general_swarmopt(initial_array, cf, initializer = init)
end

#particle swarm optimizer, taken from other code by I. Yonemoto.

function general_swarmopt{F}(a::Array{F,1}, cf::Function; reporter::Function = (c, count, a)->(), popsize = 1000,
	initializer = default_initializer!,
	termcond = (bestscore, iter) -> iter > 100,
	acceptfn = (testscore, bestscore) -> prod(testscore .< bestscore))

	#Swarm optimization function.
	#Prerequisites:
	#		a is a seeding, linear array for the process.  For a default system you should seed with a zeros array.
	#		cf is a cost function that takes an array that looks like the seeding array.
	#		parameters:
	#		reporter is a function that gets called to do reports; c should be the current cost and count should be count.
	#
	#		popsize is how big of a swarm we should use.  Defaults to 5x the vector dimension
	#		maxiter is how many iterations it should take before we quit
	#		velfactor is the initial multiplicand for the velocity.
	#		accelfactor is what factor we should accelerate our swarm by each timestep.
	#		drag is the drag applied to the velocities
	#		sthresh is the score threshold - if we are below this score then we stop.

	#TODO: Some of these parameters will be moved over to the header.
	velfactor = F(0)
	accelfactor = F(0.5)
	drag = F(0.9)

	#determine the length of the array
	l = length(a)

	#set up the population matrix (array of test vectors)
	population = zeros(F, l, popsize)

  initializer(population, a, cf) || return a

	#set up the velocity matrix (array of vector deltas)
	velocities = velfactor * randn(F,l, popsize)

	#set up the vector that stores the best score
	bestscore = cf(population[:,1])
	minidx = 1

	#set up an attractor vector.
	attractor = zeros(F, l)

	#iterate to move the swarm.
  iter = 0
	while (true)
		termcond(bestscore, iter) && break

		iter += 1
		#find scores for the entire population, store them.
		for idx = 1:popsize
			#test this particular population member.
			testscore = cf(population[:,idx])
			#check if testscore is categorically better than the best.
			if acceptfn(testscore, bestscore)
				bestscore = testscore
				minidx = idx
			end
		end

		attractor = population[:,minidx]

		#generate an acceleration matrix and apply it to the velocities
		velocities += (accelfactor) * (attractor * ones(F, 1, popsize) - population)
		velocities[:, minidx] = zeros(F, l)

		#update the swarm
		population += velocities

		#apply drag
		velocities *= drag

		reporter(bestscore, iter, attractor)
	end

	#output the last attractor
	attractor

end #function

max_stars_repo_name: interplanetary-robot/GenML

max_issues_repo_name: interplanetary-robot/GenML

max_forks_repo_name: interplanetary-robot/GenML