text:
module WignerSeitzCell
   export wigner_seitz_cell
   #
   using LinearAlgebra: norm
   using StaticArrays

   const ws_search_range = 3

   # for performance, please use SMatrix{3,3} for at.
   #return function to get length of a vector
   length_func(at) = vec -> norm( at * vec )

   function wigner_seitz_cell(rdim, lattice::AbstractMatrix, tau::AbstractMatrix)
      get_length = length_func( SMatrix{3,3}(lattice) )
      cutoff = set_cutoff(get_length, rdim)

      ntmp, ntau = size(tau)
      ntmp != 3 && throw(DimensionMismatch("tau dimensions are $(size(tau))"))

      delta_tau = [SVector{3}(view(tau,:,j) .- view(tau,:,i)) for j in 1:ntau for i in 1:j]

      ws_rvecs = Vector{ SVector{3, Float64} }()
      rvec_idx_tmp = Vector{ Vector{Int} }(undef, length(delta_tau))
      rvec_idx = Vector{ Vector{Int} }(undef, length(delta_tau))

      (nr1, nr2, nr3) = Tuple(rdim)
      for rv in Iterators.product(0:nr1-1, 0:nr2-1, 0:nr3-1)
         vecs = valid_cells(get_length, rv, rdim, cutoff)

         vecs_label = zeros(Int, length(vecs))
         for (it, dt) in enumerate(delta_tau)
            dist = [get_length( v + dt) for v in vecs]
            mindist = minimum(dist)
            idx = findall(x -> ((x-mindist) < 1.0E-6), dist)
            #
            vecs_label[idx] .+= 1
            rvec_idx_tmp[it] = idx
         end

         #collect vectors
         col_idx = findall(x->(x > 0), vecs_label)
         #re-label
         vecs_label[col_idx] = 1:length(col_idx)
         #
         prev_len = length(ws_rvecs)
         append!(ws_rvecs, vecs[col_idx])

         #update rvec_idx
         for it in eachindex(delta_tau)
            rvec_idx_tmp[it] .= prev_len .+ vecs_label[ rvec_idx_tmp[it] ]
            try
               append!(rvec_idx[it], rvec_idx_tmp[it])
            catch error
               error isa UndefRefError ? rvec_idx[it]=rvec_idx_tmp[it] : rethrow(error)
            end
         end
      end
      return ws_rvecs, rvec_idx
   end

   function set_cutoff(get_length::Function, rdim; lcutoff = false)
      ndim = lcutoff ? rdim : (rdim .÷ 2 .+ 1)

      max_len = 0.0
      for vec in Iterators.product( (-1, 1), (-1, 1), (-1, 1) )
         dist = get_length( SVector{3}(float.(vec .* ndim)) )
         dist > max_len && (max_len = dist)
      end
      return max_len
   end

   function valid_cells(get_length::Function, v0, rdim, cutoff)
      vectors = Vector{SVector{3,Float64}}()
      nr = ws_search_range
      for i in Iterators.product(-nr:nr, -nr:nr, -nr:nr)
         vec = SVector{3}(float.(i .* rdim .+ v0))
         get_length(vec) < cutoff && push!(vectors, vec)
      end
      return vectors
   end
end #module
#   #lattice:  'at' in QE
#   function wigner_seitz_cell(ndim, lattice, tau_a, tau_b)
#      #
#      get_length = length_func( SMatrix{3,3}(lattice) )
#      cutoff = set_cutoff(get_length, ndim)
#
#      r_vectors = rvec_image(get_length, ndim, cutoff)
#
#      wscells = Vector{ SVector{3, Float64} }()
#      #
#      for vecs in r_vectors
#         #
#         dist = zeros( length(vecs) )
#         for (i, v) in enumerate(vecs)
#            dist[i] = get_length( SVector{3}(float.(v .+ tau_b .- tau_a)) )
#         end
#         #
#         mindist = minimum(dist)
#         idx = findall(x -> ( (x-mindist) < 1.0E-6), dist)
#         #
#         append!(wscells, vecs[idx])
#      end
#      #
#      return wscells
#   end #wiger_seitz_cell

#   let rvecs = Vector{Vector{ SVector{3,Float64} }}(), ndim = (0, 0, 0)
#      #
#      global function rvec_image(get_length::Function, rdim, cutoff)
#         Tuple(rdim) != ndim && rvec_image!(get_length, rdim, cutoff, rvecs)
#         return copy(rvecs)
#      end
#      #
#      function rvec_image!(get_length::Function, rdim, cutoff, rvecs)
#         empty!(rvecs)
#         (nr1, nr2, nr3) = ndim = Tuple(rdim)
#
#         for i in Iterators.product(0:nr1-1, 0:nr2-1, 0:nr3-1)
#            v0 = SVector(float.(i))
#            push!(rvecs, valid_cells(get_length, v0, rdim, cutoff))
#         end
#         return nothing
#      end
#   end #let block


text:
using BCVUMPS
using BCVUMPS: model_tensor
using Test
using OMEinsum

@testset "exampletensor" begin
    β = rand()
    M = model_tensor(Ising(),β)
    @test size(M) == (1,1)
    for Nj = 1:3, Ni = 1:3
        M = model_tensor(Ising(Ni,Nj),β)
        @test size(M) == (Ni,Nj)
        @test M[1,1] ≈ M[Ni,Nj]
    end

    @test Ising22(1).Ni == 2
    @test Ising22(1).Nj == 2
    M = model_tensor(Ising22(1),β)
    @test M[1,1] ≈ M[1,2] ≈ M[2,1] ≈ M[2,2]

    M = model_tensor(Ising22(2),β)
    @test M[1,2] ≈ permutedims(M[1,1],[3,2,1,4])
    @test M[2,1] ≈ permutedims(M[1,1],[1,4,3,2])
    @test M[2,2] ≈ permutedims(M[1,1],[3,4,1,2])
end

ext: jl


text:
module LineEdit

using ..Terminals

import ..Terminals: raw!, width, height, cmove, getX,
                       getY, clear_line, beep

import Base: ensureroom, peek, show

abstract TextInterface

export run_interface, Prompt, ModalInterface, transition, reset_state, edit_insert

immutable ModalInterface <: TextInterface
    modes
end

type MIState
    interface::ModalInterface
    current_mode
    aborted::Bool
    mode_state
    kill_buffer::ByteString
    previous_key::Array{Char,1}
    key_repeats::Int
end
MIState(i, c, a, m) = MIState(i, c, a, m, "", Char[], 0)

type Mode <: TextInterface
end

type Prompt <: TextInterface
    prompt
    first_prompt
    # A string or function to be printed before the prompt. May not change the length of the prompt.
    # This may be used for changing the color, issuing other terminal escape codes, etc.
    prompt_prefix
    # Same as prefix except after the prompt
    prompt_suffix
    keymap_func
    keymap_func_data
    complete
    on_enter
    on_done
    hist
end

show(io::IO, x::Prompt) = show(io, string("Prompt(\"", x.prompt, "\",...)"))

immutable InputAreaState
    num_rows::Int64
    curs_row::Int64
end

type PromptState
    terminal::TextTerminal
    p::Prompt
    input_buffer::IOBuffer
    ias::InputAreaState
    indent::Int
end

input_string(s::PromptState) = bytestring(pointer(s.input_buffer.data), s.input_buffer.size)

input_string_newlines(s::PromptState) = count(c->(c == '\n'), input_string(s))
function input_string_newlines_aftercursor(s::PromptState)
    str = input_string(s)
    length(str) == 0 && return 0
    rest = str[nextind(str, position(s.input_buffer)):end]
    return count(c->(c == '\n'), rest)
end

abstract HistoryProvider
abstract CompletionProvider

type EmptyCompletionProvider <: CompletionProvider
end

type EmptyHistoryProvider <: HistoryProvider
end

reset_state(::EmptyHistoryProvider) = nothing

complete_line(c::EmptyCompletionProvider, s) = []

terminal(s::IO) = s
terminal(s::PromptState) = s.terminal

for f in [:terminal, :edit_insert, :on_enter, :add_history, :buffer, :edit_backspace, :(Base.isempty),
        :replace_line, :refresh_multi_line, :input_string, :edit_move_left, :edit_move_right,
        :edit_move_word_left, :edit_move_word_right, :update_display_buffer]
    @eval ($f)(s::MIState, args...) = $(f)(s.mode_state[s.current_mode], args...)
end

function common_prefix(completions)
    ret = ""
    c1 = completions[1]
    isempty(c1) && return ret
    i = 1
    cc, nexti = next(c1, i)
    while true
        for c in completions
            (i > endof(c) || c[i] != cc) && return ret
        end
        ret *= string(cc)
        i >= endof(c1) && return ret
        i = nexti
        cc, nexti = next(c1, i)
    end
end

# Show available completions
function show_completions(s::PromptState, completions)
    colmax = maximum(map(length, completions))
    num_cols = max(div(width(terminal(s)), colmax+2), 1)
    entries_per_col, r = divrem(length(completions), num_cols)
    entries_per_col += r != 0
    # skip any lines of input after the cursor
    cmove_down(terminal(s), input_string_newlines_aftercursor(s))
    println(terminal(s))
    for row = 1:entries_per_col
        for col = 0:num_cols
            idx = row + col*entries_per_col
            if idx <= length(completions)
                cmove_col(terminal(s), (colmax+2)*col)
                print(terminal(s), completions[idx])
            end
        end
        println(terminal(s))
    end
    # make space for the prompt
    for i = 1:input_string_newlines(s)
        println(terminal(s))
    end
end

# Prompt Completions
complete_line(s::MIState) = complete_line(s.mode_state[s.current_mode], s.key_repeats)
function complete_line(s::PromptState, repeats)
    completions, partial, should_complete = complete_line(s.p.complete, s)
    if length(completions) == 0
        beep(terminal(s))
    elseif !should_complete
        # should_complete is false for cases where we only want to show
        # a list of possible completions but not complete, e.g. foo(\t
        show_completions(s, completions)
    elseif length(completions) == 1
        # Replace word by completion
        prev_pos = position(s.input_buffer)
        seek(s.input_buffer, prev_pos-sizeof(partial))
        edit_replace(s, position(s.input_buffer), prev_pos, completions[1])
    else
        p = common_prefix(completions)
        if length(p) > 0 && p != partial
            # All possible completions share the same prefix, so we might as
            # well complete that
            prev_pos = position(s.input_buffer)
            seek(s.input_buffer, prev_pos-sizeof(partial))
            edit_replace(s, position(s.input_buffer), prev_pos, p)
        elseif repeats > 0
            show_completions(s, completions)
        end
    end
end

clear_input_area(terminal, s) = (_clear_input_area(terminal, s.ias); s.ias = InputAreaState(0, 0))
clear_input_area(s) = clear_input_area(s.terminal, s)
function _clear_input_area(terminal, state::InputAreaState)
    # Go to the last line
    if state.curs_row < state.num_rows
        cmove_down(terminal, state.num_rows-state.curs_row)
    end

    # Clear lines one by one going up
    for j=0:(state.num_rows-2)
        clear_line(terminal)
        cmove_up(terminal)
    end

    # Clear top line
    clear_line(terminal)
end

prompt_string(s::PromptState) = s.p.prompt
prompt_string(s::String) = s

refresh_multi_line(termbuf::TerminalBuffer, s::PromptState) = s.ias =
    refresh_multi_line(termbuf, terminal(s), buffer(s), s.ias, s, indent = s.indent)

function refresh_multi_line(termbuf::TerminalBuffer, terminal::UnixTerminal, buf, state::InputAreaState, prompt = ""; indent = 0)
    cols = width(terminal)

    _clear_input_area(termbuf, state)

    curs_row = -1 # relative to prompt
    curs_col = -1 # absolute
    curs_pos = -1 # 1-based column position of the cursor
    cur_row = 0
    buf_pos = position(buf)
    line_pos = buf_pos
    # Write out the prompt string
    write_prompt(termbuf, prompt)
    prompt = prompt_string(prompt)

    seek(buf, 0)

    llength = 0

    l = ""

    plength = strwidth(prompt)
    pslength = length(prompt.data)
    # Now go through the buffer line by line
    while cur_row == 0 || (!isempty(l) && l[end] == '\n')
        l = readline(buf)
        hasnl = !isempty(l) && l[end] == '\n'
        cur_row += 1
        # We need to deal with UTF8 characters. Since the IOBuffer is a bytearray, we just count bytes
        llength = strwidth(l)
        slength = length(l.data)
        if cur_row == 1 # First line
            if line_pos < slength
                num_chars = strwidth(l[1:line_pos])
                curs_row = div(plength+num_chars-1, cols) + 1
                curs_pos = (plength+num_chars-1) % cols + 1
            end
            # Substract -1 if there's a '\n' at the end of the line (since it doesn't take up a column)
            # The other -1, since we want 10,20 for cols=10 to still not add a row (but we want 11,21 to)
            cur_row += div(max(plength+(llength-hasnl)-1,0), cols)
            line_pos -= slength
            write(termbuf, l)
        else
            # We expect to be line after the last valid output line (due to
            # the '\n' at the end of the previous line)
            if curs_row == -1
                if line_pos < slength
                    num_chars = strwidth(l[1:line_pos])
                    curs_row = cur_row + div(indent+num_chars-1, cols)
                    curs_pos = (indent+num_chars-1) % cols + 1
                end
                line_pos -= slength # '\n' gets an extra pos
                cur_row += div(max(indent+(llength-hasnl)-1,0), cols)
                cmove_col(termbuf, indent+1)
                write(termbuf, l)
                # There's an issue if the cursor is after the very right end of the screen. In that case we need to
                # move the cursor to the next line, and emit a newline if needed
                if curs_pos == cols
                    # only emit the newline if the cursor is at the end of the line we're writing
                    if line_pos == 0
                        write(termbuf, "\n")
                        cur_row += 1
                    end
                    curs_row += 1
                    curs_pos = 0
                    cmove_col(termbuf, 1)
                end
            else
                cur_row += div(llength+indent-1, cols)
                cmove_col(termbuf, indent+1)
                write(termbuf, l)
            end
        end
    end

    seek(buf, buf_pos)

    # If we are at the end of the buffer, we need to put the cursor one past the
    # last character we have written

    if curs_row == -1
        curs_pos = ((cur_row == 1 ? plength : indent)+llength-1) % cols + 1
        curs_row = cur_row
    end

    # Same issue as above. TODO: We should figure out
    # how to refactor this to avoid duplcating functionality.
    if curs_pos == cols
        if line_pos == 0
            write(termbuf, "\n")
            cur_row += 1
        end
        curs_row += 1
        curs_pos = 0
        cmove_col(termbuf, 1)
    end

    # Let's move the cursor to the right position
    # The line first
    n = cur_row - curs_row
    if n > 0
        cmove_up(termbuf, n)
    end

    #columns are 1 based
    cmove_col(termbuf, curs_pos+1)

    # Updated cur_row,curs_row
    return InputAreaState(cur_row, curs_row)
end

# Edit functionality
is_non_word_char(c) = c in " \t\n\"\\'`@\$><=:;|&{}()[].,+-*/?%^~"

char_move_left(s::PromptState) = char_move_left(s.input_buffer)
function char_move_left(buf::IOBuffer)
    while position(buf) > 0
        seek(buf, position(buf)-1)
        c = peek(buf)
        (((c & 0x80) == 0) || ((c & 0xc0) == 0xc0)) && break
    end
    pos = position(buf)
    c = read(buf, Char)
    seek(buf, pos)
    c
end

function edit_move_left(buf::IOBuffer)
    if position(buf) > 0
        #move to the next base UTF8 character to the left
        while true
            c = char_move_left(buf)
            if charwidth(c) != 0 || c == '\n' || position(buf) == 0
                break
            end
        end
        return true
    end
    return false
end
edit_move_left(s::PromptState) = edit_move_left(s.input_buffer) && refresh_line(s)

function edit_move_word_left(s)
    if position(s.input_buffer) > 0
        char_move_word_left(s.input_buffer)
        refresh_line(s)
    end
end

char_move_right(s) = char_move_right(buffer(s))
function char_move_right(buf::IOBuffer)
    !eof(buf) && read(buf, Char)
end

function char_move_word_right(buf::IOBuffer, is_delimiter=is_non_word_char)
    while !eof(buf) && is_delimiter(char_move_right(buf))
    end
    while !eof(buf)
        pos = position(buf)
        if is_delimiter(char_move_right(buf))
            seek(buf, pos)
            break
        end
    end
end

function char_move_word_left(buf::IOBuffer, is_delimiter=is_non_word_char)
    while position(buf) > 0 && is_delimiter(char_move_left(buf))
    end
    while position(buf) > 0
        pos = position(buf)
        if is_delimiter(char_move_left(buf))
            seek(buf, pos)
            break
        end
    end
end

char_move_word_right(s) = char_move_word_right(buffer(s))
char_move_word_left(s) = char_move_word_left(buffer(s))

function edit_move_right(buf::IOBuffer)
    if !eof(buf)
        # move to the next base UTF8 character to the right
        while true
            c = char_move_right(buf)
            eof(buf) && break
            pos = position(buf)
            nextc = read(buf,Char)
            seek(buf,pos)
            (charwidth(nextc) != 0 || nextc == '\n') && break
        end
        return true
    end
    return false
end
edit_move_right(s::PromptState) = edit_move_right(s.input_buffer) && refresh_line(s)

function edit_move_word_right(s)
    if !eof(s.input_buffer)
        char_move_word_right(s)
        refresh_line(s)
    end
end

## Move line up/down
# Querying the terminal is expensive, memory access is cheap
# so to find the current column, we find the offset for the start
# of the line.

function edit_move_up(buf::IOBuffer)
    npos = rsearch(buf.data, '\n', position(buf))
    npos == 0 && return false # we're in the first line
    # We're interested in character count, not byte count
    offset = length(bytestring(buf.data[(npos+1):(position(buf))]))
    npos2 = rsearch(buf.data, '\n', npos-1)
    seek(buf, npos2)
    for _ = 1:offset
        pos = position(buf)
        if read(buf, Char) == '\n'
            seek(buf, pos)
            break
        end
    end
    return true
end
function edit_move_up(s)
    changed = edit_move_up(buffer(s))
    changed && refresh_line(s)
    changed
end

function edit_move_down(buf::IOBuffer)
    npos = rsearch(buf.data[1:buf.size], '\n', position(buf))
    # We're interested in character count, not byte count
    offset = length(bytestring(buf.data[(npos+1):(position(buf))]))
    npos2 = search(buf.data[1:buf.size], '\n', position(buf)+1)
    if npos2 == 0 #we're in the last line
        return false
    end
    seek(buf, npos2)
    for _ = 1:offset
        pos = position(buf)
        if eof(buf) || read(buf, Char) == '\n'
            seek(buf, pos)
            break
        end
    end
    return true
end
function edit_move_down(s)
    changed = edit_move_down(buffer(s))
    changed && refresh_line(s)
    changed
end

# splice! for IOBuffer: convert from 0-indexed positions, update the size,
# and keep the cursor position stable with the text
function splice_buffer!{T<:Integer}(buf::IOBuffer, r::UnitRange{T}, ins::String = "")
    pos = position(buf)
    if !isempty(r) && pos in r
        seek(buf, first(r))
    elseif pos > last(r)
        seek(buf, pos - length(r))
    end
    splice!(buf.data, r .+ 1, ins.data) # position(), etc, are 0-indexed
    buf.size = buf.size + sizeof(ins) - length(r)
    seek(buf, position(buf) + sizeof(ins))
end

function edit_replace(s, from, to, str)
    splice_buffer!(buffer(s), from:to-1, str)
end

function edit_insert(s::PromptState, c)
    str = string(c)
    edit_insert(s.input_buffer, str)
    if !('\n' in str) && eof(s.input_buffer) &&
        ((position(s.input_buffer) + length(s.p.prompt) + sizeof(str) - 1) < width(terminal(s)))
        #Avoid full update
        write(terminal(s), str)
    else
        refresh_line(s)
    end
end

function edit_insert(buf::IOBuffer, c)
    if eof(buf)
        write(buf, c)
    else
        splice_buffer!(buf, position(buf):position(buf)-1, string(c))
    end
end

function edit_backspace(s::PromptState)
    if edit_backspace(s.input_buffer)
        refresh_line(s)
    else
        beep(terminal(s))
    end
end
function edit_backspace(buf::IOBuffer)
    if position(buf) > 0
        oldpos = position(buf)
        char_move_left(buf)
        splice_buffer!(buf, position(buf):oldpos-1)
        return true
    else
        return false
    end
end

edit_delete(s) = edit_delete(buffer(s)) ? refresh_line(s) : beep(terminal(s))
function edit_delete(buf::IOBuffer)
    eof(buf) && return false
    oldpos = position(buf)
    char_move_right(buf)
    splice_buffer!(buf, oldpos:position(buf)-1)
    true
end

function edit_werase(buf::IOBuffer)
    pos1 = position(buf)
    char_move_word_left(buf, isspace)
    pos0 = position(buf)
    pos0 < pos1 || return false
    splice_buffer!(buf, pos0:pos1-1)
    true
end
function edit_werase(s)
    edit_werase(buffer(s)) && refresh_line(s)
end

function edit_delete_prev_word(buf::IOBuffer)
    pos1 = position(buf)
    char_move_word_left(buf)
    pos0 = position(buf)
    pos0 < pos1 || return false
    splice_buffer!(buf, pos0:pos1-1)
    true
end
function edit_delete_prev_word(s)
    edit_delete_prev_word(buffer(s)) && refresh_line(s)
end

function edit_delete_next_word(buf::IOBuffer)
    pos0 = position(buf)
    char_move_word_right(buf)
    pos1 = position(buf)
    pos0 < pos1 || return false
    splice_buffer!(buf, pos0:pos1-1)
    true
end
function edit_delete_next_word(s)
    edit_delete_next_word(buffer(s)) && refresh_line(s)
end

function edit_yank(s::MIState)
    edit_insert(buffer(s), s.kill_buffer)
    refresh_line(s)
end

function edit_kill_line(s::MIState)
    buf = buffer(s)
    pos = position(buf)
    killbuf = readline(buf)
    if length(killbuf) > 1 && killbuf[end] == '\n'
        killbuf = killbuf[1:end-1]
        char_move_left(buf)
    end
    s.kill_buffer = s.key_repeats > 0 ? s.kill_buffer * killbuf : killbuf

    splice_buffer!(buf, pos:position(buf)-1)
    refresh_line(s)
end

edit_transpose(s) = edit_transpose(buffer(s)) && refresh_line(s)
function edit_transpose(buf::IOBuffer)
    position(buf) == 0 && return false
    eof(buf) && char_move_left(buf)
    char_move_left(buf)
    pos = position(buf)
    a, b = read(buf, Char), read(buf, Char)
    seek(buf, pos)
    write(buf, b, a)
    return true
end

edit_clear(buf::IOBuffer) = truncate(buf, 0)

function edit_clear(s::MIState)
    edit_clear(buffer(s))
    refresh_line(s)
end

function replace_line(s::PromptState, l::IOBuffer)
    s.input_buffer = l
end

function replace_line(s::PromptState, l)
    s.input_buffer.ptr = 1
    s.input_buffer.size = 0
    write(s.input_buffer, l)
end

history_prev(::EmptyHistoryProvider) = ("", false)
history_next(::EmptyHistoryProvider) = ("", false)
history_search(::EmptyHistoryProvider, args...) = false
add_history(::EmptyHistoryProvider, s) = nothing
add_history(s::PromptState) = add_history(mode(s).hist, s)
history_next_prefix(s, hist) = false
history_prev_prefix(s, hist) = false

function history_prev(s, hist)
    l, ok = history_prev(mode(s).hist)
    if ok
        replace_line(s, l)
        move_input_start(s)
        refresh_line(s)
    else
        beep(terminal(s))
    end
end
function history_next(s, hist)
    l, ok = history_next(mode(s).hist)
    if ok
        replace_line(s, l)
        move_input_end(s)
        refresh_line(s)
    else
        beep(terminal(s))
    end
end

refresh_line(s) = refresh_multi_line(s)
refresh_line(s, termbuf) = refresh_multi_line(termbuf, s)

default_completion_cb(::IOBuffer) = []
default_enter_cb(_) = true

write_prompt(terminal, s::PromptState) = write_prompt(terminal, s, s.p.prompt)
function write_prompt(terminal, s::PromptState, prompt)
    prefix = isa(s.p.prompt_prefix,Function) ? s.p.prompt_prefix() : s.p.prompt_prefix
    suffix = isa(s.p.prompt_suffix,Function) ? s.p.prompt_suffix() : s.p.prompt_suffix
    Base.have_color && write(terminal, prefix)
    write(terminal, prompt)
    write(terminal, Base.text_colors[:normal])
    Base.have_color && write(terminal, suffix)
end
write_prompt(terminal, s::ASCIIString) = write(terminal, s)

normalize_key(key::Char) = string(key)
normalize_key(key::Integer) = normalize_key(char(key))
function normalize_key(key::String)
    '\0' in key && error("Matching \\0 not currently supported.")
    buf = IOBuffer()
    i = start(key)
    while !done(key, i)
        c, i = next(key, i)
        if c == '*'
            write(buf, '\0')
        elseif c == '^'
            c, i = next(key, i)
            write(buf, uppercase(c)-64)
        elseif c == '\\'
            c, i == next(key, i)
            if c == 'C'
                c, i == next(key, i)
                @assert c == '-'
                c, i == next(key, i)
                write(buf, uppercase(c)-64)
            elseif c == 'M'
                c, i == next(key, i)
                @assert c == '-'
                c, i == next(key, i)
                write(buf, '\e')
                write(buf, c)
            end
        else
            write(buf, c)
        end
    end
    return takebuf_string(buf)
end

# Turn a Dict{Any,Any} into a Dict{Char,Any}
# For now we use \0 to represent unknown chars so that they are sorted before everything else
# If we ever actually want to match \0 in input, this will have to be reworked
function normalize_keymap(keymap::Dict)
    ret = Dict{Char,Any}()
    for key in keys(keymap)
        newkey = normalize_key(key)
        current = ret
        i = start(newkey)
        while !done(newkey, i)
            c, i = next(newkey, i)
            if haskey(current, c)
                if !isa(current[c], Dict)
                    println(ret)
                    error("Conflicting Definitions for keyseq " * escape_string(newkey) * " within one keymap")
                end
            elseif done(newkey, i)
                if isa(keymap[key], String)
                    current[c] = normalize_key(keymap[key])
                else
                    current[c] = keymap[key]
                end
                break
            else
                current[c] = Dict{Char,Any}()
            end
            current = current[c]
        end
    end
    ret
end

keymap_gen_body(keymaps, body::Expr, level) = body
keymap_gen_body(keymaps, body::Function, level) = keymap_gen_body(keymaps, :($(body)(s)))
keymap_gen_body(keymaps, body::Char, level) = keymap_gen_body(keymaps, keymaps[body])
keymap_gen_body(keymaps, body::Nothing, level) = nothing
function keymap_gen_body(keymaps, body::String, level)
    if length(body) == 1
        return keymap_gen_body(keymaps, body[1], level)
    end
    current = keymaps
    for c in body
        if haskey(current, c)
            if isa(current[c], Dict)
                current = current[c]
            else
                return keymap_gen_body(keymaps, current[c], level)
            end
        elseif haskey(current, '\0')
            return keymap_gen_body(keymaps, current['\0'], level)
        else
            error("No match for redirected key $body")
        end
    end
    error("No exact match for redirected key $body")
end

keymap_gen_body(a, b) = keymap_gen_body(a, b, 1)
function keymap_gen_body(dict, subdict::Dict, level)
    block = Expr(:block)
    bc = symbol("c" * string(level))
    push!(block.args, :($bc = read(LineEdit.terminal(s), Char)))

    last_if = Expr(:block)
    haskey(subdict, '\0') && push!(last_if.args, keymap_gen_body(dict, subdict['\0'], level+1))
    level == 1 && push!(last_if.args, :(LineEdit.update_key_repeats(s, [$bc])))

    for c in keys(subdict)
        c == '\0' && continue
        cblock = Expr(:if, :($bc == $c))
        cthen = Expr(:block)
        if !isa(subdict[c], Dict)
            cs = [symbol("c" * string(i)) for i=1:level]
            push!(cthen.args, :(LineEdit.update_key_repeats(s, [$(cs...)])))
        end
        push!(cthen.args, keymap_gen_body(dict, subdict[c], level+1))

        push!(cblock.args, cthen)
        push!(cblock.args, last_if)
        last_if = cblock
    end

    push!(block.args, last_if)
    return block
end

update_key_repeats(s, keystroke) = nothing
function update_key_repeats(s::MIState, keystroke)
    s.key_repeats  = s.previous_key == keystroke ? s.key_repeats + 1 : 0
    s.previous_key = keystroke
    return
end

export @keymap

# deep merge where target has higher precedence
function keymap_merge!(target::Dict, source::Dict)
    for k in keys(source)
        if !haskey(target, k)
            target[k] = source[k]
        elseif isa(target[k], Dict)
            keymap_merge!(target[k], source[k])
        else
            # Ignore, target has higher precedence
        end
    end
end

fixup_keymaps!(d, l, s, sk) = nothing
function fixup_keymaps!(dict::Dict, level, s, subkeymap)
    if level > 1
        for d in dict
            fixup_keymaps!(d[2], level-1, s, subkeymap)
        end
    else
        if haskey(dict, s)
            if isa(dict[s], Dict) && isa(subkeymap, Dict)
                keymap_merge!(dict[s], subkeymap)
            end
        else
            dict[s] = deepcopy(subkeymap)
        end
    end
end

function add_specialisations(dict, subdict, level)
    default_branch = subdict['\0']
    if isa(default_branch, Dict)
        for s in keys(default_branch)
            s == '\0' && add_specialisations(dict, default_branch, level+1)
            fixup_keymaps!(dict, level, s, default_branch[s])
        end
    end
end

fix_conflicts!(x) = fix_conflicts!(x, 1)
fix_conflicts!(others, level) = nothing
function fix_conflicts!(dict::Dict, level)
    # needs to be done first for every branch
    if haskey(dict, '\0')
        add_specialisations(dict, dict, level)
    end
    for d in dict
        d[1] == '\0' && continue
        fix_conflicts!(d[2], level+1)
    end
end

keymap_prepare(keymaps::Expr) = keymap_prepare(eval(keymaps))
keymap_prepare(keymaps::Dict) = keymap_prepare([keymaps])
function keymap_prepare{D<:Dict}(keymaps::Array{D})
    push!(keymaps, {"*"=>:(error("Unrecognized input"))})
    keymaps = map(normalize_keymap, keymaps)
    map(fix_conflicts!, keymaps)
    keymaps
end

function keymap_unify(keymaps)
    length(keymaps) == 1 && return keymaps[1]
    ret = Dict{Char,Any}()
    for keymap in keymaps
        keymap_merge!(ret, keymap)
    end
    fix_conflicts!(ret)
    return ret
end

macro keymap(keymaps)
    dict = keymap_unify(keymap_prepare(keymaps))
    body = keymap_gen_body(dict, dict)
    esc(quote
        (s, data) -> begin
            $body
            return :ok
        end
    end)
end

const escape_defaults = merge!(
    {i => nothing for i=1:31}, # Ignore control characters by default
    { # And ignore other escape sequences by default
    "\e*" => nothing,
    "\e[*" => nothing,
    # Also ignore extended escape sequences
    # TODO: Support ranges of characters
    "\e[1**" => nothing,
    "\e[2**" => nothing,
    "\e[3**" => nothing,
    "\e[4**" => nothing,
    "\e[5**" => nothing,
    "\e[6**" => nothing,
    "\e[1~" => "\e[H",
    "\e[4~" => "\e[F",
    "\e[7~" => "\e[H",
    "\e[8~" => "\e[F",
    "\eOA"  => "\e[A",
    "\eOB"  => "\e[B",
    "\eOC"  => "\e[C",
    "\eOD"  => "\e[D",
    "\eOH"  => "\e[H",
    "\eOF"  => "\e[F",
})

function write_response_buffer(s::PromptState, data)
    offset = s.input_buffer.ptr
    ptr = data.response_buffer.ptr
    seek(data.response_buffer, 0)
    write(s.input_buffer, readall(data.response_buffer))
    s.input_buffer.ptr = offset + ptr - 2
    data.response_buffer.ptr = ptr
    refresh_line(s)
end

type SearchState
    terminal
    histprompt
    #rsearch (true) or ssearch (false)
    backward::Bool
    query_buffer::IOBuffer
    response_buffer::IOBuffer
    ias::InputAreaState
    #The prompt whose input will be replaced by the matched history
    parent
    SearchState(terminal, histprompt, backward, query_buffer, response_buffer) =
        new(terminal, histprompt, backward, query_buffer, response_buffer, InputAreaState(0,0))
end

terminal(s::SearchState) = s.terminal

function update_display_buffer(s::SearchState, data)
    history_search(data.histprompt.hp, data.query_buffer, data.response_buffer, data.backward, false) || beep(terminal(s))
    refresh_line(s)
end

function history_next_result(s::MIState, data::SearchState)
    history_search(data.histprompt.hp, data.query_buffer, data.response_buffer, data.backward, true) || beep(terminal(s))
    refresh_line(data)
end

function history_set_backward(s::SearchState, backward)
    s.backward = backward
end

refresh_multi_line(termbuf::TerminalBuffer, term, s::Union(SearchState,PromptState)) = (@assert term == terminal(s); refresh_multi_line(termbuf,s))
function refresh_multi_line(termbuf::TerminalBuffer, s::SearchState)
    buf = IOBuffer()
    write(buf, pointer(s.query_buffer.data), s.query_buffer.ptr-1)
    write(buf, "': ")
    offset = buf.ptr
    ptr = s.response_buffer.ptr
    seek(s.response_buffer, 0)
    write(buf, readall(s.response_buffer))
    buf.ptr = offset + ptr - 1
    s.response_buffer.ptr = ptr
    s.ias = refresh_multi_line(termbuf, s.terminal, buf, s.ias, s.backward ? "(reverse-i-search)`" : "(forward-i-search)`")
end

function refresh_multi_line(s::Union(SearchState,PromptState))
    refresh_multi_line(terminal(s), s)
end

function refresh_multi_line(terminal::UnixTerminal, args...; kwargs...)
    outbuf = IOBuffer()
    termbuf = TerminalBuffer(outbuf)
    ret = refresh_multi_line(termbuf, terminal, args...;kwargs...)
    # Output the entire refresh at once
    write(terminal, takebuf_array(outbuf))
    flush(terminal)
    return ret
end

function reset_state(s::SearchState)
    if s.query_buffer.size != 0
        s.query_buffer.size = 0
        s.query_buffer.ptr = 1
    end
    if s.response_buffer.size != 0
        s.response_buffer.size = 0
        s.response_buffer.ptr = 1
    end
    reset_state(s.histprompt.hp)
end

type HistoryPrompt <: TextInterface
    hp::HistoryProvider
    complete
    keymap_func::Function
    HistoryPrompt(hp) = new(hp, EmptyCompletionProvider())
end

init_state(terminal, p::HistoryPrompt) = SearchState(terminal, p, true, IOBuffer(), IOBuffer())

state(s::MIState, p) = s.mode_state[p]
state(s::PromptState, p) = (@assert s.p == p; s)
mode(s::MIState) = s.current_mode
mode(s::PromptState) = s.p
mode(s::SearchState) = @assert false

# Search Mode completions
function complete_line(s::SearchState, repeats)
    completions, partial, should_complete = complete_line(s.histprompt.complete, s)
    # For now only allow exact completions in search mode
    if length(completions) == 1
        prev_pos = position(s.query_buffer)
        seek(s.query_buffer, prev_pos-sizeof(partial))
        edit_replace(s, position(s.query_buffer), prev_pos, completions[1])
    end
end

function accept_result(s, p)
    parent = state(s, p).parent
    replace_line(state(s, parent), state(s, p).response_buffer)
    transition(s, parent)
end

function copybuf!(dst::IOBuffer, src::IOBuffer)
    n = src.size
    ensureroom(dst, n)
    copy!(dst.data, 1, src.data, 1, n)
    dst.size = src.size
    dst.ptr = src.ptr
end

function enter_search(s::MIState, p::HistoryPrompt, backward::Bool)
    # a bit of hack to help fix #6325
    buf = buffer(s)
    p.hp.last_mode = mode(s)
    p.hp.last_buffer = copy(buf)

    ss = state(s, p)
    ss.parent = mode(s)
    ss.backward = backward
    truncate(ss.query_buffer, 0)
    copybuf!(ss.response_buffer, buf)
    transition(s, p)
end

function setup_search_keymap(hp)
    p = HistoryPrompt(hp)
    pkeymap = {
        "^R"      => :(LineEdit.history_set_backward(data, true); LineEdit.history_next_result(s, data)),
        "^S"      => :(LineEdit.history_set_backward(data, false); LineEdit.history_next_result(s, data)),
        '\r'      => s->accept_result(s, p),
        '\n'      => '\r',
        # Limited form of tab completions
        '\t'      => :(LineEdit.complete_line(s); LineEdit.update_display_buffer(s, data)),
        "^L"      => :(Terminals.clear(LineEdit.terminal(s)); LineEdit.update_display_buffer(s, data)),

        # Backspace/^H
        '\b'      => :(LineEdit.edit_backspace(data.query_buffer) ?
                        LineEdit.update_display_buffer(s, data) : beep(LineEdit.terminal(s))),
        127       => '\b',
        # Meta Backspace
        "\e\b"    => :(LineEdit.edit_delete_prev_word(data.query_buffer) ?
                        LineEdit.update_display_buffer(s, data) : beep(LineEdit.terminal(s))),
        "\e\x7f"  => "\e\b",
        # Word erase to whitespace
        "^W"      => :(LineEdit.edit_werase(data.query_buffer) ?
                        LineEdit.update_display_buffer(s, data) : beep(LineEdit.terminal(s))),
        # ^C and ^D
        "^C"      => :(LineEdit.edit_clear(data.query_buffer);
                       LineEdit.edit_clear(data.response_buffer);
                       LineEdit.update_display_buffer(s, data);
                       LineEdit.reset_state(data.histprompt.hp);
                       LineEdit.transition(s, data.parent)),
        "^D"      => "^C",
        # Other ways to cancel search mode (it's difficult to bind \e itself)
        "^G"      => "^C",
        "\e\e"    => "^C",
        # ^K
        11        => s->transition(s, state(s, p).parent),
        # ^Y
        25        => :(LineEdit.edit_yank(s); LineEdit.update_display_buffer(s, data)),
        # ^U
        21        => :(LineEdit.edit_clear(data.query_buffer);
                       LineEdit.edit_clear(data.response_buffer);
                       LineEdit.update_display_buffer(s, data)),
        # Right Arrow
        "\e[C"    => s->(accept_result(s, p); edit_move_right(s)),
        # Left Arrow
        "\e[D"    => s->(accept_result(s, p); edit_move_left(s)),
        # Up Arrow
        "\e[A"    => s->(accept_result(s, p); edit_move_up(s)),
        # Down Arrow
        "\e[B"    => s->(accept_result(s, p); edit_move_down(s)),
        # ^B
        2         => s->(accept_result(s, p); edit_move_left(s)),
        # ^F
        6         => s->(accept_result(s, p); edit_move_right(s)),
        # Meta B
        "\eb"     => s->(accept_result(s, p); edit_move_word_left(s)),
        # Meta F
        "\ef"     => s->(accept_result(s, p); edit_move_word_right(s)),
        # Ctrl-Left Arrow
        "\e[1;5D" => "\eb",
        # Ctrl-Right Arrow
        "\e[1;5C" => "\ef",
        # ^A
        1         => s->(accept_result(s, p); move_line_start(s); refresh_line(s)),
        # ^E
        5         => s->(accept_result(s, p); move_line_end(s); refresh_line(s)),
        "^Z"      => :(return :suspend),
        # Try to catch all Home/End keys
        "\e[H"    => s->(accept_result(s, p); move_input_start(s); refresh_line(s)),
        "\e[F"    => s->(accept_result(s, p); move_input_end(s); refresh_line(s)),
        # Use ^N and ^P to change search directions and iterate through results
        "^N"      => :(LineEdit.history_set_backward(data, false); LineEdit.history_next_result(s, data)),
        "^P"      => :(LineEdit.history_set_backward(data, true); LineEdit.history_next_result(s, data)),
        # Bracketed paste mode
        "\e[200~" => quote
            ps = LineEdit.state(s, LineEdit.mode(s))
            input = readuntil(ps.terminal, "\e[201~")[1:(end-6)]
            LineEdit.edit_insert(data.query_buffer, input); LineEdit.update_display_buffer(s, data)
        end,
        "*"       => :(LineEdit.edit_insert(data.query_buffer, c1); LineEdit.update_display_buffer(s, data))
    }
    p.keymap_func = @eval @LineEdit.keymap $([pkeymap, escape_defaults])
    keymap = {
        "^R"    => s->(enter_search(s, p, true)),
        "^S"    => s->(enter_search(s, p, false)),
    }
    (p, keymap)
end

keymap(state, p::HistoryPrompt) = p.keymap_func
keymap_data(state, ::HistoryPrompt) = state

Base.isempty(s::PromptState) = s.input_buffer.size == 0

on_enter(s::PromptState) = s.p.on_enter(s)

move_input_start(s) = (seek(buffer(s), 0))
move_input_end(s) = (seekend(buffer(s)))
function move_line_start(s::MIState)
    buf = buffer(s)
    curpos = position(buf)
    curpos == 0 && return
    if s.key_repeats > 0
        move_input_start(s)
    else
        seek(buf, rsearch(buf.data, '\n', curpos))
    end
end
function move_line_end(s::MIState)
    buf = buffer(s)
    eof(buf) && return
    if s.key_repeats > 0
        move_input_end(s)
        return
    end
    pos = search(buffer(s).data, '\n', position(buf)+1)
    if pos == 0
        move_input_end(s)
        return
    end
    seek(buf, pos-1)
end

function commit_line(s)
    move_input_end(s)
    refresh_line(s)
    println(terminal(s))
    add_history(s)
    state(s, mode(s)).ias = InputAreaState(0, 0)
end

const default_keymap =
{
    # Tab
    '\t' => s->begin
        buf = buffer(s)
        # Yes, we are ignoring the possiblity
        # the we could be in the middle of a multi-byte
        # sequence, here but that's ok, since any
        # whitespace we're interested in is only one byte
        i = position(buf)
        if i != 0
            c = buf.data[i]
            if c == '\n' || c == '\t' ||
               # hack to allow path completion in cmds
               # after a space, e.g., `cd <tab>`, while still
               # allowing multiple indent levels
               (c == ' ' && i > 3 && buf.data[i-1] == ' ')
                edit_insert(s, " "^4)
                return
            end
        end
        complete_line(s)
        refresh_line(s)
    end,
    # Enter
    '\r' => quote
        if LineEdit.on_enter(s) || (eof(LineEdit.buffer(s)) && s.key_repeats > 1)
            LineEdit.commit_line(s)
            return :done
        else
            LineEdit.edit_insert(s, '\n')
        end
    end,
    '\n' => '\r',
    # Backspace/^H
    '\b' => edit_backspace,
    127 => '\b',
    # Meta Backspace
    "\e\b" => edit_delete_prev_word,
    "\e\x7f" => "\e\b",
    # ^D
    4 => quote
        if LineEdit.buffer(s).size > 0
            LineEdit.edit_delete(s)
        else
            println(LineEdit.terminal(s))
            return :abort
        end
    end,
    # ^B
    2 => edit_move_left,
    # ^F
    6 => edit_move_right,
    # Meta B
    "\eb" => edit_move_word_left,
    # Meta F
    "\ef" => edit_move_word_right,
    # Ctrl-Left Arrow
    "\e[1;5D" => "\eb",
    # Ctrl-Right Arrow
    "\e[1;5C" => "\ef",
    # Meta Enter
    "\e\r" => :(LineEdit.edit_insert(s, '\n')),
    "\e\n" => "\e\r",
    # Simply insert it into the buffer by default
    "*" => :(LineEdit.edit_insert(s, c1)),
    # ^U
    21 => edit_clear,
    # ^K
    11 => edit_kill_line,
    # ^Y
    25 => edit_yank,
    # ^A
    1 => :(LineEdit.move_line_start(s); LineEdit.refresh_line(s)),
    # ^E
    5 => :(LineEdit.move_line_end(s); LineEdit.refresh_line(s)),
    # Try to catch all Home/End keys
    "\e[H"  => :(LineEdit.move_input_start(s); LineEdit.refresh_line(s)),
    "\e[F"  => :(LineEdit.move_input_end(s); LineEdit.refresh_line(s)),
    # ^L
    12 => :(Terminals.clear(LineEdit.terminal(s)); LineEdit.refresh_line(s)),
    # ^W
    23 => edit_werase,
    # Meta D
    "\ed" => edit_delete_next_word,
    # ^C
    "^C" => s->begin
        try # raise the debugger if present
            ccall(:jl_raise_debugger, Int, ())
        end
        move_input_end(s)
        refresh_line(s)
        print(terminal(s), "^C\n\n")
        transition(s, :reset)
        refresh_line(s)
    end,
    "^Z" => :(return :suspend),
    # Right Arrow
    "\e[C" => edit_move_right,
    # Left Arrow
    "\e[D" => edit_move_left,
    # Up Arrow
    "\e[A" => edit_move_up,
    # Down Arrow
    "\e[B" => edit_move_down,
    # Delete
    "\e[3~" => edit_delete,
    # Bracketed Paste Mode
    "\e[200~" => s->begin
        ps = state(s, mode(s))
        input = readuntil(ps.terminal, "\e[201~")[1:(end-6)]
        input = replace(input, '\r', '\n')
        if position(buffer(s)) == 0
            indent = Base.indentation(input)[1]
            input = Base.unindent(input[(indent+1):end], indent)
        end
        edit_insert(s, input)
    end,
    "^T"      => edit_transpose,
}

function history_keymap(hist)
    return {
        # ^P
        16 => :(LineEdit.history_prev(s, $hist)),
        # ^N
        14 => :(LineEdit.history_next(s, $hist)),
        # Up Arrow
        "\e[A" => :(LineEdit.edit_move_up(s) || LineEdit.history_prev(s, $hist)),
        # Down Arrow
        "\e[B" => :(LineEdit.edit_move_down(s) || LineEdit.history_next(s, $hist)),
        # Page Up
        "\e[5~" => :(LineEdit.history_prev_prefix(s, $hist)),
        # Page Down
        "\e[6~" => :(LineEdit.history_next_prefix(s, $hist))
    }
end

function deactivate(p::Union(Prompt,HistoryPrompt), s::Union(SearchState,PromptState), termbuf)
    clear_input_area(termbuf, s)
    s
end

function activate(p::Union(Prompt,HistoryPrompt), s::Union(SearchState,PromptState), termbuf)
    s.ias = InputAreaState(0, 0)
    refresh_line(s, termbuf)
end

function activate(p::Union(Prompt,HistoryPrompt), s::MIState, termbuf)
    @assert p == s.current_mode
    activate(p, s.mode_state[s.current_mode], termbuf)
end
activate(m::ModalInterface, s::MIState, termbuf) = activate(s.current_mode, s, termbuf)

function transition(s::MIState, mode)
    if mode == :abort
        s.aborted = true
        return
    end
    if mode == :reset
        reset_state(s)
        return
    end
    termbuf = TerminalBuffer(IOBuffer())
    s.mode_state[s.current_mode] = deactivate(s.current_mode, s.mode_state[s.current_mode], termbuf)
    s.current_mode = mode
    activate(mode, s.mode_state[mode], termbuf)
    write(terminal(s), takebuf_array(termbuf.out_stream))
end

function reset_state(s::PromptState)
    if s.input_buffer.size != 0
        s.input_buffer.size = 0
        s.input_buffer.ptr = 1
    end
    s.ias = InputAreaState(0, 0)
end

function reset_state(s::MIState)
    for (mode,state) in s.mode_state
        reset_state(state)
    end
end

const default_keymap_func = @LineEdit.keymap [LineEdit.default_keymap, LineEdit.escape_defaults]

function Prompt(prompt;
    first_prompt = prompt,
    prompt_prefix = "",
    prompt_suffix = "",
    keymap_func = default_keymap_func,
    keymap_func_data = nothing,
    complete = EmptyCompletionProvider(),
    on_enter = default_enter_cb,
    on_done = ()->nothing,
    hist = EmptyHistoryProvider())

    Prompt(prompt, first_prompt, prompt_prefix, prompt_suffix, keymap_func, keymap_func_data,
        complete, on_enter, on_done, hist)
end

run_interface(::Prompt) = nothing

init_state(terminal, prompt::Prompt) = PromptState(terminal, prompt, IOBuffer(), InputAreaState(1, 1), length(prompt.prompt))

function init_state(terminal, m::ModalInterface)
    s = MIState(m, m.modes[1], false, Dict{Any,Any}())
    for mode in m.modes
        s.mode_state[mode] = init_state(terminal, mode)
    end
    s
end

function run_interface(terminal, m::ModalInterface)
    s = init_state(terminal, m)
    while !s.aborted
        p = s.current_mode
        buf, ok, suspend = prompt!(terminal, m, s)
        while suspend
            @unix_only ccall(:jl_repl_raise_sigtstp, Cint, ())
            buf, ok, suspend = prompt!(terminal, m, s)
        end
        s.mode_state[s.current_mode].p.on_done(s, buf, ok)
    end
end

buffer(s::PromptState) = s.input_buffer
buffer(s::SearchState) = s.query_buffer

keymap(s::PromptState, prompt::Prompt) = prompt.keymap_func
keymap_data(s::PromptState, prompt::Prompt) = prompt.keymap_func_data
keymap(ms::MIState, m::ModalInterface) = keymap(ms.mode_state[ms.current_mode], ms.current_mode)
keymap_data(ms::MIState, m::ModalInterface) = keymap_data(ms.mode_state[ms.current_mode], ms.current_mode)

function prompt!(terminal, prompt, s = init_state(terminal, prompt))
    Base.reseteof(terminal)
    raw!(terminal, true)
    enable_bracketed_paste(terminal)
    try
        start_reading(terminal)
        activate(prompt, s, terminal)
        while true
            state = keymap(s, prompt)(s, keymap_data(s, prompt))
            if state == :abort
                stop_reading(terminal)
                return buffer(s), false, false
            elseif state == :done
                stop_reading(terminal)
                return buffer(s), true, false
            elseif state == :suspend
                @unix_only begin
                    stop_reading(terminal)
                    return buffer(s), true, true
                end
            else
                @assert state == :ok
            end
        end
    finally
        raw!(terminal, false) && disable_bracketed_paste(terminal)
    end
end

end # module

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
using Statistics

medium = Medium(:acou_homo1);
ageom = AGeom(medium.mgrid,:xwell);
tgrid = range(0.0,stop=2.0,length=1000)
wav = ricker(10.0, tgrid, tpeak=0.25, );
srcwav = SrcWav(tgrid, ageom, [:P])
update!(srcwav, [:P], wav)

vp0=mean(medium[:vp])
rho0=mean(medium[:rho])
rec1 = GeoPhyInv.Born.mod(vp0=vp0,
            medium_pert=medium,
		     rho0=rho0,
			 ageom=ageom, srcwav=srcwav, tgridmod=tgrid, src_flag=2)

pa=SeisForwExpt(Fdtd(),npw=1,medium=medium,
    ageom=[ageom], srcwav=[srcwav],
        sflags=[2], rflags=[1],
	    tgrid=tgrid, verbose=true );

@time update!(pa);

# least-squares misfit
paerr=GeoPhyInv.VNamedD_misfit(rec1, pa.c.data[1])
err=GeoPhyInv.func_grad!(paerr)

# normalize error
error = err[1]/paerr.ynorm

# desired accuracy?
@test error<1e-2


text:
function TausStep(z::Unsigned, S1::Integer, S2::Integer, S3::Integer, M::Unsigned)
    b = (((z << S1) ⊻ z) >> S2)
    return (((z & M) << S3) ⊻ b)
end

LCGStep(z::Unsigned, A::Unsigned, C::Unsigned) = A * z + C

make_rand_num(::Type{Float64}, tmp) = 2.3283064365387e-10 * Float64(tmp)
make_rand_num(::Type{Float32}, tmp) = 2.3283064f-10 * Float32(tmp)

function next_rand(::Type{FT}, state::NTuple{4, T}) where {FT, T <: Unsigned}
    state = (
        TausStep(state[1], Cint(13), Cint(19), Cint(12), T(4294967294)),
        TausStep(state[2], Cint(2), Cint(25), Cint(4), T(4294967288)),
        TausStep(state[3], Cint(3), Cint(11), Cint(17), T(4294967280)),
        LCGStep(state[4], T(1664525), T(1013904223))
    )
    tmp = (state[1] ⊻ state[2] ⊻ state[3] ⊻ state[4])
    return (
        state,
        make_rand_num(FT, tmp)
    )
end

function gpu_rand(::Type{T}, state, randstate::AbstractVector{NTuple{4, UInt32}}) where T
    threadid = GPUArrays.threadidx_x(state)
    stateful_rand = next_rand(T, randstate[threadid])
    randstate[threadid] = stateful_rand[1]
    return stateful_rand[2]
end

floattype(::Type{T}) where T <: Union{Int64, UInt64} = Float64
floattype(::Type{T}) where T <: Union{Int32, UInt32} = Float32

to_number_range(x::AbstractFloat, ::Type{T}) where T <: Unsigned = T(round(x * typemax(T)))

to_number_range(x::F, ::Type{T}) where {T <: Signed, F <: AbstractFloat} = Base.unsafe_trunc(T, round(((x - F(0.5)) * typemax(T)) * T(2)))

function gpu_rand(::Type{T}, state, randstate::AbstractVector{NTuple{4, UInt32}}) where T <: Integer
    f = gpu_rand(floattype(T), state, randstate)
    return to_number_range(f, T)
end

let rand_state_dict = Dict()
    global cached_state, clear_cache
    clear_cache() = (empty!(rand_state_dict); return)
    function cached_state(x)
        dev = GPUArrays.device(x)
        get!(rand_state_dict, dev) do
            N = GPUArrays.threads(dev)
            res = similar(x, NTuple{4, UInt32}, N)
            copyto!(res, [ntuple(i-> rand(UInt32), 4) for i=1:N])
            res
        end
    end
end
function Random.rand!(A::GPUArray{T}) where T <: Number
    rstates = cached_state(A)
    gpu_call(A, (rstates, A,)) do state, randstates, a
        idx = linear_index(state)
        idx > length(a) && return
        @inbounds a[idx] = gpu_rand(T, state, randstates)
        return
    end
    A
end

Random.rand(X::Type{<: GPUArray}, i::Integer...) = rand(X, Float32, i...)
Random.rand(X::Type{<: GPUArray}, size::NTuple{N, Int}) where N = rand(X, Float32, size...)
Random.rand(X::Type{<: GPUArray{T}}, i::Integer...) where T = rand(X, T, i...)
Random.rand(X::Type{<: GPUArray{T}}, size::NTuple{N, Int}) where {T, N} = rand(X, T, size...)
Random.rand(X::Type{<: GPUArray{T, N}}, size::NTuple{N, Integer}) where {T, N} = rand(X, T, size...)
Random.rand(X::Type{<: GPUArray{T, N}}, size::NTuple{N, Int}) where {T, N} = rand(X, T, size...)

function Random.rand(X::Type{<: GPUArray}, ::Type{ET}, size::Integer...) where ET
    A = similar(X, ET, size)
    rand!(A)
end


text:
using GPUCompiler

module TestRuntime
    # dummy methods
    signal_exception() = return
    malloc(sz) = C_NULL
    report_oom(sz) = return
    report_exception(ex) = return
    report_exception_name(ex) = return
    report_exception_frame(idx, func, file, line) = return

    # for validation
    sin(x) = Base.sin(x)
end

struct TestCompilerParams <: AbstractCompilerParams end
GPUCompiler.runtime_module(::CompilerJob{<:Any,TestCompilerParams}) = TestRuntime

kernel() = nothing

function main()
    source = FunctionSpec(kernel)
    target = NativeCompilerTarget()
    params = TestCompilerParams()
    job = CompilerJob(target, source, params)

    println(GPUCompiler.compile(:asm, job)[1])
end

isinteractive() || main()


text:
#!/usr/bin/env julia

using Luxor

using Test

using Random
Random.seed!(42)

function draw_luxor_demo(fname, flag)
    Drawing("A2landscape", fname)

    origin() # move 0/0 to center
    background("purple")

    setopacity(0.7)     # opacity from 0 to 1
    sethue(0.3, 0.7, 0.9) # sethue sets the color but doesn't change the opacity
    setline(20)         # line width
    # graphics functions use :fill, :stroke, :fillstroke, :clip, or leave blank
    rect(-400, -400, 800, 800, :fill)
    Random.seed!(42)
    randomhue()
    circle(0, 0, 460, :stroke)

    # clip the following graphics to a circle positioned above the x axis
    circle(0, -200, 400, :clip)
    sethue("gold")
    setopacity(0.7)
    setline(10)

    # simple line drawing
    for i in 0:pi/36:2pi - pi/36
        move(0, 0)
        line(cos(i) * 600, sin(i) * 600)
        strokepath()
        move(cos(i) * 300, sin(i) * 300)
        rline(0, 5)
        rline(5, 0)
        rline(0, -5)
        rline(-5, 0)
        rmove(10, 0)
        rline(-10, 0)
        strokepath()
    end

    # finish clipping
    clipreset()

    # here using Mac OS X fonts
    fontsize(20)
    setcolor"turquoise"
    fontface("Optima-ExtraBlack")
    textwidth = textextents("Luxor")[5]
    # move the text by half the width
    flag ? textcentered("Luxor centered", 0, 300) :  textcentred("Luxor centred", 0, 300)
    # check that using a Point instead of x, y works
    pnt = Point(0, 330)
    flag ? textcentered("Centered", pnt) :  textcentred("Centred", pnt)

    textright("Right", 0, 360)
    textright("Right Point", Point(0, 390))

    # Check that this just takes default (:baseline)
    text("Invalid valign", Point(0, 420), valign=:foobar)

    # text on curve
    fontsize(18)
    fontface("Avenir-Black")

    textcurve("THIS IS TEXT ON A CURVE " ^ 14, 0, 550, O)
    @test finish() == true

    # Test that calling finish() twice doesn't segfault.
    @test finish() == false
end

# Try both synonyms of textcentered/textcentred
fname = "luxor-test1.png"
draw_luxor_demo(fname, false)
println("...basic test: output in $(fname)")
fname = "luxor-test2.png"
draw_luxor_demo(fname, true)
println("...finished basic test: output in $(fname)")


text:
using DoubleExponentialFormulas
using LinearAlgebra: norm
using Test

quadde32 = QuadDE(Float32)
quadde64 = QuadDE(Float64)
quaddeBF = QuadDE(BigFloat)

# Test ∫f(x)dx = 0 in [a, b] if a == b
let
    a = rand(Int)

    # Odd function
    I, E = quadde32(x -> x, a, a)
    @test I isa Float32
    @test I ≈ 0
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    I, E = quadde64(x -> x, a, a)
    @test I isa Float64
    @test I ≈ 0
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    I, E = quaddeBF(x -> x, a, a)
    @test I isa BigFloat
    @test I ≈ 0
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    # Even function
    I, E = quadde32(x -> x^2, a, a)
    @test I isa Float32
    @test I ≈ 0
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    I, E = quadde64(x -> x^2, a, a)
    @test I isa Float64
    @test I ≈ 0
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    I, E = quaddeBF(x -> x^2, a, a)
    @test I isa BigFloat
    @test I ≈ 0
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)
end

# Test the anti-symmetric equivalence when switching inteval limits with QuadDE
# ∫f(x)dx in [a, b] = -∫f(x)dx in [b, a]
let
    f(x::AbstractFloat) = exp(-x^2)

    # [-1, 1]
    I1, E1 = quadde32(f, -1, 1)
    I2, E2 = quadde32(f, 1, -1)
    @test I1 ≈ -I2
    @test E1 ≤ sqrt(eps(typeof(I1)))*norm(I1)
    @test E2 ≤ sqrt(eps(typeof(I2)))*norm(I2)

    I1, E1 = quadde64(f, -1, 1)
    I2, E2 = quadde64(f, 1, -1)
    @test I1 ≈ -I2
    @test E1 ≤ sqrt(eps(typeof(I1)))*norm(I1)
    @test E2 ≤ sqrt(eps(typeof(I2)))*norm(I2)

    I1, E1 = quaddeBF(f, -1, 1)
    I2, E2 = quaddeBF(f, 1, -1)
    @test I1 ≈ -I2
    @test E1 ≤ sqrt(eps(typeof(I1)))*norm(I1)
    @test E2 ≤ sqrt(eps(typeof(I2)))*norm(I2)

    # [a, b] (a and b are finite numbers)
    I1, E1 = quadde32(f, -2, 2)
    I2, E2 = quadde32(f, 2, -2)
    @test I1 ≈ -I2
    @test E1 ≤ sqrt(eps(typeof(I1)))*norm(I1)
    @test E2 ≤ sqrt(eps(typeof(I2)))*norm(I2)

    I1, E1 = quadde64(f, -2, 2)
    I2, E2 = quadde64(f, 2, -2)
    @test I1 ≈ -I2
    @test E1 ≤ sqrt(eps(typeof(I1)))*norm(I1)
    @test E2 ≤ sqrt(eps(typeof(I2)))*norm(I2)

    I1, E1 = quaddeBF(f, -2, 2)
    I2, E2 = quaddeBF(f, 2, -2)
    @test I1 ≈ -I2
    @test E1 ≤ sqrt(eps(typeof(I1)))*norm(I1)
    @test E2 ≤ sqrt(eps(typeof(I2)))*norm(I2)

    # [0, ∞]
    I1, E1 = quadde32(f, 0, Inf)
    I2, E2 = quadde32(f, Inf, 0)
    @test I1 ≈ -I2
    @test E1 ≤ sqrt(eps(typeof(I1)))*norm(I1)
    @test E2 ≤ sqrt(eps(typeof(I2)))*norm(I2)

    I1, E1 = quadde64(f, 0, Inf)
    I2, E2 = quadde64(f, Inf, 0)
    @test I1 ≈ -I2
    @test E1 ≤ sqrt(eps(typeof(I1)))*norm(I1)
    @test E2 ≤ sqrt(eps(typeof(I2)))*norm(I2)

    I1, E1 = quaddeBF(f, 0, Inf)
    I2, E2 = quaddeBF(f, Inf, 0)
    @test I1 ≈ -I2
    @test E1 ≤ sqrt(eps(typeof(I1)))*norm(I1)
    @test E2 ≤ sqrt(eps(typeof(I2)))*norm(I2)

    # [a, ∞]
    I1, E1 = quadde32(f, 2, Inf)
    I2, E2 = quadde32(f, Inf, 2)
    @test I1 ≈ -I2
    @test E1 ≤ sqrt(eps(typeof(I1)))*norm(I1)
    @test E2 ≤ sqrt(eps(typeof(I2)))*norm(I2)

    I1, E1 = quadde64(f, 2, Inf)
    I2, E2 = quadde64(f, Inf, 2)
    @test I1 ≈ -I2
    @test E1 ≤ sqrt(eps(typeof(I1)))*norm(I1)
    @test E2 ≤ sqrt(eps(typeof(I2)))*norm(I2)

    I1, E1 = quaddeBF(f, 2, Inf)
    I2, E2 = quaddeBF(f, Inf, 2)
    @test I1 ≈ -I2
    @test E1 ≤ sqrt(eps(typeof(I1)))*norm(I1)
    @test E2 ≤ sqrt(eps(typeof(I2)))*norm(I2)

    # [-∞, 0]
    I1, E1 = quadde32(f, -Inf, 0)
    I2, E2 = quadde32(f, 0, -Inf)
    @test I1 ≈ -I2
    @test E1 ≤ sqrt(eps(typeof(I1)))*norm(I1)
    @test E2 ≤ sqrt(eps(typeof(I2)))*norm(I2)

    I1, E1 = quadde64(f, -Inf, 0)
    I2, E2 = quadde64(f, 0, -Inf)
    @test I1 ≈ -I2
    @test E1 ≤ sqrt(eps(typeof(I1)))*norm(I1)
    @test E2 ≤ sqrt(eps(typeof(I2)))*norm(I2)

    I1, E1 = quaddeBF(f, -Inf, 0)
    I2, E2 = quaddeBF(f, 0, -Inf)
    @test I1 ≈ -I2
    @test E1 ≤ sqrt(eps(typeof(I1)))*norm(I1)
    @test E2 ≤ sqrt(eps(typeof(I2)))*norm(I2)

    # [-∞, b]
    I1, E1 = quadde32(f, -Inf, 2)
    I2, E2 = quadde32(f, 2, -Inf)
    @test I1 ≈ -I2
    @test E1 ≤ sqrt(eps(typeof(I1)))*norm(I1)
    @test E2 ≤ sqrt(eps(typeof(I2)))*norm(I2)

    I1, E1 = quadde64(f, -Inf, 2)
    I2, E2 = quadde64(f, 2, -Inf)
    @test I1 ≈ -I2
    @test E1 ≤ sqrt(eps(typeof(I1)))*norm(I1)
    @test E2 ≤ sqrt(eps(typeof(I2)))*norm(I2)

    I1, E1 = quaddeBF(f, -Inf, 2)
    I2, E2 = quaddeBF(f, 2, -Inf)
    @test I1 ≈ -I2
    @test E1 ≤ sqrt(eps(typeof(I1)))*norm(I1)
    @test E2 ≤ sqrt(eps(typeof(I2)))*norm(I2)

    # [-∞, ∞]
    I1, E1 = quadde32(f, -Inf, Inf)
    I2, E2 = quadde32(f, Inf, -Inf)
    @test I1 ≈ -I2
    @test E1 ≤ sqrt(eps(typeof(I1)))*norm(I1)
    @test E2 ≤ sqrt(eps(typeof(I2)))*norm(I2)

    I1, E1 = quadde64(f, -Inf, Inf)
    I2, E2 = quadde64(f, Inf, -Inf)
    @test I1 ≈ -I2
    @test E1 ≤ sqrt(eps(typeof(I1)))*norm(I1)
    @test E2 ≤ sqrt(eps(typeof(I2)))*norm(I2)

    I1, E1 = quaddeBF(f, -Inf, Inf)
    I2, E2 = quaddeBF(f, Inf, -Inf)
    @test I1 ≈ -I2
    @test E1 ≤ sqrt(eps(typeof(I1)))*norm(I1)
    @test E2 ≤ sqrt(eps(typeof(I2)))*norm(I2)
end

# Parity check for even function
let
    # Even function
    f(x::AbstractFloat) = exp(-x^2)

    # ∫f(x)dx in [a, b] == ∫f(x)dx in [-b, -a] if f(x) is an even function
    I1, E1 = quadde32(f,  0, 1)
    I2, E2 = quadde32(f, -1, 0)
    @test I1 ≈ I2
    @test E1 ≤ sqrt(eps(typeof(I1)))*norm(I1)
    @test E2 ≤ sqrt(eps(typeof(I2)))*norm(I2)

    I1, E1 = quadde64(f,  0, 1)
    I2, E2 = quadde64(f, -1, 0)
    @test I1 ≈ I2
    @test E1 ≤ sqrt(eps(typeof(I1)))*norm(I1)
    @test E2 ≤ sqrt(eps(typeof(I2)))*norm(I2)

    I1, E1 = quaddeBF(f,  0, 1)
    I2, E2 = quaddeBF(f, -1, 0)
    @test I1 ≈ I2
    @test E1 ≤ sqrt(eps(typeof(I1)))*norm(I1)
    @test E2 ≤ sqrt(eps(typeof(I2)))*norm(I2)

    I1, E1 = quadde32(f,  1,  2)
    I2, E2 = quadde32(f, -2, -1)
    @test I1 ≈ I2
    @test E1 ≤ sqrt(eps(typeof(I1)))*norm(I1)
    @test E2 ≤ sqrt(eps(typeof(I2)))*norm(I2)

    I1, E1 = quadde64(f,  1,  2)
    I2, E2 = quadde64(f, -2, -1)
    @test I1 ≈ I2
    @test E1 ≤ sqrt(eps(typeof(I1)))*norm(I1)
    @test E2 ≤ sqrt(eps(typeof(I2)))*norm(I2)

    I1, E1 = quaddeBF(f,  1,  2)
    I2, E2 = quaddeBF(f, -2, -1)
    @test I1 ≈ I2
    @test E1 ≤ sqrt(eps(typeof(I1)))*norm(I1)
    @test E2 ≤ sqrt(eps(typeof(I2)))*norm(I2)

    # ∫f(x)dx in [a, ∞] == ∫f(x)dx in [-∞, -a] if f(x) is an even function
    I1, E1 = quadde32(f,  0, Inf)
    I2, E2 = quadde32(f, -Inf, 0)
    @test I1 ≈ I2
    @test E1 ≤ sqrt(eps(typeof(I1)))*norm(I1)
    @test E2 ≤ sqrt(eps(typeof(I2)))*norm(I2)

    I1, E1 = quadde64(f,  0, Inf)
    I2, E2 = quadde64(f, -Inf, 0)
    @test I1 ≈ I2
    @test E1 ≤ sqrt(eps(typeof(I1)))*norm(I1)
    @test E2 ≤ sqrt(eps(typeof(I2)))*norm(I2)

    I1, E1 = quaddeBF(f,  0, Inf)
    I2, E2 = quaddeBF(f, -Inf, 0)
    @test I1 ≈ I2
    @test E1 ≤ sqrt(eps(typeof(I1)))*norm(I1)
    @test E2 ≤ sqrt(eps(typeof(I2)))*norm(I2)

    I1, E1 = quadde32(f,  Inf,  2)
    I2, E2 = quadde32(f, -2, -Inf)
    @test I1 ≈ I2
    @test E1 ≤ sqrt(eps(typeof(I1)))*norm(I1)
    @test E2 ≤ sqrt(eps(typeof(I2)))*norm(I2)

    I1, E1 = quadde64(f,  Inf,  2)
    I2, E2 = quadde64(f, -2, -Inf)
    @test I1 ≈ I2
    @test E1 ≤ sqrt(eps(typeof(I1)))*norm(I1)
    @test E2 ≤ sqrt(eps(typeof(I2)))*norm(I2)

    I1, E1 = quaddeBF(f,  Inf,  2)
    I2, E2 = quaddeBF(f, -2, -Inf)
    @test I1 ≈ I2
    @test E1 ≤ sqrt(eps(typeof(I1)))*norm(I1)
    @test E2 ≤ sqrt(eps(typeof(I2)))*norm(I2)
end

# Parity check for odd function
let
    # Odd function
    f(x::AbstractFloat) = x*exp(-x^2)

    # ∫f(x)dx in [a, b] == -∫f(x)dx in [-b, -a] if f(x) is an odd function
    I1, E1 = quadde32(f,  0, 1)
    I2, E2 = quadde32(f, -1, 0)
    @test I1 ≈ -I2
    @test E1 ≤ sqrt(eps(typeof(I1)))*norm(I1)
    @test E2 ≤ sqrt(eps(typeof(I2)))*norm(I2)

    I1, E1 = quadde64(f,  0, 1)
    I2, E2 = quadde64(f, -1, 0)
    @test I1 ≈ -I2
    @test E1 ≤ sqrt(eps(typeof(I1)))*norm(I1)
    @test E2 ≤ sqrt(eps(typeof(I2)))*norm(I2)

    I1, E1 = quaddeBF(f,  0, 1)
    I2, E2 = quaddeBF(f, -1, 0)
    @test I1 ≈ -I2
    @test E1 ≤ sqrt(eps(typeof(I1)))*norm(I1)
    @test E2 ≤ sqrt(eps(typeof(I2)))*norm(I2)

    I1, E1 = quadde32(f,  1,  2)
    I2, E2 = quadde32(f, -2, -1)
    @test I1 ≈ -I2
    @test E1 ≤ sqrt(eps(typeof(I1)))*norm(I1)
    @test E2 ≤ sqrt(eps(typeof(I2)))*norm(I2)

    I1, E1 = quadde64(f,  1,  2)
    I2, E2 = quadde64(f, -2, -1)
    @test I1 ≈ -I2
    @test E1 ≤ sqrt(eps(typeof(I1)))*norm(I1)
    @test E2 ≤ sqrt(eps(typeof(I2)))*norm(I2)

    I1, E1 = quaddeBF(f,  1,  2)
    I2, E2 = quaddeBF(f, -2, -1)
    @test I1 ≈ -I2
    @test E1 ≤ sqrt(eps(typeof(I1)))*norm(I1)
    @test E2 ≤ sqrt(eps(typeof(I2)))*norm(I2)

    # ∫f(x)dx in [a, ∞] == ∫f(x)dx in [-∞, -a] if f(x) is an odd function
    I1, E1 = quadde32(f,  0, Inf)
    I2, E2 = quadde32(f, -Inf, 0)
    @test I1 ≈ -I2
    @test E1 ≤ sqrt(eps(typeof(I1)))*norm(I1)
    @test E2 ≤ sqrt(eps(typeof(I2)))*norm(I2)

    I1, E1 = quadde64(f,  0, Inf)
    I2, E2 = quadde64(f, -Inf, 0)
    @test I1 ≈ -I2
    @test E1 ≤ sqrt(eps(typeof(I1)))*norm(I1)
    @test E2 ≤ sqrt(eps(typeof(I2)))*norm(I2)

    I1, E1 = quaddeBF(f,  0, Inf)
    I2, E2 = quaddeBF(f, -Inf, 0)
    @test I1 ≈ -I2
    @test E1 ≤ sqrt(eps(typeof(I1)))*norm(I1)
    @test E2 ≤ sqrt(eps(typeof(I2)))*norm(I2)

    I1, E1 = quadde32(f,  Inf,  2)
    I2, E2 = quadde32(f, -2, -Inf)
    @test I1 ≈ -I2
    @test E1 ≤ sqrt(eps(typeof(I1)))*norm(I1)
    @test E2 ≤ sqrt(eps(typeof(I2)))*norm(I2)

    I1, E1 = quadde64(f,  Inf,  2)
    I2, E2 = quadde64(f, -2, -Inf)
    @test I1 ≈ -I2
    @test E1 ≤ sqrt(eps(typeof(I1)))*norm(I1)
    @test E2 ≤ sqrt(eps(typeof(I2)))*norm(I2)

    I1, E1 = quaddeBF(f,  Inf,  2)
    I2, E2 = quaddeBF(f, -2, -Inf)
    @test I1 ≈ -I2
    @test E1 ≤ sqrt(eps(typeof(I1)))*norm(I1)
    @test E2 ≤ sqrt(eps(typeof(I2)))*norm(I2)
end

# Test the equivalence w/o splitting integral interval
# ∫f(x)dx in [a, b] = ∫f(x)dx in [a, c] + ∫f(x)dx in [c, b]  where a ≤ c ≤ b
let
    f(x::AbstractFloat) = exp(-x^2)

    I,  E  = quadde32(f, -1, 1)
    I1, E1 = quadde32(f, -1, 0)
    I2, E2 = quadde32(f, 0, 1)
    @test I ≈ I1 + I2
    @test E  ≤ sqrt(eps(typeof(I)))*norm(I)
    @test E1 ≤ sqrt(eps(typeof(I1)))*norm(I1)
    @test E2 ≤ sqrt(eps(typeof(I2)))*norm(I2)

    I,  E  = quadde64(f, -1, 1)
    I1, E1 = quadde64(f, -1, 0)
    I2, E2 = quadde64(f, 0, 1)
    @test I ≈ I1 + I2
    @test E  ≤ sqrt(eps(typeof(I)))*norm(I)
    @test E1 ≤ sqrt(eps(typeof(I1)))*norm(I1)
    @test E2 ≤ sqrt(eps(typeof(I2)))*norm(I2)

    I,  E  = quaddeBF(f, -1, 1)
    I1, E1 = quaddeBF(f, -1, 0)
    I2, E2 = quaddeBF(f, 0, 1)
    @test I ≈ I1 + I2
    @test E  ≤ sqrt(eps(typeof(I)))*norm(I)
    @test E1 ≤ sqrt(eps(typeof(I1)))*norm(I1)
    @test E2 ≤ sqrt(eps(typeof(I2)))*norm(I2)

    I,  E  = quadde32(f, -Inf, Inf)
    I1, E1 = quadde32(f, -Inf, 0)
    I2, E2 = quadde32(f, 0, Inf)
    @test I ≈ I1 + I2
    @test E  ≤ sqrt(eps(typeof(I)))*norm(I)
    @test E1 ≤ sqrt(eps(typeof(I1)))*norm(I1)
    @test E2 ≤ sqrt(eps(typeof(I2)))*norm(I2)

    I,  E  = quadde64(f, -Inf, Inf)
    I1, E1 = quadde64(f, -Inf, 0)
    I2, E2 = quadde64(f, 0, Inf)
    @test I ≈ I1 + I2
    @test E  ≤ sqrt(eps(typeof(I)))*norm(I)
    @test E1 ≤ sqrt(eps(typeof(I1)))*norm(I1)
    @test E2 ≤ sqrt(eps(typeof(I2)))*norm(I2)

    I,  E  = quaddeBF(f, -Inf, Inf)
    I1, E1 = quaddeBF(f, -Inf, 0)
    I2, E2 = quaddeBF(f, 0, Inf)
    @test I ≈ I1 + I2
    @test E  ≤ sqrt(eps(typeof(I)))*norm(I)
    @test E1 ≤ sqrt(eps(typeof(I1)))*norm(I1)
    @test E2 ≤ sqrt(eps(typeof(I2)))*norm(I2)
end


text:
using BinaryBuilder, Pkg

name = "HelFEM"
version = v"0.0.1"
sources = [
    GitSource("https://github.com/mortenpi/HelFEM.git", "d3b9ce4d34fcce3296a80aef0ed9b60ce45319c9")
]

script = raw"""
cp -v ${WORKSPACE}/srcdir/HelFEM/julia/CMake.system ${WORKSPACE}/srcdir/HelFEM/CMake.system

# Set up some platform specific CMake configuration. This is more or less borrowed from
# the Armadillo build_tarballs.jl script:
#   https://github.com/JuliaPackaging/Yggdrasil/blob/48d7a89b4aa46b1a8c91269bb138a660f4ee4ece/A/armadillo/build_tarballs.jl#L23-L52
#
# We need to manually set up OpenBLAS because FindOpenBLAS.cmake does not work with BB:
if [[ "${nbits}" == 64 ]] && [[ "${target}" != aarch64* ]]; then
    OPENBLAS="${libdir}/libopenblas64_.${dlext}"
else
    OPENBLAS="${libdir}/libopenblas.${dlext}"
fi

# On i686, it seems that we need to explicitly tell GCC that we want SSE2. Refs:
#   https://github.com/JuliaLang/julia/blob/v1.4.1/src/atomics.h#L8-L10
#   https://stackoverflow.com/questions/16410149/error-sse2-instruction-set-not-enabled-when-including-emmintrin-h
if [[ "${target}" == i686* ]]; then
    export CXXFLAGS="-msse -msse2"
fi

# Compile libhelfem as a static library
cd ${WORKSPACE}/srcdir/HelFEM/
cmake \
    -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} -DCMAKE_BUILD_TYPE=Release \
    -DCMAKE_POSITION_INDEPENDENT_CODE=ON \
    -DUSE_OPENMP=OFF \
    -DHELFEM_BINARIES=OFF -DHELFEM_FIND_DEPS=ON \
    -B build/ -S .
make -C build/ -j${nproc} helfem
make -C build/ install
# Copy the HelFEM license
install_license LICENSE

# Compile the CxxWrap wrapper as a shared library
cd ${WORKSPACE}/srcdir/HelFEM/julia
cmake \
    -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} -DCMAKE_BUILD_TYPE=Release \
    -DBLAS_LIBRARIES=${OPENBLAS} \
    -DJulia_PREFIX=$prefix -DCMAKE_FIND_ROOT_PATH=$prefix -DJlCxx_DIR=$prefix/lib/cmake/JlCxx \
    -B build/ -S .
make -C build/ -j${nproc}
make -C build/ install
"""

# These are the platforms the libcxxwrap_julia_jll is built on.
#
# The libgfortran constraint is necessary because the default is libgfortran 3, but the
# Julia_jll is only available for libgfortran 4.
platforms = [
    Platform("x86_64", "linux"; libc="glibc", cxxstring_abi = "cxx11", libgfortran_version=v"4"),
    Platform("i686", "linux"; libc="glibc", cxxstring_abi = "cxx11", libgfortran_version=v"4"),
    Platform("armv7l", "linux"; libc="glibc", cxxstring_abi = "cxx11", libgfortran_version=v"4"),
    Platform("aarch64", "linux"; libc="glibc", cxxstring_abi = "cxx11", libgfortran_version=v"4"),
    Platform("x86_64", "macos", libgfortran_version=v"4"),
    Platform("x86_64", "windows"; cxxstring_abi = "cxx11", libgfortran_version=v"4"),
    Platform("i686", "windows"; cxxstring_abi = "cxx11", libgfortran_version=v"4"),
    Platform("x86_64", "freebsd", libgfortran_version=v"4"),
]

products = [
    LibraryProduct("libhelfem-cxxwrap", :libhelfem),
]

dependencies = [
    BuildDependency(PackageSpec(name = "Julia_jll", version = "1.4.1")),
    Dependency(PackageSpec(name = "libcxxwrap_julia_jll", version = "0.8.0")),
    Dependency(PackageSpec(name = "armadillo_jll", version = "9.850.1")),
    Dependency(PackageSpec(name = "GSL_jll", version = "2.6.0")),
    Dependency(PackageSpec(name = "OpenBLAS_jll", version = "0.3.9")),
]

build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sevenfourtwo/Yggdrasil

max_issues_repo_name: sevenfourtwo/Yggdrasil

max_forks_repo_name: sevenfourtwo/Yggdrasil


text:
@testset "Constraint Methods" begin

    # Test Collision Avoidance
    T = Float64
    N = 20
    dt = 0.1
    p = 3
    model = DoubleIntegratorGame(p=p)
    probsize = ProblemSize(N,model)
    game_con = GameConstraintValues(probsize)
    radius = 1.0
    add_collision_avoidance!(game_con, radius)

    pu = model.pu
    @test game_con.state_conlist[1].constraints[1].x1 == pu[1]
    @test game_con.state_conlist[1].constraints[2].x1 == pu[1]
    @test game_con.state_conlist[1].constraints[1].x2 == pu[2]
    @test game_con.state_conlist[1].constraints[2].x2 == pu[3]

    @test game_con.state_conlist[2].constraints[1].x1 == pu[2]
    @test game_con.state_conlist[2].constraints[2].x1 == pu[2]
    @test game_con.state_conlist[2].constraints[1].x2 == pu[1]
    @test game_con.state_conlist[2].constraints[2].x2 == pu[3]

    @test game_con.state_conlist[3].constraints[1].x1 == pu[3]
    @test game_con.state_conlist[3].constraints[2].x1 == pu[3]
    @test game_con.state_conlist[3].constraints[1].x2 == pu[1]
    @test game_con.state_conlist[3].constraints[2].x2 == pu[2]

    # Test Spherical Collision Avoidance
    T = Float64
    N = 20
    dt = 0.1
    p = 3
    d = 3
    model = DoubleIntegratorGame(p=p, d=d)
    probsize = ProblemSize(N,model)
    game_con = GameConstraintValues(probsize)
    radius = 1.0
    add_spherical_collision_avoidance!(game_con, radius)

    pu = model.pu
    @test game_con.state_conlist[1].constraints[1].x1 == pu[1]
    @test game_con.state_conlist[1].constraints[2].x1 == pu[1]
    @test game_con.state_conlist[1].constraints[1].x2 == pu[2]
    @test game_con.state_conlist[1].constraints[2].x2 == pu[3]

    @test game_con.state_conlist[2].constraints[1].x1 == pu[2]
    @test game_con.state_conlist[2].constraints[2].x1 == pu[2]
    @test game_con.state_conlist[2].constraints[1].x2 == pu[1]
    @test game_con.state_conlist[2].constraints[2].x2 == pu[3]

    @test game_con.state_conlist[3].constraints[1].x1 == pu[3]
    @test game_con.state_conlist[3].constraints[2].x1 == pu[3]
    @test game_con.state_conlist[3].constraints[1].x2 == pu[1]
    @test game_con.state_conlist[3].constraints[2].x2 == pu[2]

    # Test Control Bounds
    T = Float64
    N = 20
    dt = 0.1
    p = 3
    model = DoubleIntegratorGame(p=p)
    probsize = ProblemSize(N,model)
    game_con = GameConstraintValues(probsize)
    u_min = -10*ones(model.m)
    u_max =  10*ones(model.m)
    u_min[1] = -Inf
    u_max[1] =  Inf
    add_control_bound!(game_con, u_max, u_min)

    @test game_con.control_conlist.constraints[1].u_min == u_min
    @test game_con.control_conlist.constraints[1].u_max == u_max

    # Test Circle Constraint
    T = Float64
    N = 20
    dt = 0.1
    p = 3
    model = DoubleIntegratorGame(p=p)
    px = model.px
    probsize = ProblemSize(N,model)
    game_con = GameConstraintValues(probsize)
    P = 5
    xc = SVector{P,T}([1.0, 2.0, 3.0, 4.0, 5.0])
    yc = SVector{P,T}([-1.0, -2.0, -3.0, -4.0, -5.0])
    radius = SVector{P,T}([0.1, 0.2, 0.3, 0.4, 0.5])
    add_circle_constraint!(game_con, xc, yc, radius)

    @test game_con.state_conlist[1].constraints[1].xi == px[1][1]
    @test game_con.state_conlist[1].constraints[1].yi == px[1][2]

    @test game_con.state_conlist[2].constraints[1].xi == px[2][1]
    @test game_con.state_conlist[2].constraints[1].yi == px[2][2]

    @test game_con.state_conlist[3].constraints[1].xi == px[3][1]
    @test game_con.state_conlist[3].constraints[1].yi == px[3][2]

    # Test WallConstraint
    T = Float64
    N = 20
    dt = 0.1
    p = 3
    model = DoubleIntegratorGame(p=p)
    px = model.px
    n = model.n
    probsize = ProblemSize(N,model)
    game_con = GameConstraintValues(probsize)
    P = 5
    x = 4
    y = 2
    x1 = SVector{P,T}([ 0.0,  0.0,  1.0,  3.0, -2.0])
    y1 = SVector{P,T}([ 1.0, -1.0,  2.0,  2.0,  0.0])
    x2 = SVector{P,T}([ 1.0,  1.0,  2.0,  2.0,  0.0])
    y2 = SVector{P,T}([ 0.0,  0.0,  1.0,  1.0,  0.0])
    xv = SVector{P,T}([ 1.0,  1.0,  1.0,  1.0,  0.0])./sqrt(2)
    yv = SVector{P,T}([ 1.0, -1.0,  1.0, -1.0,  sqrt(2)])./sqrt(2)

    con = WallConstraint(n,x1,y1,x2,y2,xv,yv,x,y)
    walls = [Wall([x1[j],y1[j]], [x2[j],y2[j]], [xv[j],yv[j]]) for j=1:P]
    add_wall_constraint!(game_con, walls )

    @test game_con.state_conlist[1].constraints[1].x == px[1][1]
    @test game_con.state_conlist[1].constraints[1].y == px[1][2]

    @test game_con.state_conlist[2].constraints[1].x == px[2][1]
    @test game_con.state_conlist[2].constraints[1].y == px[2][2]

    @test game_con.state_conlist[3].constraints[1].x == px[3][1]
    @test game_con.state_conlist[3].constraints[1].y == px[3][2]

    # Test helpers
    T = Float64
    N = 20
    dt = 0.1
    p = 3
    model = DoubleIntegratorGame(p=p)
    px = model.px
    n = model.n
    probsize = ProblemSize(N,model)
    game_con = GameConstraintValues(probsize)

    # Add control constraint
    u_min = -10*ones(model.m)
    u_max =  10*ones(model.m)
    add_control_bound!(game_con, u_max, u_min)

    # Add state constraint
    P = 5
    xc = SVector{P,T}([1.0, 2.0, 3.0, 4.0, 5.0])
    yc = SVector{P,T}([-1.0, -2.0, -3.0, -4.0, -5.0])
    radius = SVector{P,T}([0.1, 0.2, 0.3, 0.4, 0.5])
    add_circle_constraint!(game_con, xc, yc, radius)

    # Set constraint parameters
    opts = Options()
    opts.ρ_0 = 1e-3
    opts.ρ_increase = 1e1
    opts.ρ_max = 1e-1
    opts.λ_max = 1e1
    set_constraint_params!(game_con, opts)

    # Test penaly_update and reset_penalties
    @test game_con.control_conval[1].params.μ0 == 1e-3
    @test game_con.state_conval[1][1].params.μ0 == 1e-3

    @test game_con.control_conval[1].params.ϕ == 1e1
    @test game_con.state_conval[1][1].params.ϕ == 1e1

    @test game_con.control_conval[1].params.μ_max == 1e-1
    @test game_con.state_conval[1][1].params.μ_max == 1e-1

    reset_penalties!(game_con)
    @test game_con.control_conval[1].μ[1] == 1e-3*ones(2model.m)
    @test game_con.state_conval[1][1].μ[1] == 1e-3*ones(P)

    penalty_update!(game_con)
    @test game_con.control_conval[1].μ[1] == 1e-2*ones(2model.m)
    @test game_con.state_conval[1][1].μ[1] == 1e-2*ones(P)

    penalty_update!(game_con)
    @test game_con.control_conval[1].μ[1] == 1e-1*ones(2model.m)
    @test game_con.state_conval[1][1].μ[1] == 1e-1*ones(P)

    penalty_update!(game_con)
    penalty_update!(game_con)
    penalty_update!(game_con)
    penalty_update!(game_con)
    @test game_con.control_conval[1].μ[1] == 1e-1*ones(2model.m)
    @test game_con.state_conval[1][1].μ[1] == 1e-1*ones(P)

    reset_penalties!(game_con)
    @test game_con.control_conval[1].μ[1] == 1e-3*ones(2model.m)
    @test game_con.state_conval[1][1].μ[1] == 1e-3*ones(P)

    # Test dual_update and reset_duals
    @test game_con.control_conval[1].params.λ_max == 1e1
    @test game_con.state_conval[1][1].params.λ_max == 1e1

    reset_penalties!(game_con)
    reset_duals!(game_con)
    @test game_con.control_conval[1].λ[1] == 0.0*ones(2model.m)
    @test game_con.state_conval[1][1].λ[1] == 0.0*ones(P)

    dual_update!(game_con)
    @test game_con.control_conval[1].λ[1] == 0.0*ones(2model.m)
    @test game_con.state_conval[1][1].λ[1] == 0.0*ones(P)

    pdtraj = PrimalDualTraj(probsize, dt)
    init_traj!(pdtraj, f=ones, amplitude=1e2)
    evaluate!(game_con, pdtraj.pr)
    @test game_con.control_conval[1].vals[1] == [90*ones(model.m); -110*ones(model.m)]
    dual_update!(game_con)
    @test game_con.control_conval[1].λ[1] == 1e-3 * [90*ones(model.m); 0*ones(model.m)]
    @test game_con.state_conval[1][1].λ[1] == 1e-3 * max.(0, game_con.state_conval[1][1].vals[1])

    init_traj!(pdtraj, f=ones, amplitude=1e5)
    evaluate!(game_con, pdtraj.pr)
    @test game_con.control_conval[1].vals[1] == [(1e5-10)*ones(model.m); -(1e5+10)*ones(model.m)]
    dual_update!(game_con)
    @test game_con.control_conval[1].λ[1] == [opts.λ_max*ones(model.m); 0*ones(model.m)]

    reset_duals!(game_con)
    @test game_con.control_conval[1].λ[1] == 0.0*ones(2model.m)
    @test game_con.state_conval[1][1].λ[1] == 0.0*ones(P)

end


text:
using Documenter, MethodOfLines

makedocs(
    sitename="MethodOfLines.jl",
    authors="Chris Rackauckas, Alex Jones et al.",
    clean=true,
    doctest=false,
    modules=[MethodOfLines],

    format=Documenter.HTML(assets=["assets/favicon.ico"],
                           canonical="https://methodoflines.sciml.ai/stable/"),

    pages=[
        "MethodOfLines.jl: Automated Finite Difference for Phyiscs-Informed Learning" => "index.md",
        "Tutorials" => ["tutorials/brusselator.md", "tutorials/icbc_sampled.md"],
        "MOLFiniteDifference" => "MOLFiniteDifference.md",
        "Boundary Conditions" => "boundary_conditions.md",
        "How it works" => "howitworks.md",
        "Notes for developers: Implement a scheme" => "devnotes.md",
        "Generated Examples" => ["generated/bruss_code.md", "generated/bruss_ode_eqs.md"]
        # "Tutorial: Burgers" => "tutorials/burgers.md",
        # "Tutorial: 1D Linear Diffusion" => "tutorials/1d_linear_diffusion.md",
        # "Tutorial: 1D Non-Linear Diffusion" => "tutorials/1d_nonlinear_diffusion.md",
        # "Tutorial: 2D Diffusion" => "tutorials/2d_diffusion.md",
        # "Tutorial: 1D Linear Convection" => "tutorials/1d_linear_convection.md",
        # "Tutorial: 1D Higher Order" => "tutorials/1d_higher_order.md",
        # "Tutorial: Stationary Nonlinear Problems" => "tutorials/stationary_nonlinear_problems.md",
        # "Tutorial: 1D Partial DAE" => "tutorials/1d_partial_DAE.md",

     ]
)

deploydocs(
    repo="github.com/SciML/MethodOfLines.jl";
    push_preview=true
)


text:
scene, layout = layoutscene(resolution = (500, 100), font = "CMU Serif Roman");
ax = layout[1, 1] = LAxis(scene)
hidexdecorations!(ax); hideydecorations!(ax);

xs = 0:0.5:10
ys = sin.(xs)
lin = lines!(ax, xs, ys, color = :blue)
sca = scatter!(ax, xs, ys, color = :red, markersize = 15px)

leg = layout[1, 2] = LLegend(scene, [lin, sca, [lin, sca]], ["a line", "some dots", " "], labelsize = 12);
leg.grid.content[1].content[3, 2] = MakieTeX.LTeX(leg.child, raw"\dot x = \frac{\partial x}{\partial t}");
scene


text:
# ---
# title: 397. Integer Replacement
# id: problem397
# author: Indigo
# date: 2021-06-20
# difficulty: Medium
# categories: Math, Bit Manipulation
# link: <https://leetcode.com/problems/integer-replacement/description/>
# hidden: true
# ---
#
# Given a positive integer `n`, you can apply one of the following operations:
#
#   1. If `n` is even, replace `n` with `n / 2`.
#   2. If `n` is odd, replace `n` with either `n + 1` or `n - 1`.
#
# Return _the minimum number of operations needed for`n` to become `1`_.
#
#
#
# **Example 1:**
#
#
#
#     Input: n = 8
#     Output: 3
#     Explanation: 8 -> 4 -> 2 -> 1
#
#
# **Example 2:**
#
#
#
#     Input: n = 7
#     Output: 4
#     Explanation: 7 -> 8 -> 4 -> 2 -> 1
#     or 7 -> 6 -> 3 -> 2 -> 1
#
#
# **Example 3:**
#
#
#
#     Input: n = 4
#     Output: 2
#
#
#
#
# **Constraints:**
#
#   * `1 <= n <= 231 - 1`
#
#
## @lc code=start
using LeetCode

function integer_replacement(n::Int)
    res = 0
    while n != 1
        if iseven(n)
            n >>= 1
        else
            n += (n & 2 == 0 || n == 3) ? -1 : 1
        end
        res += 1
    end
    return res
end
## @lc code=end


text:
arr = IntArray(4)
push!(arr, 1)
push!(arr, 2)
push!(arr, 3)
push!(arr, 4)
use_as_list(arr, 4)
@test size(arr) == 4

@test arr[0] == 1
@test arr[1] == 2
@test arr[2] == 3
@test arr[3] == 4

a = [1,2,3,4]
b = IntArray(a)
@test b[0] == 1
@test b[1] == 2
@test b[2] == 3
@test b[3] == 4
@test size(b) == 4

c = IntArray(arr)
@test b[0] == 1
@test b[1] == 2
@test b[2] == 3
@test b[3] == 4
@test size(c) == 4

d = to_native_int_array(c)
@test d[1] == 1
@test d[2] == 2
@test d[3] == 3
@test d[4] == 4
@test length(d) == 4


text:
# using Pkg
# pkg"activate .."
using ReactionEngine
using Documenter
include("../src/Utils.jl")
include("../src/IdealGas.jl")
include("../src/Inspect.jl")
include("../src/Cstr.jl")
include("../src/Plug.jl")
include("../src/ThermoProbe.jl")
include("../src/Batch.jl")
include("../src/Equil.jl")
include("../src/Transport.jl")

DocMeta.setdocmeta!(ReactionEngine, :DocTestSetup, :(using ReactionEngine); recursive=true)

makedocs(;
    modules=[ReactionEngine, IdealGas, Inspect, Cstr, Plug, ThermoProbe, Batch, Equil, Transport],
    authors="Vinod Janardhanan",
    repo="https://github.com/vinodjanardhanan/ReactionEngine.jl/blob/{commit}{path}#{line}",
    sitename="ReactionEngine.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://vinodjanardhanan.github.io/ReactionEngine.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Models"=>[
            "Inspect" => "Models/inspect.md"
            "Thermoprobe"=>"Models/tprobe.md"
            "Equilibrate"=>"Models/equil.md"
            "Batch"=>"Models/batch.md"
            "Cstr"=>"Models/cstr.md"
            "Plug"=>"Models/plug.md"
            "Sensitivity"=>"Models/sens.md"
        ],
        "Library"=>[
            "Thermodynamic Properties"=>"Library/thermo.md"
            "Surface Chemistry"=>"Library/schem.md"
            "Transport Properties"=>"Library/trans.md"
        ]

    ],
)

deploydocs(;
    repo="github.com/vinodjanardhanan/ReactionEngine.jl",
    devbranch="main",
)


text:
#this file holds functions related to data processing

#reads in a CSV type
function readSingleTypeCSV(path::String; colType=Union{Missing, String},
    streamer::Function = I_CSV_STREAMER)::DataFrame

  #count the number of cols
  #iStream = streamer(path)
  #df::DataFrame = CSV.read(iStream, strings=:raw, categorical=false, silencewarnings=true, rows=2)
  #close(iStream)

  b::IOBuffer = gzopen(path) do f
    s::String = "$(readline(f))\n$(readline(f))\n$(readline(f))"
    println
    b = IOBuffer(s)
    (b)
  end
  df::DataFrame = CSV.read(b, strings=:raw,
    categorical=false, silencewarnings=true, strict=false)

  nCols::Int = size(df,2)
  #println("nCols: $nCols")

  #set the type
  types::Dict = Dict(i=>Union{String,Missing} for i ∈ 1:nCols)
  iStream = streamer(path)
  #df = load(iStream,  colparsers=types) |> DataFrame
  b = IOBuffer(read(iStream))
  df = CSV.read(#=iStream=#b, strings=:raw, categorical=false, types=types,
     silencewarnings=true, strict=false)
  #df = uCSV.read(iStream, coltypes=types, allowmissing=true, escape='\\', quotes='\"')
  close(iStream)

  #println("$(df[1:3, :NAME])")

  return df
end

#same as the read version only returns nothing and just checks for errors
function validateSingleTypeCSV(path::String; colType=Union{Missing, String},
  streamer::Function = I_CSV_STREAMER,
  missingStrings::Vector{String}=MISSING_STRINGS,
  verifyFields::Vector{Symbol} = VERIFY_FIELDS)::Nothing

  #count the number of cols

  #iStream = streamer(path)
  b::IOBuffer = gzopen(path) do f
    b = IOBuffer("$(readline(f))\n$(readline(f))")
    (b)
  end
  df::DataFrame = CSV.read(b, strings=:raw, categorical=false,
    silencewarnings=true, strict=false)

  dfNames::Vector{Symbol} = (s::Symbol->Symbol(lowercase(string(s)))).(names(df))
  for f::Symbol ∈ verifyFields
    if f ∉ dfNames
      @error "Field $f not found in file $path, names found: $dfNames"
    end
  end

  nCols::Int = size(df,2)

  #set the type
  types::Dict = Dict(i=>Union{String,Missing} for i ∈ 1:nCols)
  iStream = streamer(path)
  CSV.read(iStream, strings=:raw, categorical=false, types=types,
    missingstrings=missingStrings, strict=false, silencewarnings=true)
  #df = load(iStream,  colparsers=types) |> DataFrame
  close(iStream)

  return nothing

end

#this fucntion is almsot a script which analyzes data
function preProcessData(; refreshData::Bool=true, validateDownloads::Bool = true,
  nccsPath::String = NCCS_PATH,  refreshMetadata::Bool = true, reserializeData::Bool = true,
  refreshSummary = true)::Vector{MetaDictionary}

  downloadNCCSCSVs(NCCS_PATH, refreshData=refreshData, #=delay=0.001,=#
    validate=validateDownloads, mode=:parallel)
  metadata::Vector{MetaDictionary} = getMetadata(refreshMetadata = refreshMetadata)
  loadNCCSFiles(metadata, reserializeData=reserializeData, serializeOnly = true, mode=:parallel)
  consolidateSummary(metadata, refreshSummary=refreshSummary)

  #=println("fileName: $(metadata[1].fileName)")
  willFail(metadata[1])=#

  return metadata
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
using Unitful

c = Circuit(zeros(typeof(1.0/u"nH"), 3, 3),
            zeros(typeof(1.0/u"Ω"), 3, 3),
            zeros(typeof(1.0u"fF"), 3, 3), 0:2)

set_inductance!(c, 0, 1, 10u"nH")
set_inductance!(c, 0, 2, 10u"nH")
set_capacitance!(c, 0, 1, 100u"fF")
set_capacitance!(c, 0, 2, 110u"fF")
set_capacitance!(c, 1, 2, 1u"fF")

pso = PSOModel(c, [(0,1), (0,2)], ["q1", "q2"])
bbox = Blackbox(collect(1:0.1:10)u"GHz", pso)

@test dimension(eltype(eltype(bbox.Y))) == dimension(u"Ω^-1")


text:
using CRCBS

## set up the environment
vtx_grid = initialize_dense_vtx_grid(4,4) # 4 x 4 grid world
#  1   2   3   4
#  5   6   7   8
#  9  10  11  12
# 13  14  15  16
env = construct_factory_env_from_vtx_grid(vtx_grid)

## Define the initial conditions of the robots
starts = [1,4]
goals = [13,16]

cost_model = FullCostModel(maximum,TravelTime())

prob = init_mapf_problem(env,starts,goals,cost_model)

## define solver
solver = CBSSolver() # PIBTPlanner()

# solve the problem
solution, cost = solve!(solver,prob)
# check if the problem was solved to optimality
@show feasible_status(solver)
@show optimal_status(solver)


text:
module Modalities

"""
    NeuroModality

Supertype for traits corresponding to neuroinformatics modalities.
"""
abstract type NeuroModality end

"""
    NeuroNone

Indicates the complete absence of any `NeuroModality`.
"""
struct NeuroNone <: NeuroModality end

"""
    modality(x) -> NeuroModality

Returns the underlying `NeuroModality`, if present. Otherwise returns
`NeuroNone`.
"""
modality(x::T) where {T} = NeuroNone

modality(x::NeuroModality) = x

"""
    HasTimeComponent

Trait type for whether a modality has a time component.
"""
struct HasTimeComponent{T} end

HasTimeComponent(::T) where {T} = HasTimeComponent(T)

HasTimeComponent(::Type{T}) where {T} = HasTimeComponent{false}()

"""
    Imaging
"""
abstract type Imaging <: NeuroModality end

"""
    Localizer Scan
"""
struct LocalizerScan <: Imaging end

"""
    MagneticResonanceImaging (MRI)

Magnetic resonance imaging.
"""
abstract type MagneticResonanceImaging <: Imaging end

include("smri.jl")
include("dwi.jl")
include("fmri.jl")
include("xray.jl")
include("electrophysiology.jl")

###
### Diffusion
###

abstract type DiffusionMRI <: MagneticResonanceImaging end

struct SingleShellDTI <: DiffusionMRI end

struct MultiShellDTI <: DiffusionMRI end

struct Microscopy <: Imaging end

end


text:
using ImageSegmentation, ImageCore, ImageFiltering, Test, SimpleWeightedGraphs, LightGraphs, StaticArrays, DataStructures
using RegionTrees: isleaf, Cell, split!
using MetaGraphs: MetaGraph, clear_props!, get_prop, has_prop, set_prop!, props, vertices

using Documenter
Base.VERSION >= v"1.6" && doctest(ImageSegmentation, manual = false)
@test isempty(detect_ambiguities(ImageSegmentation))

include("core.jl")
include("region_growing.jl")
include("felzenszwalb.jl")
include("fast_scanning.jl")
include("flood_fill.jl")
include("watershed.jl")
include("region_merging.jl")
include("meanshift.jl")
include("merge_segments.jl")


text:
## load packages
using Plots, LaTeXStrings, VortexDistributions
gr(titlefontsize=12,size=(500,300),transpose=true,colorbar=false)

function psimovie(sol,sim)
    @unpack X,t = sim; x,y = X
    anim = @animate for i in eachindex(t)
        ψ = xspace(sol[i],sim)
        showpsi(x,y,ψ)
    end
    return anim
end

## Set up parameters and ground state
# Units of ξ for length, 1/μ for time
# related by μ ≡ ħ²/mξ² = mc^2 where
# c = ħ/mξ is the speed of sound of
# the uniform system.

# Initialize simulation
L = (200.,200.)
N = (512,512)
sim = Sim(L,N)
@unpack_Sim sim

# set simulation parameters
μ = 1.0
g = 0.01
γ = 0.5
ti = 0.0
tf = .3pi
Nt = 150
t = LinRange(ti,tf,Nt)

# Useful state functions
ψ0(x,y,μ,g) = sqrt(μ/g)*sqrt(max(1.0-V(x,y,0.0)/μ,0.0)+im*0.0)
healinglength(x,y,μ,g) = 1/sqrt(g*abs2(ψ0(x,y,μ,g)))

## Make initial state
x,y = X
kx,ky = K
ψi = ψ0.(x,y',μ,g)
ϕi = kspace(ψi,sim)
@pack_Sim! sim

## evolve
sol = runsim(sim)

## plot ground state
ϕg = sol[end]
ψg = xspace(ϕg,sim)
showpsi(x,y,ψg)

## periodic dipole phase
# Billam et al, PRL 112, 145301 (2014), Supplemental
# Methods
H(x) = x > 0. ? 1.0 : 0.0
shift(x,xi) = x - xi
tans(x,xk) = tan((shift(x,xk) - π)*0.5)
tanhs(x,xk,j) = tanh((shift(x,xk) + 2*π*j)*0.5)

function kernel(x,y,xp,yp,xn,yn,j)
    return atan(tanhs(y,yn,j)*tans(x,xn)) -
    atan(tanhs(y,yp,j)*tans(x,xp))
end

# Dimensionless form
function θd(x,y,dip)
    vp,vn = dip
    xp,yp,_ = rawData(vp)
    xn,yn,_ = rawData(vn)
    s = 0.0
    for j = -5:5
        s += kernel(x,y,xp,yp,xn,yn,j)
    end
    return s + π*(H(shift(x,xp)) - H(shift(x,xn))) - y*(xp - xn)/(2*π)
end
# arbitrary domains and dipole sizes:
function thetad(x,y,xp,yp,xn,yn)
    s = 0.0
    for j = -5:5
        s += kernel(x,y,xp,yp,xn,yn,j)
    end
    s += π*(H(shift(x,xp)) - H(shift(x,xn))) - y*(xp - xn)/(2*π)
    return s - x*H(abs(yp - yn) - π) + y*H(abs(xp - xn) - π)
end

function Thetad(x,y,xp,yp,xn,yn)
    Lx = x[end]-x[1]
    Ly = y[end]-y[1]
    return @. angle(exp(im*thetad.(x*2*pi/Lx,y'*2*pi/Ly,xp*2*pi/Lx,yp*2*pi/Ly,xn*2*pi/Lx,yn*2*pi/Ly)))
end

## initial dipole, with periodic phase
ψv = copy(ψg)
d = 14
ξv = healinglength(0.,0.,μ,g)
# can use default ξ = 1
xp,yp = -50.,d/2
xn,yn = -50.,-d/2
pv = PointVortex(xp,yp,1)
nv = PointVortex(xn,yn,-1)

dipole = [nv;pv]
psi = Torus(ψv,x,y) # set field topology for VortexDistributions
vortex!(psi,dipole) # make dipole
ψv = abs.(psi.ψ).*exp.(im*Thetad(x,y,xp,yp,xn,yn))
showpsi(x,y,ψv)

## make sure phase is periodic
using Test
testphase = angle.(ψv)
@test testphase[:,1] ≈ testphase[:,end]
@test testphase[1,:] ≈ testphase[end,:]

## evolve dipole with weak damping
# Set simulation parameters
c = sqrt(μ)
tf = 3*L[1]/c # a short run for testing evolutoin and detection
γ = 0.03

t = LinRange(ti,tf,Nt)
ϕi = kspace(ψv,sim)
alg = Vern7()

@pack_Sim! sim

# remove boundary artifacts with small γ
@time solv = runsim(sim)

## plot and animate
ψd = xspace(solv[end],sim)
showpsi(x,y,ψd)

anim = psimovie(solv,sim)
gif(anim,"./examples/dipole_damping.gif",fps=25)

## Vortex detection test
Nt = 100
d = zeros(Nt)

for i=1:Nt
    psi = Torus(xspace(solv[i],sim),x,y)
    vort = findvortices(psi) |> rawData
    d[i] = abs(vort[2,2]-vort[1,2])
end

plot(t[1:Nt],d)

## Hamiltonian evolution
γ = 0.0
tf = L[1]/c/5
t = LinRange(ti,tf,Nt)
ϕi = kspace(ψv,sim)
@pack_Sim! sim

solh = runsim(sim)

## plot
ψdh = xspace(solh[end],sim)
showpsi(x,y,ψdh)

anim = psimovie(solh,sim)
gif(anim,"./examples/dipole_hamiltonian.gif",fps=25)

## energy densities
function showenergies(ψ)
    et,ei,ec = energydecomp(ψ)
    p1 = heatmap(x,y,log10.(ei),aspectratio=1)
    xlabel!(L"x/\xi");ylabel!(L"y/\xi")
    title!("Incompressible")
    p2 = heatmap(x,y,log10.(ec),aspectratio=1)
    xlabel!(L"x/\xi");ylabel!(L"y/\xi")
    title!("Compressible")
    p = plot(p1,p2,size=(600,300))
    return p
end

K2 = k2(K)

anim = @animate for i in eachindex(t)
    ψ = xspace(solh[i],sim)
    psi = XField(ψ,X,K,K2)
    showenergies(psi)
end

gif(anim,"./examples/dipoleenergies.gif",fps=30)

## energy totals

function venergy(ϕ,sim,t)
    @unpack g,X = sim; x,y = X
    ψ = xspace(ϕ,sim)
    @. ψ *= V(x,y',t)
    return kspace(ψ,sim)
end

function ienergy(ϕ,sim)
    @unpack g,X = sim; x,y = X
    ψ = xspace(ϕ,sim)
    @. ψ *= 0.5*g*abs2(ψ)
    return kspace(ψ,sim)
end

function kinetic(ϕ,sim)
    @unpack espec,K = sim; kx,ky = K
    dkx,dky = kx[2]-kx[1],ky[2]-ky[1]
    return sum(espec.*abs2.(ϕ))*dkx*dky |> real
end

dx = diff(x)[1]; dy = diff(y)[1]
dkx = diff(kx)[1]; dky = diff(ky)[1]
n0 = μ/g

function energies(sol)
    t = sol.t
    # H = zero(t)
    Ekall = zero(t)
    Ekhy = zero(t)
    Ei = zero(t)
    Ec = zero(t)
    Ev = zero(t)
    Eint = zero(t)
    Eqp = zero(t)
    Et = zero(t)
    Natoms = zero(t)

    for (i,t) in enumerate(t)
    ϕ = sol[i]
    ψ = xspace(ϕ,sim)
    psi = XField(ψ,X,K,K2)
    ek,ei,ec = energydecomp(psi)
    Ekhy[i] = sum(ek)*dx*dy |> real
    Ei[i] = sum(ei)*dx*dy |> real
    Ec[i] = sum(ec)*dx*dy |> real
    Ekall[i] = kinetic(ϕ,sim)

    Ev[i] = sum(@. V(x,y',t)*abs2(ψ))*dx*dy |> real
    Eint[i] = sum(@. g/2*(abs2(ψ) - n0)^2)*dx*dy |> real
    Eqp[i] = Ekall[i] - Ekhy[i]

    Et[i] = Ekall[i] + Ev[i] + Eint[i]
    Natoms[i] = sum(abs2.(ψ))*dx*dy
    end
    return Ei,Ec,Ekhy,Eqp,Ev,Eint,Et,Ekall,Natoms
end

@time Ei,Ec,Ekhy,Eqp,Ev,Eint,Et,Ekall,Natoms = energies(solh)

## plot energies
# simple breakdown
plot(t,Et./Natoms,label=L"E_t",legend=:left)
plot!(t,Eint./Natoms,label=L"E_{int}")
plot!(t,Ekall./Natoms,label=L"E_{kall}")
ylims!(0,1.1*Et[1]./Natoms[1])
xlabel!(L"t")

# decomposition
plot(t,(Ekall .+ Eint)./Natoms,label=L"E_{kall}+E_{int}",w=5,c=:gray,alpha=0.3)
plot!(t,Et./Natoms,label=L"E_t",legend=:left)
plot!(t,Ei./Natoms,label=L"E_k^i")
plot!(t,Ec./Natoms,label=L"E_k^c")
plot!(t,Eint./Natoms,label=L"E_{int}")
plot!(t,Eqp./Natoms,label=L"E_{qp}")
plot!(t,Ekall./Natoms,label=L"E_{kall}")
xlabel!(L"t")

# if we need to save data:
# @save "./examples/dipoledecay.jld2" solv.u sim

# another example for saving data as individual files
# using a call back in OrdinaryDiffEq


text:
# BATs pdf, cdf, quantile, sampler defined below
# Parameter layout: (κ₀, τ₀, ϕ₀, κ₁, τ₁, ϕ₁, ν)
# (κ₀, τ₀, ϕ₀) are shape, scale, location for lower Tail
# (κ₁, τ₁, ϕ₁) are shape, scale, location for upper Tail

# Psi function, its derivative, and its inverse. Branched like this to stabilize
# numerics for large argument x.
Ψ(x)  = (x < 25.0) ? log1p(exp(x))           : x
dΨ(x) = (x < 25.0) ? exp(x)/(1.0+exp(x))     : 1.0
iΨ(x) = (x < 25.0) ? log(max(0.0, expm1(x))) : x # avoids an error with log of negative number

# One half of the H function, and the extension relevant for the density.
# Trying to be careful here about numerics, but it is kind of fishy as κ->0.
function H_part(x, κ, τ, ϕ)
  if isapprox(κ, 0.0, atol = 1.0e-8) # cutoff chosen by eye.
      a = Ψ((x - ϕ) / τ)
      return exp(a) - 0.5 * (a^2) * exp(a) * κ
  end
  (1.0 + κ * Ψ((x - ϕ) / τ))^(1.0 / κ)
end

H(x, κ₀, τ₀, ϕ₀, κ₁, τ₁, ϕ₁) = H_part(x, κ₁, τ₁, ϕ₁) - H_part(-x, κ₀, τ₀, -ϕ₀)
dH_part(x, κ, τ, ϕ) = ForwardDiff.derivative(z -> H_part(z, κ, τ, ϕ), x)

# Full PDF in Distributions.jl format. Note that the numerics
# here can be concerning: the return line out_1*(out_2+out_3) can commonly look
# like subnormal_tiny * (super_huge + super_huge), which is a floating point
# nightmare.
function Distributions.pdf(d::BulkAndTailsDist, x::Real)
  # Check support.
  ((x ≥ maximum(d)) || (x ≤ minimum(d))) && return 0.0
  # Unpack parameters for readability.
  (κ₀, τ₀, ϕ₀, κ₁, τ₁, ϕ₁, ν) = params(d)
  # Evaluate the density if we're inside the defined support.
  out_1 = Distributions.pdf(TDist(ν),H(x,κ₀,τ₀,ϕ₀,κ₁,τ₁,ϕ₁))
  out_2 = dH_part(x, κ₁, τ₁, ϕ₁)
  out_3 = dH_part(-x, κ₀, τ₀, -ϕ₀)
  max(0.0, out_1*(out_2+out_3)) # numerically dubious.
end

# Note: this t-cdf does not work with autodiff. Would need hand-coded tcdf function.
function Distributions.cdf(d::BulkAndTailsDist, x::Real)
  if (x ≤ minimum(d))
    return 0.0
  elseif (x ≥ maximum(d))
    return 1.0
  else
    (κ₀, τ₀, ϕ₀, κ₁, τ₁, ϕ₁, ν) = params(d)
    return Distributions.cdf(TDist(ν),H(x,κ₀,τ₀,ϕ₀,κ₁,τ₁,ϕ₁))
  end
end

function Distributions.quantile(d::BulkAndTailsDist, p::Real)
  Roots.find_zero(x -> Distributions.cdf(d,x) - p, (minimum(d),maximum(d)), Roots.Bisection())
end

function Distributions.logpdf(d::BulkAndTailsDist, x::Real)
  # Check support.
  ((x ≥ maximum(d)) || (x ≤ minimum(d))) && return -Inf
  # Unpack parameters for readability.
  (κ₀, τ₀, ϕ₀, κ₁, τ₁, ϕ₁, ν) = params(d)
  # Evaluate the density if we're inside the defined support.
  out_1 = Distributions.logpdf(TDist(ν),H(x,κ₀,τ₀,ϕ₀,κ₁,τ₁,ϕ₁))
  out_2 = dH_part(x, κ₁, τ₁, ϕ₁)
  out_3 = dH_part(-x, κ₀, τ₀, -ϕ₀)
  out_1 + log(out_2+out_3)
end

# Note: this t-cdf does not work with autodiff. Would need hand-coded logtcdf function.
function Distributions.logcdf(d::BulkAndTailsDist, x::Real)
  if (x ≤ minimum(d))
    return -Inf
  elseif (x ≥ maximum(d))
    return 0.0
  else
    (κ₀, τ₀, ϕ₀, κ₁, τ₁, ϕ₁, ν) = params(d)
    return Distributions.logcdf(TDist(ν),H(x,κ₀,τ₀,ϕ₀,κ₁,τ₁,ϕ₁))
  end
end

function Base.rand(d::BulkAndTailsDist; rng=Random.GLOBAL_RNG)
  quantile(d,rand(rng))
end

# Define R-friendly version.
batspdf(x, parms) = Distributions.pdf(BulkAndTailsDist(parms),x)
batscdf(x, parms) = Distributions.cdf(BulkAndTailsDist(parms),x)
batsquantile(x, parms) = Distributions.quantile(BulkAndTailsDist(parms),x)
batslogpdf(x, parms) = Distributions.logpdf(BulkAndTailsDist(parms),x)
batslogcdf(x, parms) = Distributions.logcdf(BulkAndTailsDist(parms),x)

# Overload for vector inputs, because R users like to not know when they are
# broadcasting functions.
batspdf(xv::Vector, parms) = [batspdf(x, parms) for x in xv]
batscdf(xv::Vector, parms) = [batscdf(x, parms) for x in xv]
batsquantile(xv::Vector, parms) = [batsquantile(x, parms) for x in xv]
batslogpdf(xv::Vector, parms) = [batslogpdf(x, parms) for x in xv]
batslogcdf(xv::Vector, parms) = [batslogcdf(x, parms) for x in xv]