text:
defaultViewerHelpText = """
PgUp/PgDn,
Arrow keys : standard navigation
l          : move halfway toward the start
L          : move halfway to the end
Home       : jump to the start
End        : jump to the end
"""

type TwViewerData
    messages::Array
    msglen::Int
    msgwidth::Int
    currentTop::Int
    currentLine::Int
    currentLeft::Int
    showLineInfo::Bool # e.g.1/100 1.0% at top right corner
    bottomText::UTF8String
    trackLine::Bool
    showHelp::Bool
    helpText::UTF8String
    tabWidth::Int
    TwViewerData() = new( UTF8String[], 0, 0, 1, 1, 1, true, "", false, true, defaultViewerHelpText, 4 )
end

# the ways to use it:
# exact dimensions known: h,w,y,x, content to add later
# exact dimensions unknown, but content known and content drives dimensions
function newTwViewer( scr::TwScreen; height::Real=0.5,width::Real=0.8,posy::Any=:staggered,posx::Any=:staggered,
        box=true, showLineInfo=true, showHelp=true, bottomText = "", tabWidth = 4, trackLine = false, title="" )
    obj = TwObj( TwViewerData(), Val{ :Viewer } )
    obj.box = box
    obj.title = title
    obj.borderSizeV= box ? 1 : 0
    obj.borderSizeH= box ? 2 : 0
    obj.data.showLineInfo = showLineInfo
    obj.data.showHelp = showHelp
    obj.data.tabWidth = tabWidth
    obj.data.bottomText = bottomText
    obj.data.trackLine = trackLine
    link_parent_child( scr, obj, height,width,posy,posx )
    obj.data.viewContentHeight = obj.height - ( box? obj.borderSizeV *2 : 1 )
    obj
end

function newTwViewer( scr::TwObj, msgs::Array;
        posy::Any=:staggered,posx::Any=:staggered,
        box=true, showLineInfo=true, bottomText = "", showHelp=true, tabWidth = 4, trackLine = false, title="" )
    newmsgs = map( z->escape_string( replace( z, "\t", repeat( " ", tabWidth ) ) ), msgs )
    obj = TwObj( TwViewerData(), Val{ :Viewer } )
    obj.title = title
    setTwViewerMsgs( obj, newmsgs )
    obj.box = box
    obj.borderSizeV= box ? 1 : 0
    obj.borderSizeH= box ? 2 : 0
    obj.data.showLineInfo = showLineInfo
    obj.data.showHelp = showHelp
    obj.data.tabWidth = tabWidth
    obj.data.bottomText = bottomText
    obj.data.trackLine = trackLine

    h = obj.data.msglen + obj.borderSizeV * 2 + (!box && !isempty( obj.data.bottomText )? 1 : 0 )
    w = max( 25, obj.data.msgwidth + obj.borderSizeH * 2, length(title)+6 )

    link_parent_child( scr, obj, h,w,posy,posx )
    obj
end

function newTwViewer{T<:AbstractString}( scr::TwScreen, msg::T; kwargs... )
    newTwViewer( scr, split(msg,"\n"); kwargs... )
end

function viewContentDimensions( o::TwObj{TwViewerData} )
    vh = o.height
    vstart = 0
    if o.box
        vh -= o.borderSizeV * 2
        vstart = 1
    else
        if !isempty( o.title ) || o.data.showLineInfo
            vh -= 1
            vstart = 1
        end
        if !isempty( o.data.bottomText )
            vh -= 1
        end
    end
    vw = o.width - (o.box ? o.borderSizeH * 2 : 0 )
    (vh, vw, vstart)
end

function draw( o::TwObj{TwViewerData} )
    viewContentHeight, viewContentWidth, viewStartRow = viewContentDimensions( o )

    if o.box
        box( o.window, 0,0 )
    end
    if !isempty( o.title )
        mvwprintw( o.window, 0, (@compat round(Int, ( o.width - length(o.title) )/2 )), "%s", o.title )
    end
    if o.data.showLineInfo
        if o.data.msglen <= o.height - 2 * o.borderSizeV
            msg = "ALL"
        else
            if o.data.trackLine
                msg = @sprintf( "%d/%d %5.1f%%", o.data.currentLine, o.data.msglen,
                    o.data.currentLine / o.data.msglen * 100 )
            else
                msg = @sprintf( "%d/%d %5.1f%%", o.data.currentTop, o.data.msglen,
                    o.data.currentTop / (o.data.msglen - o.height + 2 * o.borderSizeV ) * 100 )
            end
        end
        mvwprintw( o.window, 0, o.width - length(msg)-3, "%s", msg )
    end
    for r in o.data.currentTop:min( o.data.currentTop + viewContentHeight - 1, o.data.msglen )
        s = o.data.messages[r]
        endpos=o.data.currentLeft + o.width - 2 * o.borderSizeH - 1
        if endpos < length(s)
            s = s[o.data.currentLeft:endpos]
        else
            s = s[o.data.currentLeft:end]
        end
        if o.data.trackLine && r == o.data.currentLine
            wattron( o.window, A_BOLD | COLOR_PAIR( o.hasFocus? 15 : 30 ) )
            s *= repeat( " ", max(0,viewContentWidth - length(s) ) )
        end
        mvwprintw( o.window, r - o.data.currentTop + viewStartRow, o.borderSizeH, "%s", s )
        if o.data.trackLine && r == o.data.currentLine
            wattroff( o.window, A_BOLD | COLOR_PAIR( o.hasFocus? 15 : 30 ) )
        end
    end
    if length( o.data.bottomText ) != 0
        mvwprintw( o.window, o.height-1, (@compat round(Int, (o.width - length(o.data.bottomText))/2 )), "%s", o.data.bottomText )
    end
end

function inject( o::TwObj{TwViewerData}, token )
    dorefresh = false
    retcode = :got_it # default behavior is that we know what to do with it
    viewContentHeight, viewContentWidth, viewStartRow = viewContentDimensions( o )

    checkTop = () -> begin
        if o.data.currentTop > o.data.currentLine
            o.data.currentTop = o.data.currentLine
        elseif o.data.currentLine - o.data.currentTop > viewContentHeight - 1
            o.data.currentTop = o.data.currentLine - viewContentHeight + 1
        end
    end
    moveby = n -> begin
        if o.data.trackLine
            o.data.currentLine = max(1, min( o.data.msglen, o.data.currentLine + n) )
            checkTop()
        else
            o.data.currentTop = max( 1, min( o.data.msglen - viewContentHeight, o.data.currentTop + n ) )
        end
        true
    end

    if token == :esc
        retcode = :exit_nothing
    elseif token == :up
        dorefresh = moveby(-1)
    elseif token == :down
        dorefresh = moveby(1)
    elseif token == :left
        if o.data.currentLeft > 1
            o.data.currentLeft -= 1
            dorefresh = true
        else
            beep()
        end
    elseif token == :right
        if o.data.currentLeft + viewContentWidth < o.data.msgwidth
            o.data.currentLeft += 1
            dorefresh = true
        else
            beep()
        end
    elseif token == :pageup
        dorefresh = moveby( -viewContentHeight )
    elseif token == :pagedown
        dorefresh = moveby( viewContentHeight )
    elseif token == :KEY_MOUSE
        (mstate,x,y, bs ) = getmouse()
        if mstate == :scroll_up
            dorefresh = moveby( -(@compat round(Int, viewContentHeight/10 )) )
        elseif mstate == :scroll_down
            dorefresh = moveby( @compat round(Int, viewContentHeight/10 ) )
        elseif mstate == :button1_pressed && o.data.trackLine
            (rely,relx) = screen_to_relative( o.window, y, x )
            if 0<=relx<o.width && 0<=rely<o.height
                o.data.currentLine = o.data.currentTop + rely - o.borderSizeH + 1
                dorefresh = true
            else
                retcode = :pass
            end
        end
    elseif  token == :home
        if o.data.currentTop != 1 || o.data.currentLeft != 1 || o.data.currentLine != 1
            o.data.currentTop = 1
            o.data.currentLeft = 1
            o.data.currentLine = 1
            dorefresh = true
        else
            beep()
        end
    elseif in( token, Any[ Symbol("end") ] )
        if o.data.currentTop + o.height-2 < o.data.msglen
            o.data.currentTop = o.data.msglen - o.height+ 2
            dorefresh = true
        else
            beep()
        end
    elseif (token == :enter || token== Symbol( "return" )) && o.data.trackLine
        if haskey( o.listeners, :select )
            for f in o.listeners[ :select ]
                retcode = f( :select, o )
            end
        end
        dorefresh = true
    elseif token == "L" # move half-way toward the end
        if o.data.trackLine
            target = min( (@compat round(Int, ceil((o.data.currentLine + o.data.msglen)/2))), o.data.msglen )
            if target != o.data.currentLine
                o.data.currentLine = target
                checkTop()
                dorefresh = true
            else
                beep()
            end
        else
            target = min( (@compat round(Int,ceil((o.data.currentTop + o.data.msglen - o.height+2)/2))), o.data.msglen - o.height + 2 )
            if target != o.data.currentTop
                o.data.currentTop = target
                dorefresh = true
            else
                beep()
            end
        end
    elseif token == "l" # move half-way toward the beginning
        if o.data.trackLine
            target = max( (@compat round(Int,floor( o.data.currentLine /2))), 1)
            if target != o.data.currentLine
                o.data.currentLine = target
                checkTop()
                dorefresh = true
            else
                beep()
            end
        else
            target = max( (@compat round(Int, floor( o.data.currentTop /2))), 1)
            if target != o.data.currentTop
                o.data.currentTop = target
                dorefresh = true
            else
                beep()
            end
        end
    else
        retcode = :pass # I don't know what to do with it
    end

    if dorefresh
        refresh(o)
    end

    return retcode
end

function helptext( o::TwObj{TwViewerData} )
    if o.data.showHelp
        o.data.helpText
    else
        ""
    end
end

function setTwViewerMsgs( o::TwObj{TwViewerData}, msgs::Array )
    o.data.messages = msgs
    o.data.msglen = length(msgs)
    o.data.msgwidth = maximum( map( x->length(x), msgs ) )
end


text:
#######################################################################
# BoxIterator
#######################################################################

"""
    BoxIterator(seed::SVector{N,Int}; maxiterations = missing)

Cartesian iterator `iter` over `SVector{N,Int}`s (`cell`s) that starts at `seed` and
grows outwards in the form of a box of increasing sides (not necesarily equal) until
it encompasses a certain N-dimensional region. To signal that a cell is in the desired
region the user calls `acceptcell!(iter, cell)`.
"""
struct BoxIterator{N,MI<:Union{Int,Missing}}
    seed::SVector{N,Int}
    maxiter::MI
    dimdir::MVector{2,Int}
    nmoves::MVector{N,Bool}
    pmoves::MVector{N,Bool}
    npos::MVector{N,Int}
    ppos::MVector{N,Int}
end

Base.IteratorSize(::BoxIterator) = Base.SizeUnknown()

Base.IteratorEltype(::BoxIterator) = Base.HasEltype()

Base.eltype(::BoxIterator{N}) where {N} = SVector{N,Int}

Base.CartesianIndices(b::BoxIterator) =
    CartesianIndices(UnitRange.(Tuple(b.npos), Tuple(b.ppos)))

function BoxIterator(seed::SVector{N}; maxiterations = missing) where {N}
    BoxIterator(seed, maxiterations, MVector(1, 2),
        ones(MVector{N,Bool}), ones(MVector{N,Bool}), MVector{N,Int}(seed), MVector{N,Int}(seed))
end

function iteratorreset!(b::BoxIterator{N}) where {N}
    b.dimdir[1] = 1
    b.dimdir[2] = 2
    b.nmoves .= ones(MVector{N,Bool})
    b.pmoves .= ones(MVector{N,Bool})
    b.npos   .= MVector(b.seed)
    b.ppos   .= MVector(b.seed)
    return nothing
end

struct BoxIteratorState{N}
    range::CartesianIndices{N, NTuple{N,UnitRange{Int}}}
    rangestate::CartesianIndex{N}
    iteration::Int
end

Base.iterate(b::BoxIterator{0}) = (SVector{0,Int}(), nothing)
Base.iterate(b::BoxIterator{0}, state) = nothing

function Base.iterate(b::BoxIterator{N}) where {N}
    range = CartesianIndices(ntuple(i -> b.seed[i]:b.seed[i], Val(N)))
    itrange = iterate(range)
    if itrange === nothing
        return nothing
    else
        (cell, rangestate) = itrange
        return (SVector(Tuple(cell)), BoxIteratorState(range, rangestate, 1))
    end
end

function Base.iterate(b::BoxIterator{N}, s::BoxIteratorState{N}) where {N}
    itrange = iterate(s.range, s.rangestate)
    facedone = itrange === nothing
    if facedone
        alldone = !any(b.pmoves) && !any(b.nmoves) || isless(b.maxiter, s.iteration)
        if alldone  # Last shells in all directions were empty, trim from boundingboxcorners
            b.npos .+= 1
            b.ppos .-= 1
            return nothing
        else
            newrange = nextface!(b)
            # newrange === nothing && return nothing
            itrange = iterate(newrange)
            # itrange === nothing && return nothing
            (cell, rangestate) = itrange
            return (SVector(Tuple(cell)), BoxIteratorState(newrange, rangestate, s.iteration + 1))
        end
    else
        (cell, rangestate) = itrange
        return (SVector(Tuple(cell)), BoxIteratorState(s.range, rangestate, s.iteration + 1))
    end
end

function nextface!(b::BoxIterator{N}) where {N}
    @inbounds for i in 1:2N
        nextdimdir!(b)
        newdim, newdir = Tuple(b.dimdir)
        if newdir == 1
            if b.nmoves[newdim]
                b.npos[newdim] -= 1
                b.nmoves[newdim] = false
                return newrangeneg(b, newdim)
            end
        else
            if b.pmoves[newdim]
                b.ppos[newdim] += 1
                b.pmoves[newdim] = false
                return newrangepos(b, newdim)
            end
        end
    end
    return nothing
end

function nextdimdir!(b::BoxIterator{N}) where {N}
    dim, dir = Tuple(b.dimdir)
    if dim < N
        dim += 1
    else
        dim = 1
        dir = ifelse(dir == 1, 2, 1)
    end
    b.dimdir[1] = dim
    b.dimdir[2] = dir
    return nothing
end

@inline function newrangeneg(b::BoxIterator{N}, dim) where {N}
    return CartesianIndices(ntuple(
        i -> b.npos[i]:(i == dim ? b.npos[i] : b.ppos[i]),
        Val(N)))
end

@inline function newrangepos(b::BoxIterator{N}, dim) where {N}
    return CartesianIndices(ntuple(
        i -> (i == dim ? b.ppos[i] : b.npos[i]):b.ppos[i],
        Val(N)))
end

function acceptcell!(b::BoxIterator{N}, cell) where {N}
    dim, dir = Tuple(b.dimdir)
    if dir == 1
        @inbounds for i in 1:N
            (cell[i] == b.ppos[i]) && (b.pmoves[i] = true)
            (i == dim || cell[i] == b.npos[i]) && (b.nmoves[i] = true)
        end
    else
        @inbounds for i in 1:N
            (i == dim || cell[i] == b.ppos[i]) && (b.pmoves[i] = true)
            (cell[i] == b.npos[i]) && (b.nmoves[i] = true)
        end
    end
    return nothing
end

# Fallback for non-BoxIterators
acceptcell!(b, cell) = nothing

#######################################################################
# CoSort
#######################################################################

struct CoSortTup{T,Tv}
    x::T
    y::Tv
end

mutable struct CoSort{T,Tv,S<:AbstractVector{T},C<:AbstractVector{Tv}} <: AbstractVector{CoSortTup{T,Tv}}
    sortvector::S
    covector::C
    offset::Int
    function CoSort{T,Tv,S,C}(sortvector, covector, offset) where {T,Tv,S,C}
        length(covector) >= length(sortvector) ? new(sortvector, covector, offset) :
            throw(DimensionMismatch("Coarray length should exceed sorting array"))
    end
end

CoSort(sortvector::S, covector::C) where {T,Tv,S<:AbstractVector{T},C<:AbstractVector{Tv}} =
    CoSort{T,Tv,S,C}(sortvector, covector, 0)

Base.size(c::CoSort) = (size(c.sortvector, 1) - c.offset,)

Base.getindex(c::CoSort, i) =
    CoSortTup(getindex(c.sortvector, i + c.offset), getindex(c.covector, i + c.offset))

Base.setindex!(c::CoSort, t::CoSortTup, i) =
    (setindex!(c.sortvector, t.x, i + c.offset); setindex!(c.covector, t.y, i + c.offset); c)

Base.isless(a::CoSortTup, b::CoSortTup) = isless(a.x, b.x)

Base.Sort.defalg(v::C) where {T<:Union{Number, Missing}, C<:CoSort{T}} = Base.DEFAULT_UNSTABLE

isgrowing(c::CoSort) = isgrowing(c.sortvector, c.offset + 1)

#######################################################################
# SparseMatrixBuilder
#######################################################################

mutable struct SparseMatrixBuilder{T,S<:AbstractSparseMatrix{T}}
    matrix::S
    colcounter::Int
    rowvalcounter::Int
    cosorter::CoSort{Int,T,Vector{Int},Vector{T}}
end

struct UnfinalizedSparseMatrixCSC{Tv,Ti} <: AbstractSparseMatrix{Tv,Ti}
    m::Int
    n::Int
    colptr::Vector{Ti}
    rowval::Vector{Ti}
    nzval::Vector{Tv}
end

function SparseMatrixBuilder{Tv}(m, n, nnzguess = missing) where {Tv}
    colptr = [1]
    rowval = Int[]
    nzval = Tv[]
    matrix = UnfinalizedSparseMatrixCSC(m, n, colptr, rowval, nzval)
    builder = SparseMatrixBuilder(matrix, 1, 0, CoSort(rowval, nzval))
    nnzguess === missing || sizehint!(builder, nnzguess)
    return builder
end

# Unspecified size constructor
SparseMatrixBuilder{Tv}() where Tv = SparseMatrixBuilder{Tv}(0, 0)

function SparseMatrixBuilder(s::SparseMatrixCSC{Tv,Int}) where {Tv}
    colptr = getcolptr(s)
    rowval = rowvals(s)
    nzval = nonzeros(s)
    nnzguess = length(nzval)
    resize!(rowval, 0)
    resize!(nzval, 0)
    resize!(colptr, 1)
    colptr[1] = 1
    builder = SparseMatrixBuilder(s, 1, 0, CoSort(rowval, nzval))
    sizehint!(builder, nnzguess)
    return builder
end

SparseArrays.getcolptr(s::UnfinalizedSparseMatrixCSC) = s.colptr
SparseArrays.nonzeros(s::UnfinalizedSparseMatrixCSC) = s.nzval
SparseArrays.rowvals(s::UnfinalizedSparseMatrixCSC) = s.rowval
Base.size(s::UnfinalizedSparseMatrixCSC) = (s.m, s.n)
Base.size(s::UnfinalizedSparseMatrixCSC, k) = size(s)[k]

function Base.sizehint!(s::SparseMatrixBuilder, n)
    sizehint!(getcolptr(s.matrix), n + 1)
    sizehint!(nonzeros(s.matrix), n)
    sizehint!(rowvals(s.matrix), n)
    return s
end

function pushtocolumn!(s::SparseMatrixBuilder, row::Int, x, skipdupcheck::Bool = true)
    nrows = size(s.matrix, 1)
    nrows == 0 || 1 <= row <= size(s.matrix, 1) || throw(ArgumentError("tried adding a row $row out of bounds ($(size(s.matrix, 1)))"))
    if skipdupcheck || !isintail(row, rowvals(s.matrix), getcolptr(s.matrix)[s.colcounter])
        push!(rowvals(s.matrix), row)
        push!(nonzeros(s.matrix), x)
        s.rowvalcounter += 1
    end
    return s
end

function isintail(element, container, start::Int)
    for i in start:length(container)
        container[i] == element && return true
    end
    return false
end

function finalizecolumn!(s::SparseMatrixBuilder, sortcol::Bool = true)
    size(s.matrix, 2) > 0 && s.colcounter > size(s.matrix, 2) && throw(DimensionMismatch("Pushed too many columns to matrix"))
    if sortcol
        s.cosorter.offset = getcolptr(s.matrix)[s.colcounter] - 1
        sort!(s.cosorter)
        isgrowing(s.cosorter) || throw(error("Internal error: repeated rows"))
    end
    s.colcounter += 1
    push!(getcolptr(s.matrix), s.rowvalcounter + 1)
    return nothing
end

function finalizecolumn!(s::SparseMatrixBuilder, ncols::Int)
    for _ in 1:ncols
        finalizecolumn!(s)
    end
    return nothing
end

function SparseArrays.sparse(s::SparseMatrixBuilder{<:Any,<:SparseMatrixCSC})
    completecolptr!(getcolptr(s.matrix), size(s.matrix, 2), s.rowvalcounter)
    return s.matrix
end

function completecolptr!(colptr, cols, lastrowptr)
    colcounter = length(colptr)
    if colcounter < cols + 1
        resize!(colptr, cols + 1)
        for col in (colcounter + 1):(cols + 1)
            colptr[col] = lastrowptr + 1
        end
    end
    return colptr
end

function SparseArrays.sparse(s::SparseMatrixBuilder{<:Any,<:UnfinalizedSparseMatrixCSC})
    m, n = size(s.matrix)
    rowval = rowvals(s.matrix)
    colptr = getcolptr(s.matrix)
    nzval = nonzeros(s.matrix)
    if m != 0 && n != 0
        completecolptr!(colptr, n, s.rowvalcounter)
    else # determine size of matrix after the fact
        m, n = isempty(rowval) ? 0 : maximum(rowval), s.colcounter - 1
    end
    return SparseMatrixCSC(m, n, colptr, rowval, nzval)
end

#######################################################################
# SparseMatrixReader
#######################################################################

struct SparseMatrixReader{Tv,Ti}
    matrix::SparseMatrixCSC{Tv,Ti}
end

Base.IteratorSize(::SparseMatrixReader) = Base.HasLength()

Base.IteratorEltype(::SparseMatrixReader) = Base.HasEltype()

Base.eltype(::SparseMatrixReader{Tv,Ti}) where {Tv,Ti} = Tuple{Ti,Ti,Tv}

Base.length(s::SparseMatrixReader) = nnz(s.matrix)

function Base.iterate(s::SparseMatrixReader, state = (1, 1))
    (ptr, col) = state
    ptr > length(s) && return nothing
    @inbounds while ptr > s.matrix.colptr[col + 1] - 1
         col += 1
    end
    return (s.matrix.rowval[ptr], col, s.matrix.nzval[ptr], ptr), (ptr + 1, col)
end

enumerate_sparse(s::SparseMatrixCSC) = SparseMatrixReader(s)

#######################################################################
# SiteSublats
#######################################################################

struct SiteSublats{G,U}
    g::G
    u::U
end

sitesublats(u) = SiteSublats(((i, s) for s in sublats(u) for i in siterange(u, s)), u)

Base.iterate(s::SiteSublats, x...) = iterate(s.g, x...)

Base.IteratorSize(::SiteSublats) = Base.HasLength()

Base.IteratorEltype(::SiteSublats) = Base.HasEltype()

Base.eltype(::SiteSublats) = Tuple{Int,Int}

Base.length(s::SiteSublats) = nsites(s.u)

#######################################################################
# Runs
#######################################################################
# iteration yields ranges of subsequent xs elements such that istogether on them gives true
struct Runs{T,F}
    xs::Vector{T}
    istogether::F
end

equalruns(xs) = Runs(xs, ==)
approxruns(xs::Vector{T}) where {T<:Number} = Runs(xs, (x, y) -> isapprox(x, y; atol = sqrt(eps(real(T)))))

function last_in_run(xs::Vector{T}, i, istogether) where {T}
    xi = xs[i]
    for j in i:length(xs)
        (j == length(xs) || !istogether(xs[j+1], xi)) && return j
    end
    return i
end

# Base.iterate(s::Runs) = iterate(s, 1)
function Base.iterate(s::Runs, j = 1)
    j > length(s.xs) && return nothing
    lst = last_in_run(s.xs, j, s.istogether)
    return j:lst, lst + 1
end

Base.IteratorSize(::Runs) = Base.SizeUnknown()

Base.IteratorEltype(::Runs) = Base.HasEltype()

Base.eltype(::Runs) = UnitRange{Int}

#######################################################################
# MarchingNeighbors
#######################################################################
struct MarchingNeighbors{D,G}
    n::CartesianIndex{D}
    gen::G
    len::Int
end

# c is a CartesianRange over vertices, and n is a CartesianIndex from which to find neighbors
function marchingneighbors(c, n)
    rp = max(n, first(c)):min(n + oneunit(n), last(c))
    rm = max(n - oneunit(n), first(c)):min(n, last(c))
    gen = Iterators.filter(!isequal(n), Iterators.flatten((rp, rm)))
    len = length(rp) + length(rm) - 2
    return MarchingNeighbors(n, gen, len)
end

function marchingneighbors_forward(c, n)
    rp = max(n, first(c)):min(n + oneunit(n), last(c))
    gen = Iterators.filter(!isequal(n), rp)
    len = length(rp) - 2
    return MarchingNeighbors(n, gen, len)
end

Base.iterate(m::MarchingNeighbors, s...) = iterate(m.gen, s...)

Base.IteratorSize(::MarchingNeighbors) = Base.HasLength()

Base.IteratorEltype(::MarchingNeighbors) = Base.HasEltype()

Base.eltype(::MarchingNeighbors{D}) where {D} = CartesianIndex{D}

Base.length(m::MarchingNeighbors) = m.len

#######################################################################
# MarchingSimplices
#######################################################################
struct MarchingSimplices{D,D´}
    cinds::CartesianIndices{D´,NTuple{D´,Base.OneTo{Int}}}
    simps::Vector{NTuple{D´,CartesianIndex{D}}}  # simplices verts in unit cuboid
    unitperms::Vector{NTuple{D,Int}}             # unit vector order of simplices
end

# c is a CartesianIndices over all vertices
function marchingsimplices(c::CartesianIndices{D}) where {D}
    unitperms = marchingsimplices_unitperms(Val(D))
    edgeperms = map(edge -> unitvector.(CartesianIndex{D}, edge), unitperms)
    simps = marchingsimplices_verts.(edgeperms)
    srange = eachindex(simps)
    refverts = ntuple(i -> Base.OneTo(c.indices[i][1:end-1]), Val(D))
    cinds = CartesianIndices((refverts..., srange))
    return MarchingSimplices(cinds, simps, unitperms)
end

# indices of unit vectors in a unit cuboid, taken in any order
marchingsimplices_unitperms(::Val{D}) where {D} = permutations(ntuple(identity, Val(D)))

# simplex vertices as the sum of unit vectors in any permutation
marchingsimplices_verts(edges::NTuple{D}) where {D} =
    (zero(CartesianIndex{D}), edges...) |> cumsum |> orientsimp

function orientsimp(sverts)
    volume = det(hcat(SVector.(Tuple.(Base.tail(sverts)))...))
    sverts´ = ifelse(volume >= 0, sverts, switchlast(sverts))
    return sverts´
end

switchlast(s::NTuple{N,T}) where {N,T} = ntuple(i -> i < N - 1 ? s[i] : s[2N - i - 1] , Val(N))

function Base.iterate(m::MarchingSimplices, s...)
    it = iterate(m.cinds, s...)
    it === nothing && return nothing
    c, s´ = it
    t = Tuple(c)
    nsimp = last(t)
    refvert = CartesianIndex(Base.front(t))
    verts = Ref(refvert) .+ m.simps[nsimp]
    return verts, s´
end

Base.IteratorSize(::MarchingSimplices{D}) where {D} = Base.HasShape{D}()

Base.axes(m::MarchingSimplices, i...) = axes(m.cinds, i...)

Base.size(m::MarchingSimplices, i...) = size(m.cinds, i...)

Base.CartesianIndices(m::MarchingSimplices) = m.cinds

Base.IteratorEltype(::MarchingSimplices) = Base.HasEltype()

Base.eltype(::MarchingSimplices{D,D´}) where {D,D´} = NTuple{D´,CartesianIndex{D}}

Base.length(m::MarchingSimplices) = prod(length.(axes(m)))

ext: jl

max_stars_repo_name: pablosanjose/Quantica

max_issues_repo_name: pablosanjose/Quantica

max_forks_repo_name: pablosanjose/Quantica


text:
using HttpServer

function start_auth_server(client::OutlookClient)
    println("Starting authentication server.")

    http = HttpHandler() do req::Request, res::Response
        finish_authentication(req.resource, client)
        stop_auth_server(server)
        Response(200) # Unrelevant at this point if a respose is actually sent, but without this we get "MethodError"
    end

    server = Server(http)
    run(server, get_auth_server_port())
end

stop_auth_server(server) = close(server)


text:
module TestMLJUnit

using Test

import MLJTuning

using TreeParzen

struct DummyModel
    init::Bool # dummy field to initialise because we can't do DummyModel()
end
DummyModel(;kwargs...) = DummyModel(true) # we literally don't care,
                                          # we just need a constructor

const MLJTreeParzenTuning = MLJTreeParzen.MLJTreeParzenTuning

function setup(;n_simultaneous=1,
               n_startup=20,
               space=Dict{Symbol},
               suggest=Dict{Symbol}[])

    tuning = MLJTreeParzenTuning(;max_simultaneous_draws=n_simultaneous,
                                 random_trials=n_startup)
    state = (space=MLJTreeParzen.MLJTreeParzenSpace(space, suggest),
             trialhist=TreeParzen.Trials.Trial[])

    return tuning, state
end

function complete_trials(array_of_stoof)

    trials = last.(array_of_stoof)

    for trial in trials
        TreeParzen.tell!(trial, rand(Float64))
    end

    return [(DummyModel(), (trial_object=trial,)) for trial in trials]

end

@testset "MLJ: default_n" begin

    obj = MLJTreeParzen.MLJTreeParzenTuning()
    @test MLJTuning.default_n(obj, nothing) == 50

end

@testset "MLJ: setup" begin

    tuning = MLJTreeParzen.MLJTreeParzenTuning()
    model = nothing # supposed to be a model but nobody is checking
                    # and we dont need it so lets prove it
    verbosity = nothing

    space_dict = Dict(:x => HP.Uniform(:x, 0.0, 1.0))
    suggestions = Dict(:x => 0.222)
    suggestions_list = [Dict(:x => 0.222), Dict(:x => 0.444)]
    preconstructed_space = MLJTreeParzen.MLJTreeParzenSpace(space_dict)

    preconstruted_suggest =
        MLJTreeParzen.MLJTreeParzenSpace(space_dict, suggestions)

    preconstruted_suggestions =
        MLJTreeParzen.MLJTreeParzenSpace(space_dict, suggestions_list)

    # test with just a dict space:
    state = MLJTuning.setup(tuning, model, space_dict, 123, verbosity)
    @test state.trialhist isa Vector{<:TreeParzen.Trials.Trial}
    @test state.space isa MLJTreeParzen.MLJTreeParzenSpace

    # test with preconstructed space and that we get same thing:
    state = MLJTuning.setup(tuning, model, preconstructed_space, 123, verbosity)
    @test state.space isa MLJTreeParzen.MLJTreeParzenSpace
    @test state.space == preconstructed_space

    # test with preconstructed space with a sugggestion and we get the
    # same thing:
    state = MLJTuning.setup(tuning, model, preconstruted_suggest, 123, verbosity)
    @test state.space isa MLJTreeParzen.MLJTreeParzenSpace
    @test state.space == preconstruted_suggest

    # test with preconstructed space with multiple suggestions and we
    # get the same thing:
    state = MLJTuning.setup(tuning, model, preconstruted_suggestions, 123, verbosity)
    @test state.space isa MLJTreeParzen.MLJTreeParzenSpace
    @test state.space == preconstruted_suggestions

end

@testset "MLJ: models" begin

    # scenarios to cover (or properties to check):
    # 1) with suggestions < num startup (return num startup followed by num simultaenous draw)
    # 2) with suggestions > num startup (throw exception)
    # 3) with suggestions == num startup (return num startup followed by num simultaenous draw)
    # 4) no suggestions (return num startup followed by num simultaenous draw)
    # 5) round 1 always returns num startup and all subsequent rounds returns num simultaneous draw

    space = Dict(:x => HP.Uniform(:x, -5., 5.))
    testmodel = DummyModel()
    fakehist = complete_trials([(nothing, TreeParzen.ask(space)) for i in 1:100]) # do too many, we can cut it down

    @testset "no suggestions" begin

        suggestions = Dict{Symbol}[]
        tuning, state = setup(;n_startup=3,space=space, suggest=suggestions)

        output, _ = MLJTuning.models(tuning, testmodel, nothing, state, 0, 0)

        @test output isa Vector{<:Tuple{Any, TreeParzen.Trials.Trial}}
        @test length(output) == 3

    end

    @testset "less suggestions" begin

        suggestions = [Dict(:x => 1), Dict(:x => 2)]
        tuning, state = setup(;n_startup=3,space=space, suggest=suggestions)

        output, _ = MLJTuning.models(tuning, testmodel, nothing, state, 0, 0)

        @test output isa Vector{<:Tuple{Any, TreeParzen.Trials.Trial}}
        @test length(output) == 3
        @test getindex.(getproperty.(last.(output), :hyperparams), :x)[1:2] == [1, 2]

    end

        @testset "exact suggestions" begin

        suggestions = [Dict(:x => 1), Dict(:x => 2), Dict(:x => 3)]
        tuning, state = setup(;n_startup=3,space=space, suggest=suggestions)

        output, _ = MLJTuning.models(tuning, testmodel, nothing, state, 0, 0)

        @test output isa Vector{<:Tuple{Any, TreeParzen.Trials.Trial}}
        @test length(output) == 3
        @test getindex.(getproperty.(last.(output), :hyperparams), :x) == [1, 2, 3]

    end

        @testset "too many suggestions" begin

        suggestions = [Dict(:x => 1), Dict(:x => 2), Dict(:x => 3), Dict(:x => 4)]
        tuning, state = setup(;n_startup=3,space=space, suggest=suggestions)

        @test_throws ArgumentError MLJTuning.models(tuning, testmodel, nothing, state, 0, 0)

    end

    @testset "num draws" begin

        for sim_value in 1:10

            suggestions = Dict{Symbol}[]
            tuning, state = setup(;n_simultaneous=sim_value,space=space, suggest=suggestions)

            output, _ = MLJTuning.models(tuning, testmodel, nothing, state, 0, 0)

            @test output isa Vector{<:Tuple{Any, TreeParzen.Trials.Trial}}
            @test length(output) == 20

            fakeMLJhist = map(fakehist) do entry
                (model       = entry[1],
                 measure     = [MLJTuning.rmse],
                 measurement = rand(),
                 metadata    = entry[2].trial_object)
            end

            for histlen in 21:30

                output, state = MLJTuning.models(tuning,
                                           testmodel,
                                           fakeMLJhist[1:histlen],
                                           state,
                                           0,
                                           0)
                @test length(output) == sim_value
                @test length(state.trialhist) == histlen

            end

        end

    end

end

end # module
true


text:
function process_value(::Val{:control}, key, val, global_p, current_p)
    control_block_global_parameters = (
        r"n[x,y,z]",
        "nparticles",
        "nsteps",
        "t_end",
        r"[x,y,z]_(min|max|end)",
    )

    process_value(key,
        val,
        global_p,
        current_p,
        control_block_global_parameters)
end


text:
using DocumentDistances: pdf2text, withtxtext, walkfilesmatching

const ExampleFileDir = joinpath(dirname(@__FILE__()), "data", "pdffile")
const ExamplePdfFile = joinpath(ExampleFileDir, "testfile.pdf")
const ExampleTxtFile = joinpath(ExampleFileDir, "testfile.txt")

@testset "withtxtext" begin
    @test withtxtext("test.pdf") == "test.txt"
    @test withtxtext("/my/path/test_123.pdf") == "/my/path/test_123.txt"
end # @testset "withtxtext"

@testset "pdf2text" begin
    isfile(ExampleTxtFile) && rm(ExampleTxtFile)
    txtfilename = pdf2text(ExamplePdfFile)
    @test txtfilename == ExampleTxtFile
    @test isfile(ExampleTxtFile)
    isfile(ExampleTxtFile) && rm(ExampleTxtFile) # Clean up again after us
end # @testset "pdf2text"

@testset "walkfilesmatching" begin
    pdfs = String[]
    walkfilesmatching(fp -> push!(pdfs, fp), ExampleFileDir, r"\.pdf$")
    @test length(pdfs) == 1
end # @testset "pdf2text"


text:
using Flux
using MacroTools: @forward

abstract type AbstractOptimiser end

const EPS = 1e-8

# TODO: should use weak refs

"""
    Descent(η = 0.1)

Classic gradient descent optimiser with learning rate `η`.
For each parameter `p` and its gradient `δp`, this runs `p -= η*δp`

# Parameters
- Learning rate (`η`): Amount by which gradients are discounted before updating
                       the weights.

# Examples
```julia
opt = Descent()

opt = Descent(0.3)

ps = Flux.params(model)

gs = gradient(ps) do
    loss(x, y)
end

Flux.Optimise.update!(opt, ps, gs)
```
"""
mutable struct Descent <: AbstractOptimiser
  eta::Float64
end

Descent() = Descent(0.1)

function apply!(o::Descent, x, Δ)
  Δ .*= o.eta
end

"""
    Momentum(η = 0.01, ρ = 0.9)

Gradient descent optimizer with learning rate `η` and momentum `ρ`.

# Parameters
- Learning rate (`η`): Amount by which gradients are discounted before updating
                       the weights.
- Momentum (`ρ`): Controls the acceleration of gradient descent in the
                  prominent direction, in effect damping oscillations.

# Examples
```julia
opt = Momentum()

opt = Momentum(0.01, 0.99)
```
"""
mutable struct Momentum <: AbstractOptimiser
  eta::Float64
  rho::Float64
  velocity::IdDict
end

Momentum(η = 0.01, ρ = 0.9) = Momentum(η, ρ, IdDict())

function apply!(o::Momentum, x, Δ)
  η, ρ = o.eta, o.rho
  v = get!(() -> zero(x), o.velocity, x)::typeof(x)
  @. v = ρ * v - η * Δ
  @. Δ = -v
end

"""
    Nesterov(η = 0.001, ρ = 0.9)

Gradient descent optimizer with learning rate `η` and Nesterov momentum `ρ`.

# Parameters
- Learning rate (`η`): Amount by which gradients are discounted before updating
                       the weights.
- Nesterov momentum (`ρ`): Controls the acceleration of gradient descent in the
                           prominent direction, in effect damping oscillations.

# Examples
```julia
opt = Nesterov()

opt = Nesterov(0.003, 0.95)
```
"""
mutable struct Nesterov <: AbstractOptimiser
  eta::Float64
  rho::Float64
  velocity::IdDict
end

Nesterov(η = 0.001, ρ = 0.9) = Nesterov(η, ρ, IdDict())

function apply!(o::Nesterov, x, Δ)
  η, ρ = o.eta, o.rho
  v = get!(() -> zero(x), o.velocity, x)::typeof(x)
  d = @. ρ^2 * v - (1+ρ) * η * Δ
  @. v = ρ*v - η*Δ
  @. Δ = -d
end

"""
    RMSProp(η = 0.001, ρ = 0.9, ϵ = $EPS)

Optimizer using the
[RMSProp](https://www.cs.toronto.edu/~tijmen/csc321/slides/lecture_slides_lec6.pdf)
algorithm. Often a good choice for recurrent networks. Parameters other than learning rate
generally don't need tuning.

# Parameters
- Learning rate (`η`): Amount by which gradients are discounted before updating
                       the weights.
- Momentum (`ρ`): Controls the acceleration of gradient descent in the
                  prominent direction, in effect damping oscillations.

# Examples
```julia
opt = RMSProp()

opt = RMSProp(0.002, 0.95)
```
"""
mutable struct RMSProp <: AbstractOptimiser
  eta::Float64
  rho::Float64
  epsilon::Float64
  acc::IdDict
end
RMSProp(η::Real = 0.001, ρ::Real = 0.9, ϵ::Real = EPS) = RMSProp(η, ρ, ϵ, IdDict())
RMSProp(η::Real, ρ::Real, acc::IdDict) = RMSProp(η, ρ, EPS, acc)

function apply!(o::RMSProp, x, Δ)
  η, ρ = o.eta, o.rho
  acc = get!(() -> zero(x), o.acc, x)::typeof(x)
  @. acc = ρ * acc + (1 - ρ) * Δ * conj(Δ)
  @. Δ *= η / (√acc + o.epsilon)
end

"""
    ADAM(η = 0.001, β::Tuple = (0.9, 0.999), ϵ = $EPS)

[ADAM](https://arxiv.org/abs/1412.6980) optimiser.

# Parameters
- Learning rate (`η`): Amount by which gradients are discounted before updating
                       the weights.
- Decay of momentums (`β::Tuple`): Exponential decay for the first (β1) and the
                                   second (β2) momentum estimate.

# Examples
```julia
opt = ADAM()

opt = ADAM(0.001, (0.9, 0.8))
```
"""
mutable struct ADAM <: AbstractOptimiser
  eta::Float64
  beta::Tuple{Float64,Float64}
  epsilon::Float64
  state::IdDict{Any, Any}
end
ADAM(η::Real = 0.001, β::Tuple = (0.9, 0.999), ϵ::Real = EPS) = ADAM(η, β, ϵ, IdDict())
ADAM(η::Real, β::Tuple, state::IdDict) = ADAM(η, β, EPS, state)

function apply!(o::ADAM, x, Δ)
  η, β = o.eta, o.beta

  mt, vt, βp = get!(o.state, x) do
      (zero(x), zero(x), Float64[β[1], β[2]])
  end :: Tuple{typeof(x),typeof(x),Vector{Float64}}

  @. mt = β[1] * mt + (1 - β[1]) * Δ
  @. vt = β[2] * vt + (1 - β[2]) * Δ * conj(Δ)
  @. Δ =  mt / (1 - βp[1]) / (√(vt / (1 - βp[2])) + o.epsilon) * η
  βp .= βp .* β

  return Δ
end

"""
    RADAM(η = 0.001, β::Tuple = (0.9, 0.999), ϵ = $EPS)

[Rectified ADAM](https://arxiv.org/abs/1908.03265) optimizer.

# Parameters
- Learning rate (`η`): Amount by which gradients are discounted before updating
                       the weights.
- Decay of momentums (`β::Tuple`): Exponential decay for the first (β1) and the
                                   second (β2) momentum estimate.

# Examples
```julia
opt = RADAM()

opt = RADAM(0.001, (0.9, 0.8))
```
"""
mutable struct RADAM <: AbstractOptimiser
  eta::Float64
  beta::Tuple{Float64,Float64}
  epsilon::Float64
  state::IdDict{Any, Any}
end
RADAM(η::Real = 0.001, β::Tuple = (0.9, 0.999), ϵ::Real = EPS) = RADAM(η, β, ϵ, IdDict())
RADAM(η::Real, β::Tuple, state::IdDict) = RADAM(η, β, EPS, state)

function apply!(o::RADAM, x, Δ)
  η, β = o.eta, o.beta
  ρ∞ = 2/(1-β[2])-1

  mt, vt, βp, t = get!(o.state, x) do
      (zero(x), zero(x), Float64[β[1], β[2]], Ref(1))
  end :: Tuple{typeof(x),typeof(x),Vector{Float64},Base.RefValue{Int}}

  @. mt = β[1] * mt + (1 - β[1]) * Δ
  @. vt = β[2] * vt + (1 - β[2]) * Δ * conj(Δ)
  ρ = ρ∞ - 2t[] * βp[2] / (1 - βp[2])
  if ρ > 4
    r = sqrt((ρ-4)*(ρ-2)*ρ∞/((ρ∞-4)*(ρ∞-2)*ρ))
    @. Δ =  mt / (1 - βp[1]) / (√(vt / (1 - βp[2])) + o.epsilon) * η * r
  else
    @. Δ =  mt / (1 - βp[1]) * η
  end
  βp .= βp .* β
  t[] += 1

  return Δ
end

"""
    AdaMax(η = 0.001, β::Tuple = (0.9, 0.999), ϵ = $EPS)

[AdaMax](https://arxiv.org/abs/1412.6980) is a variant of ADAM based on the ∞-norm.

# Parameters
- Learning rate (`η`): Amount by which gradients are discounted before updating
                       the weights.
- Decay of momentums (`β::Tuple`): Exponential decay for the first (β1) and the
                                   second (β2) momentum estimate.

# Examples
```julia
opt = AdaMax()

opt = AdaMax(0.001, (0.9, 0.995))
```
"""
mutable struct AdaMax <: AbstractOptimiser
  eta::Float64
  beta::Tuple{Float64,Float64}
  epsilon::Float64
  state::IdDict{Any, Any}
end
AdaMax(η::Real = 0.001, β::Tuple = (0.9, 0.999), ϵ::Real = EPS) = AdaMax(η, β, ϵ, IdDict())
AdaMax(η::Real, β::Tuple, state::IdDict) = AdaMax(η, β, EPS, state)

function apply!(o::AdaMax, x, Δ)
  η, β = o.eta, o.beta

  mt, ut, βp = get!(o.state, x) do
      (zero(x), zero(x), Float64[β[1], β[2]])
  end :: Tuple{typeof(x),typeof(x),Vector{Float64}}

  @. mt = β[1] * mt + (1 - β[1]) * Δ
  @. ut = max(β[2] * ut, abs(Δ))
  @. Δ = (η/(1 - βp[1])) * mt/(ut + o.epsilon)
  βp .= βp .* β

  return Δ
end

"""
    OADAM(η = 0.0001, β::Tuple = (0.5, 0.9), ϵ = $EPS)

[OADAM](https://arxiv.org/abs/1711.00141) (Optimistic ADAM)
is a variant of ADAM adding an "optimistic" term suitable for adversarial training.

# Parameters
- Learning rate (`η`): Amount by which gradients are discounted before updating
                       the weights.
- Decay of momentums (`β::Tuple`): Exponential decay for the first (β1) and the
                                   second (β2) momentum estimate.

# Examples
```julia
opt = OADAM()

opt = OADAM(0.001, (0.9, 0.995))
```
"""
mutable struct OADAM <: AbstractOptimiser
  eta::Float64
  beta::Tuple{Float64,Float64}
  epsilon::Float64
  state::IdDict{Any, Any}
end
OADAM(η::Real = 0.001, β::Tuple = (0.5, 0.9), ϵ::Real = EPS) = OADAM(η, β, ϵ, IdDict())
OADAM(η::Real, β::Tuple, state::IdDict) = RMSProp(η, β, EPS, state)

function apply!(o::OADAM, x, Δ)
  η, β = o.eta, o.beta

  mt, vt, Δ_, βp = get!(o.state, x) do
      (zero(x), zero(x), zero(x), Float64[β[1], β[2]])
  end :: Tuple{typeof(x),typeof(x),typeof(x),Vector{Float64}}

  @. mt = β[1] * mt + (1 - β[1]) * Δ
  @. vt = β[2] * vt + (1 - β[2]) * Δ * conj(Δ)
  @. Δ = -Δ_
  @. Δ_ = η * mt / (1 - βp[1]) / (√(vt / (1 - βp[2])) + o.epsilon)
  @. Δ += 2Δ_
  βp .= βp .* β

  return Δ
end

"""
    ADAGrad(η = 0.1, ϵ = $EPS)

[ADAGrad](http://www.jmlr.org/papers/volume12/duchi11a/duchi11a.pdf) optimizer. It has
parameter specific learning rates based on how frequently it is updated.
Parameters don't need tuning.

# Parameters
- Learning rate (`η`): Amount by which gradients are discounted before updating
                       the weights.

# Examples
```julia
opt = ADAGrad()

opt = ADAGrad(0.001)
```
"""
mutable struct ADAGrad <: AbstractOptimiser
  eta::Float64
  epsilon::Float64
  acc::IdDict
end
ADAGrad(η::Real = 0.1, ϵ::Real = EPS) = ADAGrad(η, ϵ, IdDict())
ADAGrad(η::Real, state::IdDict) = ADAGrad(η, EPS, state)

function apply!(o::ADAGrad, x, Δ)
  η = o.eta
  acc = get!(() -> fill!(similar(x), o.epsilon), o.acc, x)::typeof(x)
  @. acc += Δ * conj(Δ)
  @. Δ *= η / (√acc + o.epsilon)
end

"""
    ADADelta(ρ = 0.9, ϵ = $EPS)

[ADADelta](https://arxiv.org/abs/1212.5701) is a version of ADAGrad adapting its learning
rate based on a window of past gradient updates.
Parameters don't need tuning.

# Parameters
- Rho (`ρ`): Factor by which the gradient is decayed at each time step.

# Examples
```julia
opt = ADADelta()

opt = ADADelta(0.89)
```
"""
mutable struct ADADelta <: AbstractOptimiser
  rho::Float64
  epsilon::Float64
  state::IdDict{Any, Any}
end
ADADelta(ρ::Real = 0.9, ϵ::Real = EPS) = ADADelta(ρ, ϵ, IdDict())
ADADelta(ρ::Real, state::IdDict) = ADADelta(ρ, EPS, state)

function apply!(o::ADADelta, x, Δ)
  ρ = o.rho
  acc, Δacc = get!(() -> (zero(x), zero(x)), o.state, x)::NTuple{2,typeof(x)}
  @. acc = ρ * acc + (1 - ρ) * Δ * conj(Δ)
  # DON'T remove epsilon from numerator
  # or even out of the square roots
  @. Δ *= √(Δacc + o.epsilon) / √(acc + o.epsilon)
  @. Δacc = ρ * Δacc + (1 - ρ) * Δ * conj(Δ)
  return Δ
end

"""
    AMSGrad(η = 0.001, β::Tuple = (0.9, 0.999), ϵ = $EPS)

The [AMSGrad](https://openreview.net/forum?id=ryQu7f-RZ) version of the ADAM
optimiser. Parameters don't need tuning.

# Parameters
- Learning rate (`η`): Amount by which gradients are discounted before updating
                       the weights.
- Decay of momentums (`β::Tuple`): Exponential decay for the first (β1) and the
                                   second (β2) momentum estimate.

# Examples
```julia
opt = AMSGrad()

opt = AMSGrad(0.001, (0.89, 0.995))
```
"""
mutable struct AMSGrad <: AbstractOptimiser
  eta::Float64
  beta::Tuple{Float64, Float64}
  epsilon::Float64
  state::IdDict{Any, Any}
end
AMSGrad(η::Real = 0.001, β = (0.9, 0.999), ϵ::Real = EPS) = AMSGrad(η, β, ϵ, IdDict())
AMSGrad(η::Real, β::Tuple, state::IdDict) = AMSGrad(η, β, EPS, state)

function apply!(o::AMSGrad, x, Δ)
  η, β = o.eta, o.beta

  mt, vt, v̂t = get!(o.state, x) do
    (fill!(similar(x), o.epsilon), fill!(similar(x), o.epsilon), fill!(similar(x), o.epsilon))
  end :: NTuple{3,typeof(x)}

  @. mt = β[1] * mt + (1 - β[1]) * Δ
  @. vt = β[2] * vt + (1 - β[2]) * Δ ^ 2
  @. v̂t = max(v̂t, vt)
  @. Δ = η * mt / (√v̂t + o.epsilon)
end

"""
    NADAM(η = 0.001, β::Tuple = (0.9, 0.999), ϵ = $EPS)

[NADAM](https://openreview.net/forum?id=OM0jvwB8jIp57ZJjtNEZ) is a Nesterov variant of ADAM.
Parameters don't need tuning.

# Parameters
- Learning rate (`η`): Amount by which gradients are discounted before updating
                       the weights.
- Decay of momentums (`β::Tuple`): Exponential decay for the first (β1) and the
                                   second (β2) momentum estimate.

# Examples
```julia
opt = NADAM()

opt = NADAM(0.002, (0.89, 0.995))
```
"""
mutable struct NADAM <: AbstractOptimiser
  eta::Float64
  beta::Tuple{Float64, Float64}
  epsilon::Float64
  state::IdDict{Any, Any}
end
NADAM(η::Real = 0.001, β = (0.9, 0.999), ϵ::Real = EPS) = NADAM(η, β, ϵ, IdDict())
NADAM(η::Real, β::Tuple, state::IdDict) = NADAM(η, β, EPS, state)

function apply!(o::NADAM, x, Δ)
  η, β = o.eta, o.beta

  mt, vt, βp = get!(o.state, x) do
    (zero(x), zero(x), Float64[o.beta[1], o.beta[2]])
  end :: Tuple{typeof(x),typeof(x),Vector{Float64}}
  β1p, β2p = βp

  @. mt = β[1] * mt + (1 - β[1]) * Δ
  @. vt = β[2] * vt + (1 - β[2]) * Δ * conj(Δ)
  @. Δ = (β[1] * mt / (1 - β[1] * β1p) + (1 - β[1]) * Δ / (1 - β1p)) / (√(vt * β[2] / (1 - β2p)) + o.epsilon) * η
  βp .= βp .* β

  return Δ
end

"""
    ADAMW(η = 0.001, β::Tuple = (0.9, 0.999), decay = 0)

[ADAMW](https://arxiv.org/abs/1711.05101) is a variant of ADAM fixing (as in repairing) its
weight decay regularization.

# Parameters
- Learning rate (`η`): Amount by which gradients are discounted before updating
                       the weights.
- Decay of momentums (`β::Tuple`): Exponential decay for the first (β1) and the
                                   second (β2) momentum estimate.
- `decay`: Decay applied to weights during optimisation.

# Examples
```julia
opt = ADAMW()

opt = ADAMW(0.001, (0.89, 0.995), 0.1)
```
"""
ADAMW(η = 0.001, β = (0.9, 0.999), decay = 0) =
  Optimiser(ADAM(η, β), WeightDecay(decay))

"""
    AdaBelief(η = 0.001, β::Tuple = (0.9, 0.999), ϵ = $EPS)

The [AdaBelief](https://arxiv.org/abs/2010.07468) optimiser is a variant of the well-known
ADAM optimiser.

# Parameters
- Learning rate (`η`): Amount by which gradients are discounted before updating
                       the weights.
- Decay of momentums (`β::Tuple`): Exponential decay for the first (β1) and the
                                   second (β2) momentum estimate.

# Examples
```julia
opt = AdaBelief()

opt = AdaBelief(0.001, (0.9, 0.8))
```
"""
mutable struct AdaBelief <: AbstractOptimiser
  eta::Float64
  beta::Tuple{Float64,Float64}
  epsilon::Float64
  state::IdDict{Any, Any}
end
AdaBelief(η::Real = 0.001, β = (0.9, 0.999), ϵ::Real = EPS) = AdaBelief(η, β, ϵ, IdDict())
AdaBelief(η::Real, β::Tuple, state::IdDict) = AdaBelief(η, β, EPS, state)

function apply!(o::AdaBelief, x, Δ)
  η, β = o.eta, o.beta
  mt, st = get!(() -> (zero(x), zero(x)), o.state, x)::NTuple{2,typeof(x)}
  @. mt = β[1] * mt + (1 - β[1]) * Δ
  @. st = β[2] * st + (1 - β[2]) * (Δ - mt) * conj(Δ - mt)
  @. Δ =  η * mt / (√(st) + o.epsilon)
  return Δ
end

# Compose optimizers

"""
    Optimiser(a, b, c...)

Combine several optimisers into one; each optimiser produces a modified gradient
that will be fed into the next, and this is finally applied to the parameter as
usual.
"""
mutable struct Optimiser <: AbstractOptimiser
  os::Vector{Any}
end

Optimiser(opts::AbstractOptimiser...) = Optimiser(Any[opts...])

@forward Optimiser.os Base.getindex, Base.first, Base.last, Base.lastindex, Base.push!, Base.setindex!
@forward Optimiser.os Base.iterate

Base.getindex(c::Optimiser, i::AbstractArray) = Optimiser(c.os[i]...)

function apply!(o::Optimiser, x, Δ)
  for opt in o.os
    Δ = apply!(opt, x, Δ)
  end
  return Δ
end

"""
    InvDecay(γ = 0.001)

Apply inverse time decay to an optimiser, so that the effective step size at
iteration `n` is `eta / (1 + γ * n)` where `eta` is the initial step size.
The wrapped optimiser's step size is not modified.

See also the [Scheduling Optimisers](@ref) section of the docs
for more general scheduling techniques.

# Examples

`InvDecay` is typically composed  with other optimizers
as the last transformation of the gradient:

```julia
# Inverse decay of the learning rate
# with starting value 0.001 and decay coefficient 0.01.
opt = Optimiser(Adam(1f-3), InvDecay(1f-2))
```
"""
mutable struct InvDecay <: AbstractOptimiser
  gamma::Float64
  state::IdDict{Any, Int}
end

InvDecay(γ = 0.001) = InvDecay(γ, IdDict{Any, Int}())

function apply!(o::InvDecay, x, Δ)
  γ = o.gamma
  n = get!(o.state, x, 1)
  Δ .*= 1 / (1 + γ * n)
  o.state[x] = n + 1
  return Δ
end

"""
    ExpDecay(η = 0.001, decay = 0.1, decay_step = 1000, clip = 1e-4, start = 1)

Discount the learning rate `η` by the factor `decay` every `decay_step` steps till
a minimum of `clip`.

# Parameters
- Learning rate (`η`): Amount by which gradients are discounted before updating
                       the weights.
- `decay`: Factor by which the learning rate is discounted.
- `decay_step`: Schedule decay operations by setting the number of steps between
                two decay operations.
- `clip`: Minimum value of learning rate.
- 'start': Step at which the decay starts.

See also the [Scheduling Optimisers](@ref) section of the docs
for more general scheduling techniques.

# Examples

`ExpDecay` is typically composed  with other optimizers
as the last transformation of the gradient:
```julia
opt = Optimiser(ADAM(), ExpDecay())
```
"""
mutable struct ExpDecay <: AbstractOptimiser
  eta::Float64
  decay::Float64
  step::Int64
  clip::Float64
  start::Int64
  current::IdDict
end

ExpDecay(opt = 0.001, decay = 0.1, decay_step = 1000, clip = 1e-4, start = 0) =
  ExpDecay(opt, decay, decay_step, clip, start, IdDict())

function apply!(o::ExpDecay, x, Δ)
  η, s, decay, start = o.eta, o.step, o.decay, o.start
  n = o.current[x] = get(o.current, x, 0) + 1
  if n > start && n % s == 0 && count(x -> x > start && x % s == 0, values(o.current)) == 1
    η = max(η * decay, o.clip)
    o.eta = η
  end
  @. Δ *= η
end

"""
    WeightDecay(λ = 0)

Decay weights by ``λ``.
Typically composed  with other optimizers as the first transformation to the gradient,
making it equivalent to adding ``L_2`` regularization
with coefficient  ``λ`` to the loss.

# Examples

```julia
opt = Optimiser(WeightDecay(1f-4), ADAM())
```
"""
mutable struct WeightDecay <: AbstractOptimiser
  wd::Real
end

WeightDecay() = WeightDecay(0)

function apply!(o::WeightDecay, x, Δ)
  wd = o.wd
  @. Δ += wd * x
end

"""
    ClipValue(thresh)

Clip gradients when their absolute value exceeds `thresh`.
"""
mutable struct ClipValue{T} <: AbstractOptimiser
    thresh::T
end

apply!(o::ClipValue, x, Δ) = clamp!(Δ, -o.thresh, o.thresh)

"""
    ClipNorm(thresh)

Clip gradients when their L2 norm exceeds `thresh`.
"""
mutable struct ClipNorm{T} <: AbstractOptimiser
    thresh::T
end

function apply!(o::ClipNorm, x, Δ)
    Δnrm = norm(Δ)
    if Δnrm > o.thresh
        rmul!(Δ, o.thresh / Δnrm)
    end
    return Δ
end


text:
#md# This animation shows how the initial layout is generated.
using WordCloud
stopwords = WordCloud.stopwords_en ∪ ["said"]
textfile = pkgdir(WordCloud)*"/res/alice.txt"
wc = wordcloud(
    processtext(open(textfile), stopwords=stopwords, maxnum=300),
    masksize = (300, 200),
    outline = 3,
    angles = 0:90,
    state = initwords!)
#md# ### uniform style
gifdirectory = "animation1/uniform"
setpositions!(wc, :, (-1000,-1000))
@record gifdirectory overwrite=true filter=i->i%(2^(i÷100))==0 placewords!(wc, style=:uniform)
#md# ![](animation1/uniform/animation.gif)
#md# ### gathering style
gifdirectory = "animation1/gathering"
setpositions!(wc, :, (-1000,-1000))
@record gifdirectory overwrite=true filter=i->i%(2^(i÷100))==0 placewords!(wc, style=:gathering)
#md# ![](animation1/gathering/animation.gif)
#md#
println("results are saved in animation1")
wc
#eval# runexample(:animation1)


text:
using FlowPattern

function copy_files_from_external_folder(source_dir, app_dir)
    target_externalpath = joinpath(app_dir, "external")
    source_external_path = joinpath(source_dir, "external")
    mkpath(target_externalpath)
    for file in readdir(source_external_path)
        cp(joinpath(source_external_path, file), joinpath(target_externalpath, file))
    end
end

copy_files_from_external_folder(
    @__DIR__,
    join((@__DIR__, "Compiled")),
)

push!(ARGS, "--quiet")
FlowPattern.julia_main()

max_stars_repo_name: arthursoprana/FlowPattern

max_issues_repo_name: arthursoprana/FlowPattern

max_forks_repo_name: arthursoprana/FlowPattern


text:
using TypeClasses
using Test
using DataTypesBasic
using Suppressor

splitln(str) = split(strip(str), "\n")

# Combine
# =======

a = Callable(x -> "hello $x")
b = Callable(x -> "!")

(a ⊕ b)(:Albert)

# FunctorApplicativeMonad
# =======================

g = Callable(x -> x*2)
f = Callable(x -> x*x)

# just function composition
fg = map(g) do x2
  f(x2)
end
@test fg(3) == f(g(3))

fPLUSg = mapn(f, g) do x1, x2
  x1 + x2
end

@test fPLUSg(3) == f(3) + g(3)

fPRODg = @syntax_flatmap begin
  x1 = f
  x2 = g
  @pure x1 * x2
end

@test fPRODg(3) == f(3) * g(3)

# FlipTypes
# =========

# this works because Callable implements `ap`

a = Callable.([x -> x, y -> 2y, z -> z*z])
@test flip_types(a)(3) == [3, 6, 9]


text:
#!/usr/bin/env julia

using RobotOS
using JuMP, Ipopt
using JLD
@rosimport barc.msg: ECU, Z_KinBkMdl, mpcSol
@rosimport std_msgs.msg: Float32
rostypegen()
using barc.msg
using std_msgs.msg
include("modules/mpcPathFollowing.jl")

function euclideanDistance(s1::Array{Float64},s2::Array{Float64})
	dx = s1[1] - s2[1]
	dy = s1[2] - s2[2]
	return sqrt( dx^2 + dy^2 )
end

function updateState(msg::Z_KinBkMdl,
					sCurr::Array{Float64},
					sPrev::Array{Float64},
					dist::Array{Float64},
					sRef::Array{Float64,2},
					uRef::Array{Float64,2},
					sp10Ref::Array{Float64,2},
					up10Ref::Array{Float64,2},
					dp10Ref::Array{Float64},
					N::Int64)
    # get state
    sPrev[:] = copy(sCurr)
    sCurr[:] = [msg.x , msg.y , msg.psi , msg.v]

    # compute distance travelled in previous time step
    delD  	   = euclideanDistance(sCurr, sPrev)
    dist[:]    = dist[1] + delD
    idx 	   = indmin( (dp10Ref - dist[1]).^2 )

	sRef[:,:] = copy( sp10Ref[:,idx:idx+N] )
	uRef[:,:] = copy( up10Ref[:,idx:idx+N-1] )
end

function main()
    # initiate node, set up publisher / subscriber topics
    println("loading parameters ...")
    s0 			= get_param("/initial_state")
    sF 			= get_param("/final_state")
    L 			= get_param("/vehicle_length")
    Ts 			= get_param("/simulation_time_step")
    N 			= get_param("/horizon_preview")
    Q 			= get_param("/state_penalty_matrix")
    R 			= get_param("/input_penalty_matrix")
    stateDim 	= size(s0)[1]
    inputDim 	= 2
    Q 			= reshape(Q,(stateDim, stateDim))
    R 			= reshape(R,(inputDim, inputDim))
	mdlParams   = MdlParams(L,Ts)
	mpcParams   = MpcParams(N,Q,R)
	sol         = MpcSol()

	println("running first optimization with a horizon preview of N = ", N, " steps ...")
	mdl    		= MpcModel(mpcParams,mdlParams)

	println("loading pre-computed solution ...")
    home        = ENV["HOME"]
	filePath 	= string(home,"/barc/workspace/src/barc/src/control/parking/data")
	data 		= load("$filePath/parkPath.jld")
	sp10 		= data["sp10"]
	up10 		= data["up10"]
	scaleTime10 = data["scaleTime10"]
	sp10Ref 	= data["sp10Ref"]
	up10Ref 	= data["up10Ref"]
	dp10Ref 	= data["dp10Ref"]

	sp10Rest 	= sp10[:,end]*ones(1,N+1)
	up10Rest 	= [0,0]*ones(1,N)
	sp10Ref 	= hcat(sp10Ref, sp10Rest)
	up10Ref 	= hcat(up10Ref, up10Rest)
	sCurr 		= copy(s0)
	sPrev 		= copy(s0)
	sRef 		= copy(sp10Ref[:,1:N+1] )
	uRef 		= copy(up10Ref[:,1:N] )
	dist  		= [0.0]
	solveMPC 	= true

	println("initilizing node ...")
    init_node("controller")
    p1 = Publisher("controller", ECU, queue_size=1)
    p2 = Publisher("mpcSolution", mpcSol, queue_size=1)
    s1 = Subscriber{Z_KinBkMdl}("plant", updateState, (sCurr,sPrev,dist,sRef,uRef,sp10Ref,up10Ref,dp10Ref,N), queue_size=1)
    loop_rate = Rate( Int(1.0/Ts) )

	println("running mpc ...")
    while ! is_shutdown()

    	if euclideanDistance(sCurr, sF) > 0.1 && solveMPC
    		SolveMpcProblem(mdl, sol, sCurr, sRef, uRef)
    		publish(p1, ECU( sol.acc, sol.df) )
    		publish(p2, mpcSol( sol.zOL[1,:], sol.zOL[2,:], sol.zOL[3,:] , sol.zOL[4,:], sol.uOL[1,:], sol.uOL[2,:]) )
    		println(sol.solverStatus)
    	else
    		solveMPC = false
    		cmd = ECU(-sCurr[4], 0)
    		publish(p1, cmd )
    	end

        rossleep(loop_rate)
    end
end

if ! isinteractive()
    main()
end

max_stars_repo_name: treyfortmuller/barc

max_issues_repo_name: treyfortmuller/barc


text:
# # Paper examples

using Convex, ECOS

# Summation.
println("Summation example")
x = Variable();
e = 0;
@time begin
  for i = 1:1000
    global e
    e += x;
  end
  p = minimize(e, x>=1);
end
@time solve!(p, () -> ECOS.Optimizer(verbose=0))

# Indexing.
println("Indexing example")
x = Variable(1000, 1);
e = 0;
@time begin
  for i = 1:1000
    global e
    e += x[i];
  end
  p = minimize(e, x >= ones(1000, 1));
end
@time solve!(p, () -> ECOS.Optimizer(verbose=0))

# Matrix constraints.
println("Matrix constraint example")
n, m, p = 100, 100, 100
X = Variable(m, n);
A = randn(p, m);
b = randn(p, n);
@time begin
  p = minimize(norm(vec(X)), A * X == b);
end
@time solve!(p, () -> ECOS.Optimizer(verbose=0))

# Transpose.
println("Transpose example")
X = Variable(5, 5);
A = randn(5, 5);
@time begin
  p = minimize(norm2(X - A), X' == X);
end
@time solve!(p, () -> ECOS.Optimizer(verbose=0))

n = 3
A = randn(n, n);
#@time begin
  X = Variable(n, n);
  p = minimize(norm(vec(X' - A)), X[1,1] == 1);
  solve!(p, () -> ECOS.Optimizer(verbose=0))
#end


text:
function Request(service::Service, operation_name::AbstractString)
    request_handle_ref = Ref{Ptr{Cvoid}}(C_NULL)
    err = blpapi_Service_createRequest(service.handle, request_handle_ref, operation_name)
    error_check(err, "Failed to create request")
    return Request(request_handle_ref[], service)
end

Request(session::Session, service_name::AbstractString, operation_name::AbstractString) = Request(Service(session, service_name), operation_name)

function Element(request::Request)
    return Element(blpapi_Request_elements(request.handle), request)
end

function Base.getindex(request::Request, element_name::AbstractString)
    elements = Element(request)
    return elements[element_name]
end

function Base.setindex!(request::Request, val, element_name::AbstractString)
    elements = Element(request)
    elements[element_name] = val
end

function send_request(request::Request, queue::EventQueue=EventQueue()) :: Tuple{EventQueue, CorrelationId}
    correlation_id_ref = Ref(CorrelationId())
    session = request.service.session
    err = blpapi_Session_sendRequest(session.handle, request.handle, correlation_id_ref, queue.handle)
    error_check(err, "Failed to send request")
    return queue, correlation_id_ref[]
end

function send_request(f::Function, session::Session, service_name::AbstractString, operation_name::AbstractString, queue::EventQueue=EventQueue()) :: Tuple{EventQueue, CorrelationId}
    request = Request(session, service_name, operation_name)

    # inspect result schema
    # elements = Element(req)
    # elements_schema = SchemaElementDefinition(elements)

    f(request)
    return send_request(request, queue)
end

function for_each_response_message_element(f::Function, queue::EventQueue, corr_id::CorrelationId; timeout_milliseconds::Integer=UInt32(0), verbose::Bool=false)
    while true

        response_event = next_event(queue, timeout_milliseconds=timeout_milliseconds)

        try
            if response_event.event_type == BLPAPI_EVENTTYPE_TIMEOUT
                throw(BLPTimeoutException(timeout_milliseconds))

            elseif response_event.event_type == BLPAPI_EVENTTYPE_REQUEST_STATUS
                throw(BLPResponseException("Got response with type $(response_event.event_type): $response_event"))
            end

            if !(response_event.event_type == BLPAPI_EVENTTYPE_RESPONSE || response_event.event_type == BLPAPI_EVENTTYPE_PARTIAL_RESPONSE)
                throw(BLPUnknownException("Tried to handle non-response event of type $(response_event.event_type)."))
            end

            # process BLPAPI_EVENTTYPE_RESPONSE or BLPAPI_EVENTTYPE_PARTIAL_RESPONSE
            verbose && @info("Reading messages from event $(response_event.event_type)")
            for message in each_message(response_event)
                check_has_correlation_id(message, corr_id)
                element = Element(message)

                if verbose
                    println("Reponse Element Schema")
                    println(SchemaElementDefinition(element))
                end

                f(element)
            end

            # check if response is complete
            if response_event.event_type == BLPAPI_EVENTTYPE_RESPONSE
                verbose && @info("Finished reading events")
                break
            end

        finally
            # destroy this event early, before GC
            destroy!(response_event)
        end
    end

    nothing
end

"""
    parse_response_as(::Type{T}, queue::EventQueue, corr_id::CorrelationId; timeout_milliseconds::Integer=UInt32(0), verbose::Bool=false) :: Vector{T} where {T}

Parses all response messages as `Vector{T}`, applying `T(element)` for each `element` read from the response.
"""
function parse_response_as(::Type{T}, queue::EventQueue, corr_id::CorrelationId; timeout_milliseconds::Integer=UInt32(0), verbose::Bool=false) where {T}
    result = Vector{T}()
    for_each_response_message_element(queue, corr_id; timeout_milliseconds=timeout_milliseconds, verbose=verbose) do element
        push!(result, T(element))
    end
    return result
end

max_stars_repo_licenses: OpenSSL MIT

max_issues_repo_licenses: OpenSSL MIT

max_forks_repo_licenses: OpenSSL MIT


text:
module BoundaryDiscreteModelsTests

using Test
using Gridap.Geometry
using Gridap.ReferenceFEs
using Gridap.FESpaces
using Gridap.CellData
using Gridap.Arrays
using Gridap.TensorValues

domain = (0,1,0,1,0,1)
cells = (3,3,3)
bgmodel = CartesianDiscreteModel(domain,cells)

labels = get_face_labeling(bgmodel)
bgface_to_mask = get_face_mask(labels,[22,23],2)

model = BoundaryDiscreteModel(Polytope{2},bgmodel,bgface_to_mask)

trian = Triangulation(model)
btrian = BoundaryTriangulation(model)
strian = SkeletonTriangulation(model)

@test get_background_triangulation(trian) === Triangulation(bgmodel)
@test get_background_triangulation(btrian) === trian
@test get_background_triangulation(strian) === trian
@test get_background_triangulation(Triangulation(model.model)) === Triangulation(model.model)

reffe = ReferenceFE(lagrangian,Float64,1)
V = FESpace(model,reffe)
bgV = FESpace(bgmodel,reffe)

vh = FEFunction(V,rand(num_free_dofs(V)))
bgvh = FEFunction(bgV,rand(num_free_dofs(bgV)))

@test get_triangulation(vh) === Triangulation(model)
@test get_background_triangulation(get_triangulation(vh)) === Triangulation(bgmodel)
@test get_background_triangulation(get_triangulation(vh)) === get_triangulation(bgvh)

Ω = Triangulation(bgmodel)
Γ = trian
∂Γ = btrian
sΓ = strian

dΩ = Measure(Ω,3)
dΓ = Measure(Γ,3)
d∂Γ = Measure(∂Γ,3)
dsΓ = Measure(sΓ,3)

n_Γ = get_normal_vector(Γ)
n_∂Γ = get_normal_vector(∂Γ)

@test ∑(∫(1)dΓ) ≈ 2.0
@test ∑(∫(1)d∂Γ) ≈ 6.0
@test ∑(∫(1)dsΓ) ≈ 9.0

@test ∑(∫(n_Γ⋅VectorValue(1.0,0.0,0.0))dΓ) ≈ 0.0
@test ∑(∫(n_Γ⋅VectorValue(0.0,1.0,0.0))dΓ) ≈ -1.0
@test ∑(∫(n_Γ⋅VectorValue(0.0,0.0,1.0))dΓ) ≈ 1.0
@test ∑(∫(n_Γ⋅VectorValue(0.0,0.0,1.0))d∂Γ) ≈ 3.0
@test ∑(∫(n_∂Γ⋅VectorValue(0.0,0.0,1.0))d∂Γ) ≈ -1.0

domain   = (-1,1,-1,1,-1,1)
cells    = (1,1,1)
bgmodel  = CartesianDiscreteModel(domain,cells)
model    = BoundaryDiscreteModel(Polytope{2},bgmodel,Int32[3,4,5,6,1,2])
topo     = get_grid_topology(model)
cell_vertices = get_faces(topo,2,0)
@test cell_vertices[5] == [1,2,3,4]
@test cell_vertices[6] == [5,6,7,8]
@test cell_vertices[1] == [1,2,5,6]
@test cell_vertices[2] == [3,4,7,8]
@test cell_vertices[3] == [1,3,5,7]
@test cell_vertices[4] == [2,4,6,8]
end


text:
using SemiparametricAUC
using DataArrays
using Base.Test

# write your own tests here
@test 1 == 1

println("Testing Som's Package")
include("calculate-auc.jl")

@test calculate_auc(ya = DataArray([2,3]), yb = DataArray([3,2,1])) == (1.5468749999999993,0.6931471805599452)


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

struct __enum_EnumMatrixDeviceType <: ProtoEnum
    UNDETECTED::Int32
    CREATOR::Int32
    VOICE::Int32
    __enum_EnumMatrixDeviceType() = new(0,1,2)
end #struct __enum_EnumMatrixDeviceType
const EnumMatrixDeviceType = __enum_EnumMatrixDeviceType()

mutable struct LedValue <: ProtoType
    red::UInt32
    green::UInt32
    blue::UInt32
    white::UInt32
    LedValue(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct LedValue

mutable struct EverloopImage <: ProtoType
    led::Base.Vector{LedValue}
    everloop_length::Int32
    EverloopImage(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct EverloopImage

mutable struct ServoParams <: ProtoType
    pin::UInt32
    angle::UInt32
    ServoParams(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ServoParams

struct __enum_GpioParams_EnumMode <: ProtoEnum
    INPUT::Int32
    OUTPUT::Int32
    __enum_GpioParams_EnumMode() = new(0,1)
end #struct __enum_GpioParams_EnumMode
const GpioParams_EnumMode = __enum_GpioParams_EnumMode()

mutable struct GpioParams <: ProtoType
    pin::UInt32
    mode::Int32
    value::UInt32
    values::UInt32
    GpioParams(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct GpioParams

mutable struct MicArrayParams <: ProtoType
    gain::Int32
    azimutal_angle::Float32
    polar_angle::Float32
    radial_distance_mm::Float32
    sound_speed_mmseg::Float32
    sampling_frequency_hz::UInt32
    MicArrayParams(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct MicArrayParams

struct __enum_WakeWordParams_MicChannel <: ProtoEnum
    channel0::Int32
    channel1::Int32
    channel2::Int32
    channel3::Int32
    channel4::Int32
    channel5::Int32
    channel6::Int32
    channel7::Int32
    channel8::Int32
    __enum_WakeWordParams_MicChannel() = new(0,1,2,3,4,5,6,7,8)
end #struct __enum_WakeWordParams_MicChannel
const WakeWordParams_MicChannel = __enum_WakeWordParams_MicChannel()

mutable struct WakeWordParams <: ProtoType
    wake_word::AbstractString
    channel::Int32
    lm_path::AbstractString
    dic_path::AbstractString
    enable_verbose::Bool
    stop_recognition::Bool
    WakeWordParams(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct WakeWordParams

mutable struct MatrixDeviceParams <: ProtoType
    device_type::Int32
    MatrixDeviceParams(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct MatrixDeviceParams

export EnumMatrixDeviceType, LedValue, EverloopImage, ServoParams, GpioParams_EnumMode, GpioParams, MicArrayParams, WakeWordParams_MicChannel, WakeWordParams, MatrixDeviceParams


text:
module JavaCoreNLP

export
    Annotation,
    StanfordCoreNLP,
    annotate!

include("core.jl")

end # module


text:
for t in (:ICD10AM, :ICD10CA, :ICD10GM)
  @eval begin
    struct $t <: AbstractICD10
      level::Int8
      data::NTuple{5,UInt8}
    end

    function $t(str::AbstractString, validateinput = false)
      punct = occursin(".", str) ? true : false
      if validateinput
        validICD10input(str) || throw(
          DomainError(
            str,
            "ICD-10 codes should have format `ANN[.][N[N]]` where `A` is letter A-Z, `N` is decimal digit, and parts in brackets are optional",
          ),
        )
      end
      level = punct ? length(str) - 1 : length(str)
      ch4 = punct ? 5 : 4

      if level == 3
        return $t(level, (NTuple{3,UInt8}.(str)..., UInt8(0), UInt8(0)))
      elseif level == 4
        return $t(level, (NTuple{3,UInt8}.(str)..., UInt8(str[ch4]), UInt8(0)))
      elseif level == 5
        return $t(level, (NTuple{3,UInt8}.(str)..., UInt8(str[ch4]), UInt8(str[ch4+1])))
      else
        throw(
          DomainError(str, "ICD-10 codes should be 3-5 characters excluding the period"),
        )
      end
    end

    function $t(icd::T) where {T<:AbstractICD10}
      $t(icd.level, icd.data)
    end
  end
end

"""
    ICD10AM

Australian version of ICD-10
"""
ICD10AM

"""
    ICD10CA

Canadian version of ICD-10
"""
ICD10CA

# """
#     ICD10CM

# United States version of ICD-10
# """
# ICD10CM

"""
    ICD10GM

German version of ICD-10
"""
ICD10GM

"""
    ICD10AM(str::String, validateinput = false)

Create an ICD10AM code from a string.

ICD10AM codes have the form `ANN[.][N[N]]` where `A` is a letter A-Z and
`N` is
a digit 0-9. Parts in brackets are optional.

If `validateinput=true` then the input string is checked for valid format using
a regex.
"""
ICD10AM(str::AbstractString, validateinput = false)

"""
    ICD10CA(str::String, validateinput = false)

Create an ICD10CA code from a string.

ICD10CA codes have the form `ANN[.][N[N]]` where `A` is a letter A-Z and
`N` is a digit 0-9. Parts in brackets are optional.

If `validateinput=true` then the input string is checked for valid format using
a regex.
"""
ICD10CA(str::AbstractString, validateinput = false)

"""
    ICD10GM(str::String, validateinput = false)

Create an ICD10GM code from a string.

ICD10GM codes have the form `ANN[.][N[N]]` where `A` is a letter A-Z and
`N` is a digit 0-9. Parts in brackets are optional.

If `validateinput=true` then the input string is checked for valid format using
a regex.
"""
ICD10GM(str::AbstractString, validateinput = false)

"""
    ICD10AM(icd::T) where {T<:AbstractICD10}

Convert ICD-10 code to Australian ICD-10-AM.

Note this just changes the type of the code. It does not do any checking or
translating of concepts between versions.
"""
ICD10AM(icd::T) where {T<:AbstractICD10}
"""
    ICD10CA(icd::T) where {T<:AbstractICD10}

Convert ICD-10 code to Canadian ICD-10-CA.

Note this just changes the type of the code. It does not do any checking or
translating of concepts between versions.
"""
ICD10CA(icd::T) where {T<:AbstractICD10}
# """
#     ICD10CM(icd::T) where {T<:AbstractICD10}

# Convert ICD-10 code to United States ICD-10-CM.

# Note this just changes the type of the code. It does not do any checking or
# translating of concepts between versions.
# """
# ICD10CM(icd::T) where {T<:AbstractICD10}
"""
    ICD10GM(icd::T) where {T<:AbstractICD10}

Convert ICD-10 code to German ICD-10-GM.

Note this just changes the type of the code. It does not do any checking or
translating of concepts between versions.
"""
ICD10GM(icd::T) where {T<:AbstractICD10}


text:
module Print_config
    import ..LTK_universe:Universe
    function write_config(filename::String,univ::Universe)
        fp = open(filename,"w")
        NT = univ.NT
        NZ = univ.NZ
        NY = univ.NY
        NX = univ.NX
        NC = univ.NC

        for it=1:NT
            for iz=1:NZ
                for iy=1:NY
                    for ix=1:NX
                        for μ=1:4
                            for a=1:NC
                                for b=1:NC
                                    #[:,:,ix,iy,iz,it]
                                    write(fp,real(univ.U[μ][a,b,ix,iy,iz,it]))
                                    write(fp,imag(univ.U[μ][a,b,ix,iy,iz,it]))
                                end
                            end

                        end
                    end
                end
            end
        end
        close(fp)

    end
end


text:
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the compressible ideal GLM-MHD equations

equations = IdealGlmMhdEquations3D(1.4)

"""
    initial_condition_blast_wave(x, t, equations::IdealGlmMhdEquations3D)

Weak magnetic blast wave setup taken from Section 6.1 of the paper:
- A. M. Rueda-Ramírez, S. Hennemann, F. J. Hindenlang, A. R. Winters, G. J. Gassner (2021)
  An entropy stable nodal discontinuous Galerkin method for the resistive MHD
  equations. Part II: Subcell finite volume shock capturing
  [doi: 10.1016/j.jcp.2021.110580](https://doi.org/10.1016/j.jcp.2021.110580)
"""
function initial_condition_blast_wave(x, t, equations::IdealGlmMhdEquations3D)
  # Center of the blast wave is selected for the domain [0, 3]^3
  inicenter = (1.5, 1.5, 1.5)
  x_norm = x[1] - inicenter[1]
  y_norm = x[2] - inicenter[2]
  z_norm = x[3] - inicenter[3]
  r = sqrt(x_norm^2 + y_norm^2 + z_norm^2)

  delta_0 = 0.1
  r_0 = 0.3
  lambda = exp(5.0 / delta_0 * (r - r_0))

  prim_inner = SVector(1.2, 0.1,  0.0, 0.1, 0.9, 1.0, 1.0, 1.0, 0.0)
  prim_outer = SVector(1.2, 0.2, -0.4, 0.2, 0.3, 1.0, 1.0, 1.0, 0.0)
  prim_vars = (prim_inner + lambda * prim_outer) / (1.0 + lambda)

  return prim2cons(prim_vars, equations)
end
initial_condition = initial_condition_blast_wave

surface_flux = (flux_lax_friedrichs, flux_nonconservative_powell)
volume_flux  = (flux_hindenlang_gassner, flux_nonconservative_powell)
polydeg = 3
basis = LobattoLegendreBasis(polydeg)
indicator_sc = IndicatorHennemannGassner(equations, basis,
                                         alpha_max=0.5,
                                         alpha_min=0.001,
                                         alpha_smooth=true,
                                         variable=density_pressure)
volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
                                                 volume_flux_dg=volume_flux,
                                                 volume_flux_fv=surface_flux)

solver = DGSEM(polydeg=polydeg, surface_flux=surface_flux, volume_integral=volume_integral)

# Mapping as described in https://arxiv.org/abs/2012.12040 but with slightly less warping.
# The mapping will be interpolated at tree level, and then refined without changing
# the geometry interpolant.
function mapping(xi_, eta_, zeta_)
  # Transform input variables between -1 and 1 onto [0,3]
  xi = 1.5 * xi_ + 1.5
  eta = 1.5 * eta_ + 1.5
  zeta = 1.5 * zeta_ + 1.5

  y = eta + 3/11 * (cos(1.5 * pi * (2 * xi - 3)/3) *
                    cos(0.5 * pi * (2 * eta - 3)/3) *
                    cos(0.5 * pi * (2 * zeta - 3)/3))

  x = xi + 3/11 * (cos(0.5 * pi * (2 * xi - 3)/3) *
                   cos(2 * pi * (2 * y - 3)/3) *
                   cos(0.5 * pi * (2 * zeta - 3)/3))

  z = zeta + 3/11 * (cos(0.5 * pi * (2 * x - 3)/3) *
                     cos(pi * (2 * y - 3)/3) *
                     cos(0.5 * pi * (2 * zeta - 3)/3))

  return SVector(x, y, z)
end

trees_per_dimension = (2, 2, 2)
mesh = P4estMesh(trees_per_dimension,
                 polydeg=3,
                 mapping=mapping,
                 initial_refinement_level=2,
                 periodicity=true)

# create the semi discretization object
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)

###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 0.5)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval)

alive_callback = AliveCallback(analysis_interval=analysis_interval)

amr_indicator = IndicatorLöhner(semi,
                                variable=density_pressure)
amr_controller = ControllerThreeLevel(semi, amr_indicator,
                                      base_level=2,
                                      max_level =4, max_threshold=0.15)
amr_callback = AMRCallback(semi, amr_controller,
                           interval=5,
                           adapt_initial_condition=true,
                           adapt_initial_condition_only_refine=true)

cfl = 1.4
stepsize_callback = StepsizeCallback(cfl=cfl)

glm_speed_callback = GlmSpeedCallback(glm_scale=0.5, cfl=cfl)

callbacks = CallbackSet(summary_callback,
                        analysis_callback,
                        alive_callback,
                        amr_callback,
                        stepsize_callback,
                        glm_speed_callback)

###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary


text:
using Gridap

u_2d(x) = x[1]+x[2]
u_3d(x) = x[1]+x[2]+x[3]
domain_2d = (0,1,0,1)
domain_3d = (0,1,0,1,0,1)
cells_2d = (4,4)
cells_3d = (4,4,4)

for cells in (cells_2d,cells_3d)
  domain = length(cells) == 3 ? domain_3d : domain_2d
  u = length(cells) == 3 ? u_3d : u_2d
  model = CartesianDiscreteModel(domain,cells)
  order = 1
  reffe = ReferenceFE(lagrangian,Float64,order)
  V = TestFESpace(model,reffe,dirichlet_tags="boundary")
  U = TrialFESpace(u,V)
  Ω = Triangulation(model)
  dΩ = Measure(Ω,2*order)
  a(u,v) = ∫( ∇(u)⋅∇(v) )dΩ
  l(v) = 0
  op = AffineFEOperator(a,l,U,V)
  uh = solve(op)
  sum(∫( abs2(u-uh) )dΩ)
end


text:
#md# Test the performance of different trainers
using WordCloud
using Random

println("This test will take several minutes")
@show Threads.nthreads()
words = [Random.randstring(rand(1:8)) for i in 1:200]
weights = randexp(length(words)) .* 2000 .+ rand(20:100, length(words));
wc1 = wordcloud(words, weights, mask=shape(ellipse, 500, 500, color=0.15), angles=(0,90,45), density=0.55)

words = [Random.randstring(rand(1:8)) for i in 1:400]
weights = randexp(length(words)) .* 2000 .+ rand(20:100, length(words));
wc2 = wordcloud(words, weights, mask=shape(ellipse, 500, 500, color=0.15), angles=(0,90,45))

words = [Random.randstring(rand(1:8)) for i in 1:5000]
weights = randexp(length(words)) .* 2000 .+ rand(20:100, length(words));
wc3 = wordcloud(words, weights, mask=shape(box, 2000, 2000, 100, color=0.15), angles=(0,90,45))

wcs = [wc1, wc1, wc2, wc3] #repeat wc1 to trigger compiling
ts = [WordCloud.Stuffing.trainepoch_E!,WordCloud.Stuffing.trainepoch_EM!,
WordCloud.Stuffing.trainepoch_EM2!,WordCloud.Stuffing.trainepoch_EM3!,
WordCloud.Stuffing.trainepoch_P!,WordCloud.Stuffing.trainepoch_P2!,WordCloud.Stuffing.trainepoch_Px!]
es = [[] for i in 1:length(wcs)]
for (i,wc) in enumerate(wcs)
    println("\n\n", "*"^10, "wordcloud - $(length(wc.words)) words on mask$(size(wc.mask))", "*"^10)
    for (j,t) in enumerate(ts)
        println("\n", i-1, "==== ", j, "/", length(ts), " ", nameof(t))
        placewords!(wc)
        @time e = @elapsed generate!(wc, trainer=t, retry=1)
        push!(es[i], string(nameof(t)) * (getstate(wc)==:generate! ? " ✔" : " ✘")=>e)
    end
end
println("SUMMARY")
for (i,(wc,e)) in enumerate(zip(wcs, es))
    println("##$(i-1) $(length(wc.words))@$(size(wc.mask)):")
    println(repr("text/plain", e))
end


text:
using TransformVariables,
      LogDensityProblems,
      DynamicHMC,
      Distributions,
      Statistics,
      StatsFuns,
      ForwardDiff
import LogDensityProblems: ADgradient

using Random

export dynamicHMC

function dynamicHMC(
    rng::AbstractRNG,
    m::JointDistribution,
    _data,
    N::Int = 1000;
    method = logpdf,
    ad_backend = Val(:ForwardDiff),
    reporter = DynamicHMC.NoProgressReport(),
    kwargs...,
)
    ℓ(pars) = logpdf(m, merge(pars, _data), method)
    t = xform(m, _data)
    P = LogDensityProblems.TransformedLogDensity(t, ℓ)
    ∇P = LogDensityProblems.ADgradient(ad_backend, P)

    results = DynamicHMC.mcmc_with_warmup(
        rng,
        ∇P,
        N;
        reporter = reporter,
        kwargs...,
    )
    samples = TransformVariables.transform.(t, results.chain)
    return samples
end

function dynamicHMC(
    rng::AbstractRNG,
    m::JointDistribution,
    _data,
    ::Val{Inf};
    method = logpdf,
    ad_backend = Val(:ForwardDiff),
    reporter = DynamicHMC.NoProgressReport(),
    kwargs...,
)
    ℓ(pars) = logpdf(m, merge(pars, _data), method)
    t = xform(m, _data)
    P = LogDensityProblems.TransformedLogDensity(t, ℓ)
    ∇P = LogDensityProblems.ADgradient(ad_backend, P)

    results = DynamicHMC.mcmc_keep_warmup(
        rng,
        ∇P,
        0;
        reporter = reporter,
        kwargs...,
    )
    steps = DynamicHMC.mcmc_steps(
        results.sampling_logdensity,
        results.final_warmup_state,
    )
    return results, steps
end

function dynamicHMC(m::JointDistribution, args...; kwargs...)
    return dynamicHMC(Random.GLOBAL_RNG, m, args...; kwargs...)
end

using ResumableFunctions

export stream

@resumable function stream(
    rng::AbstractRNG,
    f::typeof(dynamicHMC),
    m::JointDistribution,
    _data::NamedTuple,
)
    t = xform(m, _data)
    (results, steps) = dynamicHMC(rng, m, _data, Val(Inf))
    Q = results.final_warmup_state.Q
    while true
        Q, tree_stats = DynamicHMC.mcmc_next_step(steps, Q)
        @yield (merge(t(Q.q), (_ℓ = Q.ℓq,)), tree_stats)
    end
end

function stream(
    f::typeof(dynamicHMC),
    m::JointDistribution,
    _data::NamedTuple;
    kwargs...,
)
    return stream(Random.GLOBAL_RNG, f, m, _data; kwargs...)
end


text:
using JFVM
import PyPlot, JSON

include("../functions/rel_perms_real.jl")

function polymer_viscosity(cs, cp)
    a1 = 4.0
    a2 = 0.0
    a3 = 6.0
    sp = -0.24
    μ_w = 0.001 # Pa.s
    return μ_w*(1.0+cp.*cs.^sp.*(a1+a2*cp+a3*cp.*cp))
end

function dpolymer_viscosity_dcs(cs, cp)
    a1 = 4.0
    a2 = 0.0
    a3 = 6.0
    sp = -0.24
    μ_w = 0.001 # Pa.s
    return sp*μ_w*cp.*cs.^(sp-1).*(a1+a2*cp+a3*cp.*cp)
end

function dpolymer_viscosity_dcp(cs, cp)
    a1 = 4.0
    a2 = 0.0
    a3 = 6.0
    sp = -0.24
    μ_w = 0.001 # Pa.s
    return μ_w*cs.^sp.*(a1+2*a2*cp+3*a3*cp.*cp)
end

# domain
Lx = 1.0
Ly = 1.0
Nx = 300
Ny = 10
m = createMesh1D(Nx, Lx)

# (petro)physical properties
mu_oil = 10e-3 # Pa.s
perm = 0.1e-12
poros = 0.3

# rel-perms
krw0 = 0.2
kro0 = 0.8
n_o  = 2.0
n_w  = 2.0
swc  = 0.08
sor = 0.2

KRW  = sw -> krw.(sw, krw0, sor, swc, n_w)
dKRW = sw -> dkrwdsw.(sw, krw0, sor, swc, n_w)
KRO  = sw -> kro.(sw, kro0, sor, swc, n_o)
dKRO = sw -> dkrodsw.(sw, kro0, sor, swc, n_o)

# initial and boundary conditions
p_res = 100e5 # Pa reservoir (and right boundary pressure)
sw_0 = swc # inital water saturation
sw_inj = 1.0 # injection water saturation
cp_0 = 0.0 # initial polymer concentration kg/m^3
cp_inj = 1.0 # injected polymer concentration
cs_0 = 10.0  # initial salt concentration kg/m^3
cs_inj = 5.0 # injected salt concentration

φ = createCellVariable(m, poros) # porosity
k = createCellVariable(m, perm) # m^2 permeability
k_face   = harmonicMean(k) # harmonic averaging for perm
u_inj = 1e-5 # m/s injection Darcy velocity
pv_inj = 0.3 # injected pore volumes
t_final = pv_inj*Lx/(u_inj/poros) # s final time
dt0 = t_final/Nx # s time step
dt  = t_final/Nx # variable time step

BCp = createBC(m)   # pressure boundary
BCp.left.a[:] = perm/polymer_viscosity(cs_inj, cp_inj)
BCp.left.b[:] = 0.0
BCp.left.c[:] = -u_inj
BCp.right.a[:] = 0.0
BCp.right.b[:] = 1.0
BCp.right.c[:] = p_res
M_bc_p, RHS_bc_p = boundaryConditionTerm(BCp)

# saturation boundary
BCs = createBC(m)
BCs.left.a[:] = 0.0
BCs.left.b[:] = 1.0
BCs.left.c[:] = sw_inj # injected water saturation
M_bc_s, RHS_bc_s = boundaryConditionTerm(BCs)

# salt concentration boundary
BCcs = createBC(m)
BCcs.left.a[:] = 0.0
BCcs.left.b[:] = 1.0
BCcs.left.c[:] = cs_inj # injected salt concentration
M_bc_cs, RHS_bc_cs = boundaryConditionTerm(BCcs)

# polymer concentration boundary
BCcp = createBC(m)
BCcp.left.a[:] = 0.0
BCcp.left.b[:] = 1.0
BCcp.left.c[:] = cp_inj # injected polymer concentration
M_bc_cp, RHS_bc_cp = boundaryConditionTerm(BCcp)

# initial conditions
p_init  = createCellVariable(m, p_res)
sw_init = createCellVariable(m, sw_0)
cs_init = createCellVariable(m, cs_0)
cp_init = createCellVariable(m, cp_0)
p_val  = createCellVariable(m, p_res)
sw_val = createCellVariable(m, sw_0)
cs_val = createCellVariable(m, cs_0)
cp_val = createCellVariable(m, cp_0)

uw       = -gradientTerm(p_val)
# M_bc_p, RHS_bc_p = boundaryConditionTerm(BCp)
# M_bc_s, RHS_bc_s = boundaryConditionTerm(BCs)
# M_bc_c, RHS_bc_c = boundaryConditionTerm(BCcs)
# M_bc_t, RHS_bc_t = boundaryConditionTerm(BCt)
tol_s = 1e-7
tol_c = 1e-7
max_int_loop = 40
t = 0.0
FL = fluxLimiter("SUPERBEE") # flux limiter

while t<t_final
    error_s = 2*tol_s
    error_cp = 2*tol_c
    loop_countor = 0
    while (error_s>tol_s) || (error_cp>tol_c)
        loop_countor += 1
        if loop_countor > max_int_loop
            sw_val = copyCell(sw_init)
            cp_val  = copyCell(cp_init)
            cs_val  = copyCell(cs_init)
            p_val  = copyCell(p_init)
            dt = dt/3.0
            break
        end
        ∇p0      = gradientTerm(p_val)
        ∇s0      = gradientTerm(sw_val)

        cs_face   = upwindMean(cs_val, uw)
        cp_face   = upwindMean(cp_val, uw)
        sw_face       = upwindMean(sw_val, uw)

        mu_water_face = faceEval(polymer_viscosity, cs_face, cp_face)
        d_mu_cs_face = faceEval(dpolymer_viscosity_dcs, cs_face, cp_face)
        d_mu_cp_face = faceEval(dpolymer_viscosity_dcp, cs_face, cp_face)
        d_cs_mu_dcs = (mu_water_face-d_mu_cs_face.*cs_face)./(mu_water_face.*mu_water_face)
        d_cs_mu_dcp = (-d_mu_cp_face.*cs_face)./(mu_water_face.*mu_water_face)
        d_cp_mu_dcp = (mu_water_face-d_mu_cp_face.*cp_face)./(mu_water_face.*mu_water_face)
        d_cp_mu_dcs = (-d_mu_cs_face.*cp_face)./(mu_water_face.*mu_water_face)
        d_mu_dcs = -d_mu_cs_face./(mu_water_face.*mu_water_face)
        d_mu_dcp = -d_mu_cp_face./(mu_water_face.*mu_water_face)

        # println("visc")
        krw_face      = faceEval(KRW,sw_face)
        kro_face      = faceEval(KRO,sw_face)

        λ_w_face      = k_face.*krw_face./mu_water_face
        λ_o_face      = k_face.*kro_face/mu_oil
        uw            = -λ_w_face.*∇p0
        uo            = -λ_o_face.*∇p0
        ut            = uw + uo

        dλ_w_face     = k_face.*faceEval(dKRW,sw_face)./mu_water_face
        dλ_o_face     = k_face.*faceEval(dKRO,sw_face)./mu_oil

        # water mass balance (wmb)
        M_t_s_wmb, RHS_t_s_wmb = transientTerm(sw_init, dt, φ)
        M_d_p_wmb = diffusionTerm(-λ_w_face)
        # M_a_s_wmb = convectionUpwindTerm(-dλ_w_face.*∇p0, ut)
        M_a_s_wmb, RHS_a_s_wmb = convectionTvdTerm(-dλ_w_face.*∇p0, sw_val, FL, ut)
        # M_a_cs_wmb = convectionUpwindTerm(-k_face.*krw_face.*d_mu_dcs.*∇p0, ut)
        M_a_cs_wmb, RHS_a_cs_wmb = convectionTvdTerm(-k_face.*krw_face.*d_mu_dcs.*∇p0, cs_val, FL, ut)
        # M_a_cp_wmb = convectionUpwindTerm(-k_face.*krw_face.*d_mu_dcp.*∇p0, ut)
        M_a_cp_wmb, RHS_a_cp_wmb = convectionTvdTerm(-k_face.*krw_face.*d_mu_dcp.*∇p0, cp_val, FL, ut)

        # oil mass balance (omb)
        M_t_s_omb, RHS_t_s_omb = transientTerm(sw_init, dt, -φ)
        M_d_p_omb = diffusionTerm(-λ_o_face)
        # M_a_s_omb = convectionUpwindTerm(-(dλ_o_face.*∇p0), ut)
        M_a_s_omb, RHS_a_s_omb = convectionTvdTerm(-(dλ_o_face.*∇p0), sw_val, FL, ut)

        # salt mass balance
        M_t_cs_smb, RHS_t_cs_smb = transientTerm(cs_init, dt, sw_val.*φ)
        M_t_s_smb, RHS_t_s_smb = transientTerm(sw_init, dt, cs_val.*φ)
        # M_a_cs_smb, RHS_a_cs_smb = convectionTvdTerm(uw, cs_val, FL, ut)
        # M_a_s_smb = convectionUpwindTerm(-dλ_w_face.*cs_face.*∇p0, ut)
        M_a_s_smb, RHS_a_s_smb = convectionTvdTerm(-dλ_w_face.*cs_face.*∇p0, sw_val, FL, ut)
        M_d_p_smb = diffusionTerm(-λ_w_face.*cs_face)
        # M_a_cs_smb = convectionUpwindTerm(-k_face.*krw_face.*d_cs_mu_dcs.*∇p0, ut)
        # M_a_cp_smb = convectionUpwindTerm(-k_face.*krw_face.*d_cs_mu_dcp.*∇p0, ut)
        M_a_cs_smb, RHS_a_cs_smb = convectionTvdTerm(-k_face.*krw_face.*d_cs_mu_dcs.*∇p0, cs_val, FL, ut)
        M_a_cp_smb, RHS_a_cp_smb = convectionTvdTerm(-k_face.*krw_face.*d_cs_mu_dcp.*∇p0, cp_val, FL, ut)

        # polymer mass balance
        M_t_cp_pmb, RHS_t_cp_pmb = transientTerm(cp_init, dt, sw_val.*φ)
        M_t_s_pmb, RHS_t_s_pmb = transientTerm(sw_init, dt, cp_val.*φ)
        # M_a_s_pmb = convectionUpwindTerm(-dλ_w_face.*cp_face.*∇p0, ut)
        M_a_s_pmb, RHS_a_s_pmb = convectionTvdTerm(-dλ_w_face.*cp_face.*∇p0, sw_val, FL, ut)
        # M_a_cp, RHS_a_cp = convectionTvdTerm(uw, cp_val, FL, ut)
        M_d_p_pmb = diffusionTerm(-λ_w_face.*cp_face)
        # M_a_cs_pmb = convectionUpwindTerm(-k_face.*krw_face.*d_cp_mu_dcs.*∇p0, ut)
        # M_a_cp_pmb = convectionUpwindTerm(-k_face.*krw_face.*d_cp_mu_dcp.*∇p0, ut)
        M_a_cs_pmb, RHS_a_cs_pmb = convectionTvdTerm(-k_face.*krw_face.*d_cp_mu_dcs.*∇p0, cs_val, FL, ut)
        M_a_cp_pmb, RHS_a_cp_pmb = convectionTvdTerm(-k_face.*krw_face.*d_cp_mu_dcp.*∇p0, cp_val, FL, ut)

        # create the PDE system M [p;s;c]=RHS
        # x_val = [p_val.value[:]; sw_val.value[:]; c_val.value[:]]
        M = [M_bc_p+M_d_p_wmb   M_t_s_wmb+M_a_s_wmb  M_a_cs_wmb  M_a_cp_wmb;
            M_d_p_omb   M_bc_s+M_t_s_omb+M_a_s_omb  zeros(Nx+2,Nx+2)  zeros(Nx+2,Nx+2);
            M_d_p_smb  M_a_s_smb+M_t_s_smb  M_a_cs_smb+M_t_cs_smb+M_bc_cs  M_a_cp_smb;
            M_d_p_pmb  M_a_s_pmb+M_t_s_pmb  M_a_cs_pmb  M_t_cp_pmb+M_bc_cp+M_a_cp_pmb]

        RHS = [RHS_bc_p+RHS_t_s_wmb+RHS_a_s_wmb+RHS_a_cs_wmb+RHS_a_cp_wmb+M_a_s_wmb*sw_val.value[:]+M_a_cs_wmb*cs_val.value[:]+M_a_cp_wmb*cp_val.value[:];
               RHS_bc_s+RHS_t_s_omb+RHS_a_s_omb+(M_a_s_omb)*sw_val.value[:];
               RHS_t_cs_smb+RHS_t_s_smb+RHS_bc_cs+RHS_a_s_smb+RHS_a_cs_smb+RHS_a_cp_smb+M_a_s_smb*sw_val.value[:]+M_a_cs_smb*cs_val.value[:]+M_a_cp_smb*cp_val.value[:];
               RHS_t_cp_pmb+RHS_t_s_pmb+RHS_bc_cp+RHS_a_s_pmb+RHS_a_cs_pmb+RHS_a_cp_pmb+M_a_s_pmb*sw_val.value[:]+M_a_cs_pmb*cs_val.value[:]+M_a_cp_pmb*cp_val.value[:]]

        # x_sol = solveLinearPDE(m, M, RHS)
        x_sol = 0
        try
            x_sol = M\RHS
        catch
            loop_countor = max_int_loop
            continue
        end

        p_new = x_sol[1:Nx+2]
        s_new = x_sol[(1:Nx+2)+(Nx+2)]
        cs_new = x_sol[(1:Nx+2)+2*(Nx+2)]
        cp_new = x_sol[(1:Nx+2)+3*(Nx+2)]

        error_s = sum(abs, s_new[2:end-1]-internalCells(sw_val))
        error_cp = sum(abs, cp_new[2:end-1]-internalCells(cp_val))

        #   error_c = sum(abs, internalCells(cp_val)-internalCells(cp_val))
        # println(error_s)
        # println(error_c)
        # w_p = 0.9
        # w_sw = 0.8
        # w_c = 0.8
        # p_val.value[:] = w_p*p_new[:]+(1-w_p)*p_val.value[:]
        # sw_val.value[:] = w_sw*s_new[:]+(1-w_sw)*sw_val.value[:]
        # c_val.value[:] = w_c*c_new[:]+(1-w_c)*c_val.value[:]
        p_val.value[:] = p_new[:]
        dsw_apple = 0.2
        eps_apple = sqrt(eps()) # 1e-5
        for i in eachindex(s_new)
            if s_new[i]>=(1-sor)
                if sw_val.value[i]<(1-sor-eps_apple)
                    sw_val.value[i] = 1-sor-eps_apple
                else
                    sw_val.value[i] = 1-sor
                end
            elseif s_new[i]<=swc
                if sw_val.value[i]> swc+eps_apple
                    sw_val.value[i] = swc+eps_apple
                else
                    sw_val.value[i] = swc
                end
            elseif abs(s_new[i]-sw_val.value[i])>dsw_apple
                sw_val.value[i] += dsw_apple*sign(s_new[i]-sw_val.value[i])
            else
                sw_val.value[i] = s_new[i]
            end
        end
        sw_val = createCellVariable(m, internalCells(sw_val), BCs)
        cs_new[cs_new.<0.0] = 0.0
        cp_new[cp_new.<0.0] = 0.0

        cs_val.value[:] = cs_new[:]
        cp_val.value[:] = cp_new[:]

        cs_val = createCellVariable(m, cs_new[2:end-1], BCcs)
        cp_val = createCellVariable(m, cp_new[2:end-1], BCcp)

        # GR.plot(sw_val.value)
        # sw_val.value[:] = w_sw*s_new[:]+(1-w_sw)*sw_val.value[:]
    end
    if loop_countor<max_int_loop
        p_init = copyCell(p_val)
        sw_init = copyCell(sw_val)
        cp_init = copyCell(cp_val)
        cs_init = copyCell(cs_val)

        t +=dt
        dt = dt0
        println("progress: $(t/(t_final)*100) [%]")
    end
end
# JLD.save("results2D/all_data_heterogen_wf.jld", "p", p_val, "c", c_val,
# "c_oil", c_oil, "sw", sw_val, "c_tracer", c_t_val)
# figure(1)
figure()
visualizeCells(sw_init)
# title("Sw")
visualizeCells(cs_init)
# title("c_DME")
visualizeCells(cp_init)
# title("c_tracer")
# savefig("profile.png")
# figure(2)
# plot(t_s, rec_fact)
# savefig("recovery.png")
# figure(3)
# plot(t_s, dp_hist)
# savefig("dp.png")
# PyPlot.pcolormesh(sw_val.value[2:end-1, 2:end-1]', cmap="YlGnBu", shading = "gouraud", vmin=0, vmax=1)

max_stars_repo_name: simulkade/peteng

max_issues_repo_name: simulkade/peteng

max_forks_repo_name: simulkade/peteng


text:
function unit_test()
    test_helpers()
    test_filters()
    Logging.info( "Manager Unit Tests passing. Good job.")
end

max_stars_repo_name: shdann/autoshift

max_issues_repo_name: shdann/autoshift

max_forks_repo_name: shdann/autoshift


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

"""
Deprecated.
Converts CAS csvs to a single dataframe.
"""
module CSVFeatures

export feature_frame, csv_to_dataframe

import Compat.ASCIIString
using RLESUtils, LookupCallbacks, FileUtils
using DataFrames
using Iterators

function feature_frame(csvfile::AbstractString, feature_map::Vector{LookupCallback},
                        colnames::Vector{Symbol}=Symbol[])
  csv = readcsv(csvfile)
  headers = csv[1, :] |> vec
  units = csv[2, :] |> vec
  dat = csv[3:end, :]
  tmax = size(dat, 1)
  V = Array(Any, length(feature_map))

  lookup_ids = map(feature_map) do lcb
    ids = map(lcb.lookups) do l
      findfirst(x -> x == l, headers)
    end
    any(x -> x == 0, ids) && error("Lookup not found: $lcb") #all lookup ids should be found
    return ids
  end
  for i = 1:length(V)
    input_vars = dat[:, lookup_ids[i]]
    f = feature_map[i].callback
    V[i] = mapslices(x -> f(x...), input_vars, 2) |> vec
  end
  return DataFrame(V, colnames)
end

append{T<:AbstractString}(V::Vector{T}, s::AbstractString) = map(x -> "$(x)$(s)", V)

#map a directory of trajSave_aircraft.csv files to trajSave_dataframe.csv file
function csv_to_dataframe{T<:AbstractString}(files::Vector{T}, feature_map::Vector{LookupCallback}, feature_names::Vector{T}; outdir::AbstractString="./")

  @assert length(feature_map) == length(feature_names)
  sort!(files)
  grouped_files = Iterators.groupby(x -> split(x, "_aircraft")[1], files) |> collect
  outfiles = Array(ASCIIString, length(grouped_files))
  for i = 1:length(grouped_files)
    D_ = map(enumerate(grouped_files[i])) do jf
      (j, f) = jf
      feature_frame(f, feature_map, convert(Vector{Symbol}, append(feature_names, "_$j")))
    end
    D = hcat(D_...)
    df_file = split(grouped_files[i][1], "_aircraft")[1] * "_dataframe.csv.gz"
    df_file = joinpath(outdir, basename(df_file))
    outfiles[i] = df_file
    writetable(df_file, D)
  end
  return outfiles
end

end #module


text:
using Resonance

methylation = CSV.read("data/targets_passed.csv", DataFrame)
wrangled = CSV.read("data/wrangled.csv", DataFrame)

jnd = leftjoin(methylation, wrangled, on=["Sample_ID"=>"subject", "Timepoint"=>"timepoint"], matchmissing=:notequal)

jnd.Sample_Name = jnd.sample

files = filter(f-> any(s-> !ismissing(s) && contains(f, s), jnd.sample), readdir("/augusta/echo/analysis/biobakery3/links/metaphlan/"))
open("data/methylation_files.txt", "w") do io
    println.(Ref(io), files)
end

met = metaphlan_profiles(joinpath.(Ref("/augusta/echo/analysis/biobakery3/links/metaphlan/"), files))
CSV.write("/home/kevin/Desktop/methylation_taxa.csv", met)
CSV.write("/home/kevin/Desktop/methylation_taxa_link.csv", select(jnd, ["Sample_ID", "Timepoint", "Sample_Name"]))

max_stars_repo_name: Klepac-Ceraj-Lab/ResonanceAnalysis

max_issues_repo_name: Klepac-Ceraj-Lab/ResonanceAnalysis

max_forks_repo_name: Klepac-Ceraj-Lab/ResonanceAnalysis