text:
# Imports
import Base: start, next, done, getindex, length, show
import Base.LineEdit: KeyAlias, AnyDict

export Border, ListWidget, WidgetStack, IOBufferView, MenuWidget

const LazyWidget = Union{Widget,Signal}
value(w::Widget) = w

# Defaults
children(w::Widget) = ()
focused(w::Widget) = nothing
defocused(w::Widget) = nothing
keymap(w::Widget) = Dict{Any,Any}()
isscrollable(w::Widget) = false
scroll_up(w::Widget, it) = nothing
scroll_down(w::Widget, it) = nothing
mouse_release(w::Widget, pos) = nothing
mouse_press(w::Widget, pos) = nothing
show(io::IO, w::Widget) = print(io,typeof(w))
focusable(w::Widget) = false
optheight(w::Widget) = 1
optwidth(w::Widget) = 10

focusable{T<:Widget}(l::Signal{T}) = focusable(value(l))
optheight{T<:Widget}(l::Signal{T}) = optheight(value(l))
focus{T<:Widget}(l::Signal{T}) = focus(value(l))

# ScrollableWidget
@compat abstract type ScrollableWidget <: Widget end
isscrollable(w::ScrollableWidget) = true
cur_start_line(w::ScrollableWidget) = 1

function draw_scrollable(s::Screen, w::Widget)
    it = start_it = cur_top(w)
    clear(s)
    cur_row = 1
    done = false
    while cur_row <= height(s) && !element_done(w,it)
        nrows = draw_element_line(subscreen(s,cur_row:height(s),1:width(s)), w, it)
        cur_row += nrows
        it = element_next(w,it)
    end
    draw_scroll_indicators(s,w, up = it == element_start(w), down = !element_done(w,it))
end

function draw(s::Screen, w::Widget)
    if isvisible(w)
        if isscrollable(w)
            draw_scrollable(s, w)
        else
            draw_regular(s, w)
        end
    else
        clear(s)
    end
end

# Label
immutable Label <: Widget
    text::AbstractString
    # Doesn't get a widget context, because it's only used as an aux widget
    # by other widgets
end
optwidth(l::Label) = strwidth(l.text)

function draw_regular(s::Screen, l::Label; fg = :default, bg = :default)
    swrite(s,1,1,Base._truncate_at_width_or_chars(l.text,width(s)); fg = fg, bg = bg)
end

# SimpleSlider
type SimpleSlider{T} <: Widget
    min::T
    max::T
    val::Signal{T}
    last_sliderwidth::Int
    ctx::WidgetContext
    (::Type{SimpleSlider{T}}){T}(min::T,max::T,val::T) = new{T}(min,max,Signal(T, val),0)
end
optheight(s::SimpleSlider) = 1

# ─';'─';'⎼
function draw_regular(s::Screen,slider::SimpleSlider)
    textwidth = max(strwidth(repr(slider.min)),strwidth(repr(slider.max)))
    sliderwidth = max(0,width(s)-textwidth-1)
    fraction = ((slider.val.value - slider.min)/(slider.max - slider.min))
    fwidth = floor(Int,fraction*sliderwidth)
    if fwidth >= 1
        swrite(s,1,1:(fwidth-1),'═', fg = :red)
        swrite(s,1, fwidth+1         ,' ')
        swrite(s,1, fwidth           ,'⬤', fg = :red)
        swrite(s,1,(fwidth+2):sliderwidth,'═', fg = :default)
    end
    slider.last_sliderwidth = sliderwidth
    draw(subscreen(s,1:1,(1+sliderwidth+1):width(s)),Label(repr(slider.val.value)))
end

function slider_inc(s::SimpleSlider)
    step = (s.max-s.min)/s.last_sliderwidth
    push!(s.val,min(s.val.value + step, s.max))
    #invalidate(s)
end

function slider_dec(s::SimpleSlider)
    step = (s.max-s.min)/s.last_sliderwidth
    push!(s.val,max(s.val.value - step, s.min))
    #invalidate(s)
end

function keymap(s::SimpleSlider)
    Dict(
    # Right Arrow
    "\e[C" => (o...)->(slider_inc(s); invalidate(s)),
    # Left Arrow
    "\e[D" => (o...)->(slider_dec(s); invalidate(s))
    )
end

function mouse_release(s::SimpleSlider, pos)
    row, col = pos
    if col > s.last_sliderwidth
        push!(s.val,s.max)
    else
        push!(s.val,s.min + (s.max-s.min)*(col-1)/s.last_sliderwidth)
    end
    invalidate(s)
end
mouse_press(s::SimpleSlider, pos) = mouse_release(s, pos)
scroll_up(s::SimpleSlider, _) = slider_inc(s)
scroll_down(s::SimpleSlider, _) = slider_dec(s)

# Gauge
type Gauge{T} <: Widget
    min::T
    max::T
    val::T
    ctx::Void
    (::Type{Gauge{T}}){T}(min,max,val) = new{T}(min,max,val)
end
Gauge{T}(min::T,max::T,val::T) = Gauge{T}(min,max,val)

minheight(g::Gauge) = 1
optheight(g::Gauge) = 3

function draw_regular(s::Screen, g::Gauge)
    fraction = ((g.val - g.min)/(g.max - g.min))
    swrite(s,1:height(s),1:round(Integer,fraction*width(s)),' '; bg = :magenta)
    percentage = round(Integer,fraction*100)
    swrite(s,ceil(Integer,height(s)/2),floor(Integer,width(s)/2),"$(percentage) %";
        bg = nothing)
end

# CheckBox
type CheckBox <: Widget
    checked::Bool
    label::Label
    highlighted::Bool
    ctx::WidgetContext
    CheckBox(checked,label,highlighted) = new(checked,label,highlighted)
end
CheckBox(s::AbstractString) = CheckBox(false,Label(s),false)
focused(c::CheckBox) = c.highlighted = true
defocused(c::CheckBox) = c.highlighted = false

function draw_regular(s::Screen, c::CheckBox)
    swrite(s,1,1,c.checked ? "◉ " : "◯ ")
    draw(subscreen(s,1:height(s),3:width(s)), c.label; fg = c.highlighted ? :yellow : :default)
end

minheight(c::CheckBox) = 1

keymap(c::CheckBox) = Dict(
    ' '     => (args...)->(c.checked = !c.checked; invalidate(c))
)

# WidgetStack
type WidgetStack <: Widget
    widgets::Vector{Widget}
    curidx
    ctx::WidgetContext
    WidgetStack(widgets,curidx=1) = new(widgets,curidx,WidgetContext())
end
WidgetStack(widgets::Vector) = WidgetStack(Widget[w for w in widgets], 1)
children(w::WidgetStack) = w.widgets
focus(w::WidgetStack) = focus_child(w)

function draw_regular(s::Screen,w::WidgetStack)
    cumheight = 0
    for w in w.widgets
        height = optheight(w)
        starth = 1+cumheight
        draw(subscreen(s,starth:(starth+height-1),1:width(s), w; scroll = true),w)
        cumheight += height
    end
end

optheight(w::WidgetStack) = sum(map(optheight,w.widgets))

function update_index(w::WidgetStack,i)
    if i == w.curidx
        return
    end
    focus(w.widgets[i])
    w.curidx = i
    invalidate(w)
end

keymap(w::WidgetStack) = Dict(
    "\e[A"      => (s,p,c)->(update_index(w,mod1(w.curidx - 1,length(w.widgets)))),
    "\e[B"      => (s,p,c)->(update_index(w,mod1(w.curidx + 1,length(w.widgets)))),
    '\t'        => "\e[A"
)

# ScrollableChain
# TODO: Use the chain iterator somehow?
type ScrollableChain <: ScrollableWidget
    children::Vector{Widget}
    cur_widget::Int
    ctx::WidgetContext
    function ScrollableChain(children)
        for c in children
            if !isscrollable(c)
                error("Child of type $(typeof(c)) is not scrollable")
            end
        end
        new(children,1,WidgetContext())
    end
end
children(s::ScrollableChain) = s.children
optheight(s::ScrollableChain) = sum(map(optheight,s.children))
focusable(s::ScrollableChain) = any(map(focusable,s.children))
focus(s::ScrollableChain) = focus_child(s)

function next_valid(s::ScrollableChain, it)
    while element_done(s.children[it[1]],it[2])
        if it[1] == length(s.children)
          return (it[1]+1,)
        end
        it = (it[1]+1,element_start(s.children[it[1]+1]))
    end
    it
end

cur_top(s::ScrollableChain) = next_valid(s,(s.cur_widget,cur_top(s.children[s.cur_widget])))

function element_start(s::ScrollableChain)
    next_valid(s,(1,element_start(s.children[1])))
end
function element_next(s::ScrollableChain, idx)
    next_valid(s,(idx[1],element_next(s.children[idx[1]],idx[2])))
end
function element_prev(s::ScrollableChain, idx)
    it = (idx[1],element_prev(s.children[idx[1]],idx[2]))
    while it[2] != element_start(s.children[it[1]])
        if it[1] == 1
          return it
        end
        it = (it[1]-1,element_start(s.children[it[1]-1]))
    end
    it
end
element_done(s::ScrollableChain, idx) = idx[1] > length(s.children)

function set_cur_top(s::ScrollableChain, it)
    s.cur_widget = it[1]
    set_cur_top(s.children[it[1]],it[2])
end

draw_element_line(s::Screen, c::ScrollableChain, it) =
    isvisible(c.children[it[1]]) && draw_element_line(s, c.children[it[1]], it[2])

function scroll_up(s::ScrollableChain, pos)
    cur = cur_top(s)
    if element_done(s, cur)
        return
    end
    n = element_next(s, cur)
    set_cur_top(s, n)
end
function scroll_down(s::ScrollableChain, pos)
    cur = cur_top(s)
    if element_start(s) == cur
        return
    end
    n = element_prev(s, cur)
    set_cur_top(s, n)
end

function currently_focused(w::ScrollableChain)
    # Search the the focus stack for our widget
    i = findfirst(w.ctx.focuss.stack,w)
    @assert i != 0
    if i == length(w.ctx.focuss.stack)
        return (0,w)
    else
        child = w.ctx.focuss.stack[i+1]
        idx = findfirst(w.children,child)
        @assert idx != 0
        return (idx,child)
    end
end
function focus_next(w::ScrollableChain; show_invisible = false)
    # show_visible = false, should skip over hidden, not implemented
    next = currently_focused(w)[1]+1
    if next > length(w.children)
        return false
    end
    if show_invisible
        makevisible(w.children[next])
    end
    focus(w.children[next])
    invalidate(w)
    return true
end

keymap(w::ScrollableChain) = Dict(
    '\t' => (args...)->focus_next(w)
)

child_ok(w::Widget) = true
child_ok{T}(w::Signal{T}) = T <: Widget

init_child(w::Widget,this::Widget) = this
function init_child{T<:Widget}(child::Signal{T},this::Widget)
    lift(child) do w
        w.ctx.parent = this
        if isdefined(this.ctx,:focuss)
            w.ctx.focuss = this.ctx.focuss
        end
        invalidate(this)
        nothing
    end
    return this
end

# Border
type Border <: Widget
    child::LazyWidget
    label::Label
    ctx::WidgetContext
    function Border(child::LazyWidget,label::Label)
        @assert child_ok(child)
        init_child(child,new(child,label,WidgetContext()))
    end
end
# first one is deprecated
Border(child::LazyWidget,label::AbstractString="") = Border(child,Label(label))
Border(label::AbstractString,child::LazyWidget) = Border(child,Label(label))
children(b::Border) = (b.child,)
focusable(b::Border) = focusable(b.child)
focus(b::Border) = focus(b.child)

minheight(b::Border) = 2
optheight(b::Border) = 2 + optheight(b.child)
optwidth(b::Border) = max(4 + optwidth(b.label), optwidth(b.child))

function draw_regular(s::Screen,b::Border)
    color = isfocused(b) ? :green : :default
    swrite(s,1              ,1             ,'┌', attrs = IsACS, fg = color) # Top Left
    swrite(s,1              ,2:(width(s)-1),'─', attrs = IsACS, fg = color) # Top
    swrite(s,1              ,   width(s)   ,'┐', attrs = IsACS, fg = color) # Top Right
    swrite(s,2:(height(s)-1),1             ,'│', attrs = IsACS, fg = color) # Left
    swrite(s,2:(height(s)-1),   width(s)   ,'│', attrs = IsACS, fg = color) # Right
    swrite(s,   height(s)   ,1             ,'└', attrs = IsACS, fg = color) # Bottom Left
    swrite(s,  (height(s)  ),2:(width(s)-1),'─', attrs = IsACS, fg = color) # Bot
    swrite(s,   height(s)   ,   width(s)   ,'┘', attrs = IsACS, fg = color) # Bottom Right
    # Label
    opt = optwidth(b.label)
    draw(subscreen(s,1:1,3:(min(2+opt,width(s)-1))),b.label)
    # Child
    draw(subscreen(s,2:(height(s)-1),2:(width(s)-1), value(b.child); scroll = true),value(b.child))
end

# TightCentering
# Centers the widget to the size given by its `minbounds`
type TightCentering <: Widget
    w::Widget
    ctx::WidgetContext
    TightCentering(w) = new(w,WidgetContext())
end
children(t::TightCentering) = (t.w,)
focus(t::TightCentering) = focus(t.w)

function draw_regular(s::Screen, w::TightCentering)
    cw = w.w
    bounds = (optheight(cw),optwidth(cw))
    starth = div(height(s)-bounds[1],2)
    startw = div(width(s)-bounds[2],2)
    draw(subscreen(s,starth:(starth+bounds[1]),startw:(startw+bounds[2])),cw)
end

# Layout
type RowLayout <: Widget
    cols::Vector{Widget}
    spans::Vector{Int}
    fixedwidths::Vector{Int}
    span::Int
    ctx::WidgetContext
    function RowLayout(cols,spans,fixedwidths,span)
        new(cols,spans,fixedwidths,span,WidgetContext())
    end
end
children(r::RowLayout) = r.cols
cur_top(r::RowLayout) = map(cur_top,r.cols)
focusable(r::RowLayout) = any(map(focusable,r.cols))
focus(r::RowLayout) = focus_child(r)

element_start(r::RowLayout) = map(element_start,r.cols)
function element_next(r::RowLayout,i)
    map(zip(i,r.cols)) do x
        j, w = x
        if !element_done(w,j)
            return element_next(w,j)
        else
            return j
        end
    end
end
function element_done(r::RowLayout,i)
    all(zip(i,r.cols)) do x
        j, w = x
        element_done(w,j)::Bool
    end
end

# A row layout is scrollable if all its children are
isscrollable(r::RowLayout) = all(map(isscrollable,children(r)))

function optheight(rl::RowLayout)
    maximum([optheight(w) for w in rl.cols])
end

function compute_widths(rl::RowLayout,totalwidth)
    # The extra width will be given to the last column
    curwidth = 0
    numfixed = 0
    totalfixed = 0
    for i in eachindex(rl.fixedwidths)
        if rl.fixedwidths[i] != 0
            numfixed += rl.spans[i]
            totalfixed += rl.fixedwidths[i]
        end
    end
    widths = Array{UnitRange}(0)
    sizehint!(widths,length(rl.cols))
    for i = eachindex(rl.cols)
        c = rl.cols[i]
        span = rl.spans[i]
        fixed = rl.fixedwidths[i]
        if fixed == 0
            w = floor(Integer,span*(totalwidth-totalfixed)/(rl.span-numfixed))
        else
            w = fixed
        end
        startw = 1+curwidth
        curwidth += w
        if i == endof(rl.cols)
            endw = totalwidth
        else
            endw = startw + w - 1
        end
        push!(widths,startw:endw)
    end
    widths
end

function draw_regular(s::Screen,rl::RowLayout)
    for (w,r) in zip(rl.cols,compute_widths(rl,width(s)))
        draw(subscreen(s,1:height(s),r, w; scroll = true),w)
    end
end

function draw_element_line(s::Screen,rl::RowLayout,it)
    for (w,r,i) in zip(rl.cols,compute_widths(rl,width(s)),it)
        r = draw_element_line(subscreen(s,1:height(s),r),w,i)::Int
        @assert r == 1
    end
    return 1
end

# CellDisplay
# Just a bunch of terminal cells
type CellDisplay <: ScrollableWidget
    rows::Vector{Vector{Cell}}
    cur_top::Int
    ctx::WidgetContext
    CellDisplay(rows) = new(rows,1,WidgetContext())
end
optheight(cd::CellDisplay) = length(cd.rows)
cur_top(cd::CellDisplay) = cd.cur_top
element_start(cd::CellDisplay) = 1
element_next(cd::CellDisplay, i) = i+1
element_done(cd::CellDisplay, i) = i > length(cd.rows)

function draw_element_line(s::Screen, d::CellDisplay, row)
    ncells = min(width(s),length(d.rows[row]))
    for col in 1:ncells
        c = d.rows[row][col]
        swrite(s, 1, col, c)
    end
    return 1
end

function draw_regular(s::Screen, d::CellDisplay)
    for row in 1:min(height(s),length(d.rows))
        draw_element_line(s, d, row)
    end
end

#IOBufferView
type IOBufferView <: ScrollableWidget
    buf::IOBuffer
    do_show_cursor::Bool
    may_show_cursor::Bool
    draw_empty_last_line::Bool
    cur_top::UnitRange
    ctx::WidgetContext
    function IOBufferView(buf,may_show_cursor=true,draw_empty_last_line=true)
        this = new(buf,false,may_show_cursor,draw_empty_last_line,0:0,WidgetContext())
        this.cur_top = element_start(this)
        this
    end
end
optheight(view::IOBufferView) = 10
cur_top(view::IOBufferView) = view.cur_top
set_cur_top(view::IOBufferView, it::UnitRange) = view.cur_top = it
focusable(view::IOBufferView) = true

function scroll_up(l::IOBufferView, pos)
    it = element_next(l, l.cur_top)
    if !element_done(l, it)
        l.cur_top = it
        invalidate(l)
    end
end
function scroll_down(l::IOBufferView, pos)
    if l.cur_top != element_start(l)
        l.cur_top = element_prev(l, l.cur_top)
        invalidate(l)
    end
end

function didwrite(view::IOBufferView,at,nmoved)
    at = at+1
    #@show (at,view.cur_top)
    if at < first(view.cur_top)
        # Need to move the start of the range by the number of bytes inserted
        view.cur_top = (first(view.cur_top)+nmoved):(last(view.cur_top)+nmoved)
    elseif first(view.cur_top) <= at && (
        # If we inserted into the current line
        at <= last(view.cur_top) ||
        # Or the buffer was empty
        last(view.cur_top) == -1 ||
        # Or this is the last line
        sizeof(view.buf.data)-nmoved == last(view.cur_top))
        # Naively, we'd just have to add nmoved to the end, but we may have
        # writeen one or more newlines, so we need to look for that
        idx = findnext(view.buf.data,UInt8('\n'),at)
        start = first(view.cur_top)
        view.cur_top = (start == 0 ? 1 : start):(idx == 0 ? endof(view.buf.data) : idx-1)
    end
    #@show view.cur_top
end

buffer(view::IOBufferView) = view.buf
focused(w::IOBufferView) = w.may_show_cursor && (w.do_show_cursor = true)
defocused(w::IOBufferView) = w.may_show_cursor && (w.do_show_cursor = false)

# 0:-1 is used as the invalid marker, but next is arranged in such
# a way that next of invalid is the first valid line. Since done is called
# before next in normal application code, this is not a problem.
element_start(view::IOBufferView) = element_next(view,0:-1)
function element_next(view::IOBufferView, cur_range)
    # Special case: An empty buffer with draw_empty_last_line
    # will show an empty line with a cursor. This state is represented as 1:0
    if view.draw_empty_last_line && isempty(view.buf.data) && cur_range == 0:-1
        return 1:0
    end
    if cur_range == 1:0 || last(cur_range) == sizeof(view.buf.data) ||
        isempty(view.buf.data)
        return 0:-1
    end
    # Cur range will have ended on the character before a \n, so the point
    # to start searching is the character after the \n (i.e. +2)
    idx = findnext(view.buf.data,UInt8('\n'),last(cur_range)+2)
    if idx == 0
        # The last line is empty
        if !view.draw_empty_last_line &&
            (last(cur_range)+2 >= sizeof(view.buf.data))
            return 0:-1
        else
            return  (last(cur_range)+2):sizeof(view.buf.data)
        end
    else
        return (last(cur_range)+2):(idx-1)
    end
end
function element_prev(view::IOBufferView, cur_range)
    idx = findprev(view.buf.data,'\n',first(cur_range)-2)
    if idx == 0
        return element_start(view)
    end
    return (idx+1):(first(cur_range)-2)
end
element_done(view::IOBufferView,range) = range == 0:-1

# This implementation is significantly simpler than the LineEdit.jl
# implementation, because it doesn't have to deal with a)
#   A) Calculating movement commands
#   B) Worry About the Prompt
#   C) Position the actual cursor correctly
# `offsets` are offsets into the IOBuffer
function draw_element_line(s::Screen, input::IOBufferView, offsets)
    cols = width(s)
    buf = buffer(input)
    text = @compat(String)(buf.data)
    buf_pos = position(buf)
    curs_pos = buf_pos+1
    row = 1
    col = 1
    offset = first(offsets)
    justwrapped = false
    em = LineEmulator(Cell(Cell('\0'),fg=9,bg=9))
    seek(buf,first(offsets)-1)
    while position(buf)+1 <= last(offsets) && text[position(buf)+1] != '\n'
        iscurs = (offset == curs_pos) && input.do_show_cursor
        #@assert !iscurs
        c = parse_cell!(em,buf)
        if (c.content == '\0')
            break
        end
        # For now, just consider a tabstop every 4 cells
        if c.content == '\t'
            nexttab = 4*(div(col-1,4)+1)+1
            swrite(s, row, col:(nexttab-1), Cell(c,content=' '))
            col = nexttab
        else
            swrite(s, row, col, Cell(Cell(c),
                fg = iscurs ? :black : :default,
                bg = iscurs ? :white : :default))
            col += 1
        end
        if col == cols+1
            row = row + 1
            col = 1
            justwrapped = true
        end
        offset = nextind(text,offset)
    end
    seek(buffer(input),buf_pos)
    lastoffset = offset
    ((offset <= last(offsets) && text[offset] == '\n') ||
        (offset > endof(text))) && (offset += 1)
    # Render the cursor if it's after the end of the text
    if curs_pos == lastoffset && input.do_show_cursor
        swrite(s, row, col, ' ', bg = :blue)
    end
    return row
end

function draw_scroll_indicators(s::Screen, input; up = true, down = true)
    up   && swrite(s,1        , width(s), '▲', fg=:blue)
    down && swrite(s,height(s), width(s), '▼', fg=:blue)
end

# TextInput
# Base on LineEdit.jl's multi line input widget
import Base.LineEdit: edit_insert, edit_backspace,
    edit_move_right, edit_move_left, move_input_end, move_line_end

type TextInput <: ScrollableWidget
    view::IOBufferView
    on_enter::Function
    on_done::Function
    multiline::Bool
    tabcomplete::Bool
    optheight::Int
    ctx::WidgetContext
    TextInput(buf=IOBuffer(),multiline=false,tabcomplete=false,draw_empty_last_line=true) =
        new(IOBufferView(buf,true,draw_empty_last_line),
            (args...)->true,(args...)->nothing,tabcomplete,multiline,10,
            WidgetContext())
end

children(i::TextInput) = (i.view,)
optheight(e::TextInput) = e.optheight
draw(s::Screen, input::TextInput) = draw(s,input.view)
focusable(w::TextInput) = true
focus(w::TextInput) = focus(w.view)

element_start(w::TextInput)     = element_start(w.view)
element_next(w::TextInput, idx) = element_next(w.view, idx)
element_done(w::TextInput, idx) = element_done(w.view, idx)
cur_top(w::TextInput)           = cur_top(w.view)
draw_element_line(s::Screen, w::TextInput, idx) =
    draw_element_line(s, w.view, idx)

# LineEdit support
buffer(w::TextInput) = w.view.buf
function edit_insert(w::TextInput,args...)
    pos = position(w.view.buf)
    n = edit_insert(buffer(w),args...)
    didwrite(w.view,pos,n)
    invalidate(w)
end
function edit_backspace(w::TextInput)
    pos = position(w.view.buf)
    edit_backspace(buffer(w))
    didwrite(w.view, pos, -1)
    invalidate(w)
end
edit_move_left(w::TextInput) = edit_move_left(buffer(w)) && invalidate(w)
edit_move_right(w::TextInput) = edit_move_right(buffer(w)) && invalidate(w)
move_input_end(w::TextInput) = (move_input_end(buffer(w)); invalidate(w))
move_line_end(w::TextInput) = (move_line_end(buffer(w)); invalidate(w))

function keymap(w::TextInput)
    d = AnyDict(
    # Enter
    '\r' => (s,o...)->begin
        line = String(take!(copy(buffer(w))))
        if w.on_enter(line)
            w.on_done(line)
        end
    end,
    '\n' => KeyAlias('\r'),
    # Backspace/^H
    '\b' => (s,o...)->edit_backspace(w),
    127 => KeyAlias('\b'),
    # Meta Backspace
    "\e\b" => (s,o...)->edit_delete_prev_word(w),
    # ^D
    "^D" => (s,o...)->begin
        if buffer(s).size > 0
            edit_delete(w)
        else
            return :abort
        end
    end,
    "^B" => (s,o...)->edit_move_left(w),
    "^F" => (s,o...)->edit_move_right(w),
    # Meta B
    "\eb" => (s,o...)->edit_move_word_left(w),
    # Meta F
    "\ef" => (s,o...)->edit_move_word_right(w),
    # Meta Enter
    "\e\r" => (s,o...)->(edit_insert(w, '\n')),
    "\e\n" => "\e\r",
    # Simply insert it into the buffer by default
    "*" => (s,data,c)->(edit_insert(w, c)),
    "^U" => (s,o...)->edit_clear(w),
    "^K" => (s,o...)->edit_kill_line(w),
    "^Y" => (s,o...)->edit_yank(w),
    "^A" => (s,o...)->(move_line_start(w); invalidate(w)),
    "^E" => (s,o...)->(move_line_end(w); invalidate(w)),
    # Try to catch all Home/End keys
    "\e[H"  => (s,o...)->(move_input_start(w); invalidate(w)),
    "\e[F"  => (s,o...)->(move_input_end(w); invalidate(w)),
    #"^L" => (s,o...)->(Terminals.clear(terminal(s)); refresh_line(s)),
    "^W" => (s,o...)->edit_werase(w),
    # Meta D
    "\ed" => (s,o...)->edit_delete_next_word(w),
    "^C" => (s,o...)->begin
        try # raise the debugger if present
            ccall(:jl_raise_debugger, Int, ())
        end
        move_input_end(w)
        invalidate(w)
        #transition(s, :reset)
        #refresh_line(s)
        return :done
    end,
    "^Z" => (s,o...)->(return :suspend),
    # Right Arrow
    "\e[C" => (s,o...)->edit_move_right(w),
    # Left Arrow
    "\e[D" => (s,o...)->edit_move_left(w),
    # Delete
    "\e[3~" => (s,o...)->edit_delete(w),
    "^T" => (s,o...)->edit_transpose(w),
)
    if w.tabcomplete
        d['\t'] = (s,o...)->begin
            buf = buffer(w)
            # 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(w, " "^4)
                    return
                end
            end
            complete_line(w)
            invalidate(w)
        end
    end

    if w.multiline
        merge!(d,AnyDict(
            # Up Arrow
            "\e[A" => (s,o...)->edit_move_up(w),
            # Down Arrow
            "\e[B" => (s,o...)->edit_move_down(w),
        ))
    end

    d
end

Base.Terminals.beep(t::TextInput) = nothing
Base.println(t::TextInput) = error()

# REPLWidget
import Base: LineEdit
import Base.LineEdit: ModalInterface
import Base.REPL: AbstractREPL

if isdefined(Base.REPL,:MIRepl)
    import Base.REPL.MIRepl
else
    const MIRepl = Base.REPL.LineEditREPL
end

#=
type REPLWidget
    repl::ModalInterface
    widget::Widget
    ctx::WidgetContext
    REPLWidget(repl,widget) = new(repl,widget)
end=#

type WidgetREPL <: AbstractREPL
    widget::Widget
    inputarea::Widget
    historybuf::IOBuffer
    historyview::IOBufferView
    mi::MIRepl
    WidgetREPL() = new()
end
Base.REPL.reset(repl::WidgetREPL) = nothing
Base.REPL.backend(repl::WidgetREPL) = Base.REPL.backend(repl.mi)
Base.LineEdit.deactivate(a,s,c, ::TerminalUI.TextInput) = s
Base.LineEdit.activate(a,b,c, ::TerminalUI.TextInput) = nothing
Base.LineEdit.commit_changes(::TerminalUI.TextInput, _) = nothing
function Base.REPL.prepare_next(repl::WidgetREPL)
    makevisible(repl.inputarea)
end

immutable WidgetREPLDisplay
    repl::WidgetREPL
end

function write_history(repl::WidgetREPL, tempbuf::IOBuffer)
    pos = position(repl.historybuf)+1
    didwrite(repl.historyview,pos,write(repl.historybuf,take!(tempbuf)))
end

function display(d::WidgetREPLDisplay, ::MIME"text/plain", x)
    statebuf = d.repl.historybuf
    tempbuf = IOBuffer()
    writemime(tempbuf, MIME("text/plain"), x)
    println(tempbuf)
    write_history(d.repl,tempbuf)
end
display(d::WidgetREPLDisplay, x) = display(d, MIME("text/plain"), x)

function Base.REPL.print_response(repl::WidgetREPL, val::ANY, bt, show_value::Bool, have_color::Bool)
    while true
        try
            if bt !== nothing
                tempbuf = IOBuffer()
                Base.REPL.display_error(tempbuf, val, bt)
                println(tempbuf)
                write_history(repl,tempbuf)
                iserr, lasterr = false, ()
            else
                if val !== nothing && show_value
                    try
                        display(WidgetREPLDisplay(repl),val)
                    catch err
                        println(repl.historybuf, "Error showing value of type ", typeof(val), ":")
                        rethrow(err)
                    end
                end
            end
            break
        catch err
            if bt !== nothing
                tempbuf = IOBuffer()
                println(tempbuf, "SYSTEM: show(lasterr) caused an error")
                println(tempbuf, err)
                write_history(repl,tempbuf)
                break
            end
            val = err
            bt = catch_backtrace()
        end
    end
end

function WidgetForREPL(repl::ModalInterface,wr::WidgetREPL)
    inputbuf = IOBuffer(b"", true, true)
    wr.historybuf = IOBuffer(b"",true,true)
    wr.historyview = IOBufferView(wr.historybuf,false,false)
    input = TextInput(inputbuf,true,true,false)

    mistate = LineEdit.init_state(input, repl)

    lift(input.content) do line
        if LineEdit.on_enter(mistate) && line != ""
            tempbuf = IOBuffer()
            write(tempbuf,"\033[1m\033[32mjulia> \033[1m\033[37m",line,"\033[0m\n")
            write_history(wr, tempbuf)
            buf = input.view.buf
            b = IOBuffer()
            input.view.buf = b
            input.view.cur_top = element_start(input.view)
            hide(wr.inputarea)
            isdefined(Main,:redraw) && Main.redraw()
            mistate.mode_state[mistate.current_mode].p.on_done(mistate,buf,true)
            isdefined(Main,:redraw) && Main.redraw()
            nothing
        end
    end

    cd = TerminalUI.CellDisplay(
        Vector{VT100.Cell}[[
            VT100.Cell(TerminalUI.dummycell(:green,:default,0),content = c)
            for c in "julia> "]])

    wr.inputarea = RowLayout([cd,input],[1,1],[7,0],2)
    b = ScrollableChain([wr.historyview,wr.inputarea])
end

import Base.LineEdit: refresh_multi_line
function refresh_multi_line(termbuf::LineEdit.TerminalBuffer,
    w::TextInput, s::LineEdit.ModeState)
    #invalidate(w)
end
function refresh_multi_line(w::TextInput, args...)
    #invalidate(w)
end

# ImageWidget
type ImageWidget <: Widget
    data::Vector{UInt8}
    optheight::Int
    # Screen char if allocated
    screenchar::Cell
    ctx::WidgetContext
    ImageWidget() = new(Vector{UInt8}(),10,Cell(Char(0)))
end
optheight(w::ImageWidget) = w.optheight

function draw(s::Screen,img::ImageWidget)
    clear(s)
    if img.screenchar.content != Char(0)
        swrite_image(s, 1:height(s), 1:width(s), img.screenchar)
    elseif !isempty(img.data)
        img.screenchar = swrite_image(s, 1:height(s), 1:width(s), img.data)
    end
end

# ListWidget
type ListWidget <: Widget
    item::Any
    draw_indicies::Bool
    highlight_idx::Int
    highlighted::Signal{Any}
    cur_top
    ctx::WidgetContext
    function ListWidget(item)
        this = new(item,true,start(item),Signal(Any, start(item)),start(item),WidgetContext())
        map(this.highlighted) do it
            invalidate(this)
        end
        this
    end
end
optheight(l::ListWidget) = isa(Base.iteratorsize(l.item),Union{Base.HasLength,Base.HasShape}) ? length(l.item) : 10
cur_top(l::ListWidget) = l.cur_top
isscrollable(l::ListWidget) = true

element_start(l::ListWidget) = start(l.item)
element_next(l::ListWidget,it) = next(l.item,it)[2]
element_prev(l::ListWidget,it) = prevind(l.item,it)
element_done(l::ListWidget,it) = done(l.item,it)
function draw_element_line(s::Screen, l::ListWidget, it)
    if l.draw_indicies
        # Try to figure out how much space we'll need for writing the index
        lidx = endof(l.item)
        maxwidth = strwidth(sprint(print,lidx))
        swrite(s, 1, 1, '[')
        swrite(s, 1, 2 + maxwidth, ']')
        draw(subscreen(s,1:height(s),2:(2+maxwidth)),Label(sprint(print,it)))
        s = subscreen(s,1:height(s),(4+maxwidth):width(s))
    end
    highlight = l.highlight_idx == it
    clear(s)
    if highlight
        swrite(s, 1:height(s), 1:width(s), ' ', fg = :black, bg = :yellow)
    end
    draw_regular(s,Label(sprint(print,l.item[it]));
        fg = highlight ? nothing : :default, bg = highlight ? nothing : :default)
    1
end

function scroll_up(l::ListWidget, pos)
    it = element_next(l, l.cur_top)
    if !element_done(l, it)
        l.cur_top = it
        invalidate(l)
    end
end
function scroll_down(l::ListWidget, pos)
    if l.cur_top != element_start(l)
        l.cur_top = element_prev(l, l.cur_top)
        invalidate(l)
    end
end

keymap(l::ListWidget) = Dict(
    "\e[A"      => (s,p,c)->begin
        if l.highlight_idx != element_start(l)
            l.highlight_idx = element_prev(l, l.highlight_idx)
            push!(l.highlighted, l.highlight_idx)
            invalidate(l)
        end
    end,
    "\e[B"      => (s,p,c)->begin
        it = element_next(l, l.highlight_idx)
        if !element_done(l, it)
            l.highlight_idx = it
            push!(l.highlighted, it)
            invalidate(l)
        end
    end
)

# Like a list Widget, but one of the items can be selected (has an on_done callback),
# and selections can be made by numeric entry
type MenuWidget <: Widget
    on_done::Any
    list_widget::ListWidget
    ctx::WidgetContext
    function MenuWidget(on_done, item)
        new(on_done, ListWidget(item),WidgetContext())
    end
end
children(w::MenuWidget) = (w.list_widget,)
focusable(w::MenuWidget) = true
focus(w::MenuWidget) = focus(w.list_widget)
isscrollable(W::MenuWidget) = true
cur_top(w::MenuWidget) = cur_top(w.list_widget)
optheight(w::MenuWidget) = optheight(w.list_widget)

element_start(l::MenuWidget) = element_start(l.list_widget)
element_next(l::MenuWidget,it) = element_next(l.list_widget,it)
element_prev(l::MenuWidget,it) = element_prev(l.list_widget,it)
element_done(l::MenuWidget,it) = element_done(l.list_widget,it)

draw_element_line(s::Screen, l::MenuWidget, it) =
    draw_element_line(s, l.list_widget, it)

keymap(widget::MenuWidget) = Dict(
    '\r' => (s,p,c) -> widget.on_done(widget.list_widget.highlight_idx)
)


text:
function subtree_size(g::DirectedGraph, n::Integer,
                      visited=falses(nodecount(g)),
                      stack=zeros(Int, nodecount(g)))
    first, last = 1, 2
    stack[first] = n
    visited[n] = true
    @inbounds while first < last
        i = stack[first]
        @inbounds for j = inneighbors(g, i)
            if !visited[j]
                stack[last] = Int(j)
                last += 1
                visited[j] = true
            end
        end
        first += 1
    end
    @inbounds for i = 1:last - 1
        visited[stack[i]] = false
    end
    return first - 2
end

"""
    subtree_size(g::DirectedGraph[, n::Integer])

Number of nodes that can reach node `n` following directed links.
If `n` not provided, computes this number for each node.
"""
function subtree_size(g::DirectedGraph)
    vis = falses(nodecount(g))
    stack = zeros(Int, nodecount(g))
    res = zeros(Float64, nodecount(g))
    @inbounds for i = nodes(g)
        res[i] = Float64(subtree_size(g, i, vis, stack))
    end
    return res
end

ext: jl


text:
import NaiveNASlib:CompGraph, CompVertex, InputVertex, SimpleDiGraph
import LightGraphs:adjacency_matrix,is_cyclic

@testset "Computation graph tests" begin

    @testset "Scalar computation graphs" begin

        # Setup a simple scalar graph which sums two numbers
        ins = InputVertex.(1:3)
        sumvert = CompVertex(+, ins[1], ins[2])
        scalevert = CompVertex(x -> 2x, sumvert)
        graph = CompGraph(inputs(sumvert), [scalevert])
        sumvert2 = CompVertex((x,y) -> x+y+2, ins[1], ins[3])
        graph2out = CompGraph(ins, [scalevert, sumvert2])

        @testset "Structural tests" begin
            @test adjacency_matrix(SimpleDiGraph(graph)) == [0 0 1 0; 0 0 1 0; 0 0 0 1; 0 0 0 0]
            @test adjacency_matrix(SimpleDiGraph(graph2out)) == [0 0 1 0 0 1;
            0 0 1 0 0 0; 0 0 0 1 0 0; 0 0 0 0 0 0; 0 0 0 0 0 1; 0 0 0 0 0 0]

            @test nv(graph) == 4
            @test nv(graph2out) == 6

            @test vertices(graph) == [ins[1], ins[2], sumvert, scalevert]
            @test vertices(graph2out) == [ins[1], ins[2], sumvert, scalevert, ins[3], sumvert2]
        end

        @testset "Computation tests" begin
            @test graph(2,3) == 10
            @test graph([2], [3]) == [10]
            @test graph(0.5, 1.3) ≈ 3.6
            @test graph2out(4,5,8) == (18, 14)
        end

    end

    @testset "Array computation graphs" begin

        # Setup a graph which scales one of the inputs by 3 and then merges is with the other
        ins = InputVertex.(1:2)
        scalevert = CompVertex(x -> 3 .* x, ins[1])
        mergegraph1 = CompGraph(ins, [CompVertex(hcat, ins[2], scalevert)])
        mergegraph2 = CompGraph(ins, [CompVertex(hcat, scalevert, ins[2])])

        @testset "Computation tests" begin
            @test CompGraph([ins[1]], [scalevert])(ones(Int64, 1, 2)) == [3 3]

            @test mergegraph1(ones(Int64, 3,2), ones(Int64, 3,3)) == [1 1 1 3 3; 1 1 1 3 3; 1 1 1 3 3]
            @test mergegraph2(ones(Int64, 3,2), ones(Int64, 3,3)) == [3 3 1 1 1; 3 3 1 1 1; 3 3 1 1 1]
        end

    end

    @testset "Simple graph copy" begin
        ins = InputVertex.(1:3)
        v1 = CompVertex(+, ins[1], ins[2])
        v2 = CompVertex(vcat, v1, ins[3])
        v3 = CompVertex(vcat, ins[1], v1)
        v4 = CompVertex(-, v3, v2)
        v5 = CompVertex(/, ins[1], v1)
        graph = CompGraph(ins, [v5, v4])

        gcopy = copy(graph)

        @test issame(graph, gcopy)
        @test graph(3,4,10) == gcopy(3,4,10)
    end

    @testset "Mutation graph copy" begin
        ins = InputSizeVertex.(InputVertex.(1:3), 1)
        v1 = ins[1] + ins[2]
        v2 = conc(v1, ins[3], dims=1)
        v3 = conc(ins[1], v1, dims=1)
        v4 = v3 - v2
        v5 = ins[1] / v1
        v6 = absorbvertex(identity, 4, v5)
        graph = CompGraph(ins, [v4, v6])

        gcopy = copy(graph)

        @test issame(graph, gcopy)
        @test graph(3,4,10) == gcopy(3,4,10)

        newop(x...; cf=clone) = clone(x...; cf=cf)
        newop(v::MutationVertex, ins...; cf=clone) = MutationVertex(cf(v.base, ins..., cf=cf), newop(trait(v), v), cf(v.trait, cf=cf))
        newop(::MutationTrait, v::MutationVertex) = clone(op(v))
        newop(::SizeAbsorb, v::MutationVertex) = IoIndices(nin(v), nout(v))
        graph_inds = copy(graph, newop)

        @test !issame(graph_inds, graph)
        @test graph(3,4,10) == graph_inds(3,4,10)

        # Nothing should have changed with original
        function testop(v) end
        testop(v::MutationVertex) = testop(trait(v), v)
        function testop(::MutationTrait, v) end
        testop(::SizeAbsorb, v) = @test typeof(op(v)) == IoChange
        foreach(testop, mapreduce(flatten, vcat, graph.outputs))

        # But new graph shall use IoIndices
        testop(::SizeAbsorb, v) = @test typeof(op(v)) == IoIndices
        foreach(testop, mapreduce(flatten, vcat, graph_inds.outputs))
    end

    @testset "Graph add trait" begin
        struct MockTrait <: DecoratingTrait
            t::MutationTrait
        end

        inver = inputvertex("in", 3)
        v1 = absorbvertex(+, 5, inver)
        v2 = conc(inver, v1, dims=1)
        graph = CompGraph(inver, v2)

        addtrait(x...;cf) = clone(x...;cf=cf)
        addtrait(t::MutationTrait; cf=cf) = MockTrait(clone(t, cf=clone))

        graphnew = copy(graph, addtrait)

        function testvert(::InputSizeVertex) end
        testvert(v) = @test v.trait isa MockTrait

        foreach(testvert, vertices(graphnew))
    end

    @testset "Graph rename" begin
        v0 = inputvertex("in", 3)
        v1 = absorbvertex(+, 5, v0, traitdecoration = t -> NamedTrait(t, "v1"))
        v2 = conc(v0, v1, dims=1, traitdecoration = t -> NamedTrait(t, "v2"))
        graph = CompGraph(v0, v2)

        rename(x...;cf) = clone(x...;cf=cf)
        rename(s::String; cf) = s * "new"

        graphnew = copy(graph, rename)

        @test name.(vertices(graphnew)) == ["innew", "v1new", "v2new"]
    end

    @testset "Topological sort" begin
        in1,in2 = InputVertex.(("in1", "in2"))
        v1 = CompVertex(+, in1, in2)
        v2 = CompVertex(-, in2, in1)
        v3 = CompVertex(vcat, v1,v2,v1)
        g = CompGraph([in1, in2], v3)

        @test vertices(g) == [in1, in2, v1, v2, v3]
    end
end


text:
depsjl = joinpath(dirname(@__FILE__), "../..", "deps", "deps.jl")
include(depsjl)

function ctriangulate(inTri::TriangulateIO, options::String)
  outTri = TriangulateIO()
  voronoiTri = TriangulateIO()

  ccall(
    (:call_triangulate, libtriangle),
    Nothing,
    (Ptr{UInt8}, Ref{TriangulateIO}, Ref{TriangulateIO}, Ref{TriangulateIO}),
    options, Ref(inTri), Ref(outTri), Ref(voronoiTri)
  )

  (outTri, voronoiTri)
end


text:
# ------------------------------------------------------------------------------------------
# # Introduction to DataFrames
# **[Bogumił Kamiński](http://bogumilkaminski.pl/about/), May 23, 2018**
# ------------------------------------------------------------------------------------------

using DataFrames # load package

# ------------------------------------------------------------------------------------------
# ## Load and save DataFrames
# We do not cover all features of the packages. Please refer to their documentation to learn
# them.
#
# Here we'll load `CSV` to read and write CSV files and `JLD`, which allows us to work with
# a Julia native binary format.
# ------------------------------------------------------------------------------------------

using CSV
using JLD

# ------------------------------------------------------------------------------------------
# Let's create a simple `DataFrame` for testing purposes,
# ------------------------------------------------------------------------------------------

x = DataFrame(A=[true, false, true], B=[1, 2, missing],
              C=[missing, "b", "c"], D=['a', missing, 'c'])

# ------------------------------------------------------------------------------------------
# and use `eltypes` to look at the columnwise types.
# ------------------------------------------------------------------------------------------

eltypes(x)

# ------------------------------------------------------------------------------------------
# Let's use `CSV` to save `x` to disk; make sure `x.csv` does not conflict with some file in
# your working directory.
# ------------------------------------------------------------------------------------------

CSV.write("x.csv", x)

# ------------------------------------------------------------------------------------------
# Now we can see how it was saved by reading `x.csv`.
# ------------------------------------------------------------------------------------------

print(read("x.csv", String))

# ------------------------------------------------------------------------------------------
# We can also load it back. `use_mmap=false` disables memory mapping so that on Windows the
# file can be deleted in the same session.
# ------------------------------------------------------------------------------------------

y = CSV.read("x.csv", use_mmap=false)

# ------------------------------------------------------------------------------------------
# When loading in a `DataFrame` from a `CSV`, all columns allow `Missing` by default. Note
# that the column types have changed!
# ------------------------------------------------------------------------------------------

eltypes(y)

# ------------------------------------------------------------------------------------------
# Now let's save `x` to a file in a binary format; make sure that `x.jld` does not exist in
# your working directory.
# ------------------------------------------------------------------------------------------

save("x.jld", "x", x)

# ------------------------------------------------------------------------------------------
# After loading in `x.jld` as `y`, `y` is identical to `x`.
# ------------------------------------------------------------------------------------------

y = load("x.jld", "x")

# ------------------------------------------------------------------------------------------
# Note that the column types of `y` are the same as those of `x`!
# ------------------------------------------------------------------------------------------

eltypes(y)

# ------------------------------------------------------------------------------------------
# Next, we'll create the files `bigdf.csv` and `bigdf.jld`, so be careful that you don't
# already have these files on disc!
#
# In particular, we'll time how long it takes us to write a `DataFrame` with 10^3 rows and
# 10^5 columns to `.csv` and `.jld` files.  *You can expect JLD to be faster!* Use
# `compress=true` to reduce file sizes.
# ------------------------------------------------------------------------------------------

bigdf = DataFrame(Bool, 10^3, 10^2)
@time CSV.write("bigdf.csv", bigdf)
@time save("bigdf.jld", "bigdf", bigdf)
getfield.(stat.(["bigdf.csv", "bigdf.jld"]), :size)

# ------------------------------------------------------------------------------------------
# Finally, let's clean up. Do not run the next cell unless you are sure that it will not
# erase your important files.
# ------------------------------------------------------------------------------------------

foreach(rm, ["x.csv", "x.jld", "bigdf.csv", "bigdf.jld"])

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
# This file was generated, do not modify it. # hide
cap_sum_ctry_gd = groupby(capacity, [:country]);
cap_sum_ctry = combine(cap_sum_ctry_gd, :capacity_mw => sum);


text:
using Documenter, Colors

makedocs(
    modules = [Colors],
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    sitename = "Colors",
    pages    = Any[
        "Introduction"             => "index.md",
        "Colorspaces"              => "colorspaces.md",
        "Colorscales"              => "colorscales.md",
        "Colormaps"                => "colormaps.md",
        "Named colors"             => "namedcolors.md",
        "References"               => "references.md",
        "Migrating from Color.jl"  => "migratingfromcolor.md",
        "Index"                    => "functionindex.md",
        ]
    )

deploydocs(
    repo = "github.com/JuliaGraphics/Colors.jl.git",
    target = "build")


text:
abstract type SMExpr end
abstract type Substitutable end

Numeric = Union{Integer, AbstractFloat, Nothing}

struct Num <: SMExpr
    literal::String
    num::Numeric
end

struct Var <: Substitutable
    id::String
end

struct FunExp <: Substitutable
end


text:
using GenInvGaussian
using Distributions
using Base.Test


text:
const Some{AType} = Tuple{AType, Vararg{AType}}

@inline function flatten_unrolled()
    ()
end
@inline function flatten_unrolled(item, rest...)
    item..., flatten_unrolled(rest...)...
end

@inline function reduce_unrolled(call, item)
    item
end
@inline function reduce_unrolled(call, item1, item2, rest...)
    reduce_unrolled(call, call(item1, item2), rest...)
end

@inline function map_unrolled(call, variables::Tuple{})
    ()
end
@inline function map_unrolled(call, variables)
    call(first(variables)), map_unrolled(call, tail(variables))...
end

@inline function map_unrolled(call, variables1::Tuple{}, variables2::Tuple{})
    ()
end
@inline function map_unrolled(call, variables1, variables2)
    call(first(variables1), first(variables2)),
    map_unrolled(call, tail(variables1), tail(variables2))...
end

@inline function partial_map(call, fixed, variables::Tuple{})
    ()
end
@inline function partial_map(call, fixed, variables)
    call(fixed, first(variables)), partial_map(call, fixed, tail(variables))...
end

@inline function partial_map(call, fixed, variables1::Tuple{}, variables2::Tuple{})
    ()
end
@inline function partial_map(call, fixed, variables1, variables2)
    call(fixed, first(variables1), first(variables2)),
    partial_map(call, fixed, tail(variables1), tail(variables2))...
end

@inline function filter_unrolled(call)
    ()
end
@inline function filter_unrolled(call, item, rest...)
    if call(item)
        item, filter_unrolled(call, rest...)...
    else
        rest
    end
end

@inline in_unrolled(needle) = false
@inline in_unrolled(needle, hay, stack...) =
    if needle === hay
        true
    else
        in_unrolled(needle, stack...)
    end
@inline diff_unrolled(less) = ()
@inline diff_unrolled(less, first_more, mores...) =
    if in_unrolled(first_more, less...)
        diff_unrolled(less, mores...)
    else
        first_more, diff_unrolled(less, mores...)...
    end

"""
    over(iterator, call)

Lazy `map` with the reverse argument order.

```jldoctest
julia> using LightQuery

julia> using Test: @inferred

julia> @inferred collect(over([1, -2, -3, 4], abs))
4-element Array{Int64,1}:
 1
 2
 3
 4
```
"""
@inline function over(iterator, call)
    Generator(call, iterator)
end
export over

"""
    when(iterator, call)

Lazy `filter` with the reverse argument order.

```jldoctest
julia> using LightQuery

julia> using Test: @inferred

julia> @inferred collect(when(1:4, iseven))
2-element Array{Int64,1}:
 2
 4
```
"""
@inline function when(iterator, call)
    Iterators.filter(call, iterator)
end
export when

"""
    key(pair)

The `key` in a `key => value` `pair`.

```jldoctest
julia> using LightQuery

julia> using Test: @inferred

julia> @inferred key(:a => 1)
:a

julia> @inferred key((:a, 1))
:a
```
"""
@inline function key((a_key, a_value))
    a_key
end
@inline function key(pair::Pair)
    pair.first
end
export key

"""
    value(pair)

The `value` in a `key => value` `pair`.

```jldoctest
julia> using LightQuery

julia> using Test: @inferred

julia> @inferred value(:a => 1)
1

julia> @inferred value((:a, 1))
1
```
"""
@inline function value((a_key, a_value))
    a_value
end
@inline function value(pair::Pair)
    pair.second
end
export value

"""
    @if_known(something)

If `something` is `missing`, return `missing`, otherwise, `something`.

```jldoctest
julia> using LightQuery

julia> using Test: @inferred

julia> function test(x)
            first(@if_known(x))
        end;

julia> @inferred test((1, 2))
1

julia> @inferred test(missing)
missing
```
"""
macro if_known(something)
    quote
        something = $(esc(something))
        if something === missing
            return missing
        else
            something
        end
    end
end
export @if_known


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct RenewableEnergyConverterTypeA <: Converter
        T_g::Float64
        Rrpwr::Float64
        Brkpt::Float64
        Zerox::Float64
        Lvpl1::Float64
        Vo_lim::Float64
        Lv_pnts::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
        Io_lim::Float64
        T_fltr::Float64
        K_hv::Float64
        Iqr_lims::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
        Accel::Float64
        Lvpl_sw::Int
        Q_ref::Float64
        R_source::Float64
        X_source::Float64
        ext::Dict{String, Any}
        states::Vector{Symbol}
        n_states::Int
    end

Parameters of a renewable energy generator/converter model, this model corresponds to REGCA1 in PSSE

# Arguments
- `T_g::Float64`: Converter time constant (s), validation range: `(0, nothing)`
- `Rrpwr::Float64`: Low Voltage Power Logic (LVPL) ramp rate limit (pu/s), validation range: `(0, nothing)`
- `Brkpt::Float64`: LVPL characteristic voltage 2 (pu), validation range: `(0, nothing)`
- `Zerox::Float64`: LVPL characteristic voltage 1 (pu), validation range: `(0, nothing)`
- `Lvpl1::Float64`: LVPL gain (pu), validation range: `(0, nothing)`
- `Vo_lim::Float64`: Voltage limit for high voltage reactive current management (pu), validation range: `(0, nothing)`
- `Lv_pnts::NamedTuple{(:min, :max), Tuple{Float64, Float64}}`: Voltage points for low voltage active current management (pu) (Lvpnt0, Lvpnt1)
- `Io_lim::Float64`: Current limit (pu) for high voltage reactive current management (specified as a negative value), validation range: `(nothing, 0)`
- `T_fltr::Float64`: Voltage filter time constant for low voltage active current management (s), validation range: `(0, nothing)`
- `K_hv::Float64`: Overvoltage compensation gain used in the high voltage reactive current management, validation range: `(0, nothing)`
- `Iqr_lims::NamedTuple{(:min, :max), Tuple{Float64, Float64}}`: Limit on rate of change for reactive current (pu/s) (Iqr_min, Iqr_max)
- `Accel::Float64`: Acceleration factor, validation range: `(0, 1)`
- `Lvpl_sw::Int`: Low voltage power logic (LVPL) switch. (0: LVPL not present, 1: LVPL present), validation range: `(0, 1)`
- `Q_ref::Float64`: Initial machine reactive power from power flow, validation range: `(0, nothing)`
- `R_source::Float64`: Output resistor used for the Thevenin Equivalent, validation range: `(0, nothing)`
- `X_source::Float64`: Output resistor used for the Thevenin Equivalent, validation range: `(0, nothing)`
- `ext::Dict{String, Any}`
- `states::Vector{Symbol}`: The states are:	Ip: Converter lag for Ipcmd,	Iq: Converter lag for Iqcmd,	Vmeas: Voltage filter for low voltage active current management
- `n_states::Int`: RenewableEnergyConverterTypeA has 3 states
"""
mutable struct RenewableEnergyConverterTypeA <: Converter
    "Converter time constant (s)"
    T_g::Float64
    "Low Voltage Power Logic (LVPL) ramp rate limit (pu/s)"
    Rrpwr::Float64
    "LVPL characteristic voltage 2 (pu)"
    Brkpt::Float64
    "LVPL characteristic voltage 1 (pu)"
    Zerox::Float64
    "LVPL gain (pu)"
    Lvpl1::Float64
    "Voltage limit for high voltage reactive current management (pu)"
    Vo_lim::Float64
    "Voltage points for low voltage active current management (pu) (Lvpnt0, Lvpnt1)"
    Lv_pnts::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
    "Current limit (pu) for high voltage reactive current management (specified as a negative value)"
    Io_lim::Float64
    "Voltage filter time constant for low voltage active current management (s)"
    T_fltr::Float64
    "Overvoltage compensation gain used in the high voltage reactive current management"
    K_hv::Float64
    "Limit on rate of change for reactive current (pu/s) (Iqr_min, Iqr_max)"
    Iqr_lims::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
    "Acceleration factor"
    Accel::Float64
    "Low voltage power logic (LVPL) switch. (0: LVPL not present, 1: LVPL present)"
    Lvpl_sw::Int
    "Initial machine reactive power from power flow"
    Q_ref::Float64
    "Output resistor used for the Thevenin Equivalent"
    R_source::Float64
    "Output resistor used for the Thevenin Equivalent"
    X_source::Float64
    ext::Dict{String, Any}
    "The states are:	Ip: Converter lag for Ipcmd,	Iq: Converter lag for Iqcmd,	Vmeas: Voltage filter for low voltage active current management"
    states::Vector{Symbol}
    "RenewableEnergyConverterTypeA has 3 states"
    n_states::Int
end

function RenewableEnergyConverterTypeA(T_g, Rrpwr, Brkpt, Zerox, Lvpl1, Vo_lim, Lv_pnts, Io_lim, T_fltr, K_hv, Iqr_lims, Accel, Lvpl_sw, Q_ref=1.0, R_source=0.0, X_source=1.0e5, ext=Dict{String, Any}(), )
    RenewableEnergyConverterTypeA(T_g, Rrpwr, Brkpt, Zerox, Lvpl1, Vo_lim, Lv_pnts, Io_lim, T_fltr, K_hv, Iqr_lims, Accel, Lvpl_sw, Q_ref, R_source, X_source, ext, [:Ip, :Iq, :Vmeas], 3, )
end

function RenewableEnergyConverterTypeA(; T_g, Rrpwr, Brkpt, Zerox, Lvpl1, Vo_lim, Lv_pnts, Io_lim, T_fltr, K_hv, Iqr_lims, Accel, Lvpl_sw, Q_ref=1.0, R_source=0.0, X_source=1.0e5, ext=Dict{String, Any}(), states=[:Ip, :Iq, :Vmeas], n_states=3, )
    RenewableEnergyConverterTypeA(T_g, Rrpwr, Brkpt, Zerox, Lvpl1, Vo_lim, Lv_pnts, Io_lim, T_fltr, K_hv, Iqr_lims, Accel, Lvpl_sw, Q_ref, R_source, X_source, ext, states, n_states, )
end

# Constructor for demo purposes; non-functional.
function RenewableEnergyConverterTypeA(::Nothing)
    RenewableEnergyConverterTypeA(;
        T_g=0,
        Rrpwr=0,
        Brkpt=0,
        Zerox=0,
        Lvpl1=0,
        Vo_lim=0,
        Lv_pnts=(min=0.0, max=0.0),
        Io_lim=0,
        T_fltr=0,
        K_hv=0,
        Iqr_lims=(min=0.0, max=0.0),
        Accel=0,
        Lvpl_sw=0,
        Q_ref=0,
        R_source=0,
        X_source=0,
        ext=Dict{String, Any}(),
    )
end

"""Get [`RenewableEnergyConverterTypeA`](@ref) `T_g`."""
get_T_g(value::RenewableEnergyConverterTypeA) = value.T_g
"""Get [`RenewableEnergyConverterTypeA`](@ref) `Rrpwr`."""
get_Rrpwr(value::RenewableEnergyConverterTypeA) = value.Rrpwr
"""Get [`RenewableEnergyConverterTypeA`](@ref) `Brkpt`."""
get_Brkpt(value::RenewableEnergyConverterTypeA) = value.Brkpt
"""Get [`RenewableEnergyConverterTypeA`](@ref) `Zerox`."""
get_Zerox(value::RenewableEnergyConverterTypeA) = value.Zerox
"""Get [`RenewableEnergyConverterTypeA`](@ref) `Lvpl1`."""
get_Lvpl1(value::RenewableEnergyConverterTypeA) = value.Lvpl1
"""Get [`RenewableEnergyConverterTypeA`](@ref) `Vo_lim`."""
get_Vo_lim(value::RenewableEnergyConverterTypeA) = value.Vo_lim
"""Get [`RenewableEnergyConverterTypeA`](@ref) `Lv_pnts`."""
get_Lv_pnts(value::RenewableEnergyConverterTypeA) = value.Lv_pnts
"""Get [`RenewableEnergyConverterTypeA`](@ref) `Io_lim`."""
get_Io_lim(value::RenewableEnergyConverterTypeA) = value.Io_lim
"""Get [`RenewableEnergyConverterTypeA`](@ref) `T_fltr`."""
get_T_fltr(value::RenewableEnergyConverterTypeA) = value.T_fltr
"""Get [`RenewableEnergyConverterTypeA`](@ref) `K_hv`."""
get_K_hv(value::RenewableEnergyConverterTypeA) = value.K_hv
"""Get [`RenewableEnergyConverterTypeA`](@ref) `Iqr_lims`."""
get_Iqr_lims(value::RenewableEnergyConverterTypeA) = value.Iqr_lims
"""Get [`RenewableEnergyConverterTypeA`](@ref) `Accel`."""
get_Accel(value::RenewableEnergyConverterTypeA) = value.Accel
"""Get [`RenewableEnergyConverterTypeA`](@ref) `Lvpl_sw`."""
get_Lvpl_sw(value::RenewableEnergyConverterTypeA) = value.Lvpl_sw
"""Get [`RenewableEnergyConverterTypeA`](@ref) `Q_ref`."""
get_Q_ref(value::RenewableEnergyConverterTypeA) = value.Q_ref
"""Get [`RenewableEnergyConverterTypeA`](@ref) `R_source`."""
get_R_source(value::RenewableEnergyConverterTypeA) = value.R_source
"""Get [`RenewableEnergyConverterTypeA`](@ref) `X_source`."""
get_X_source(value::RenewableEnergyConverterTypeA) = value.X_source
"""Get [`RenewableEnergyConverterTypeA`](@ref) `ext`."""
get_ext(value::RenewableEnergyConverterTypeA) = value.ext
"""Get [`RenewableEnergyConverterTypeA`](@ref) `states`."""
get_states(value::RenewableEnergyConverterTypeA) = value.states
"""Get [`RenewableEnergyConverterTypeA`](@ref) `n_states`."""
get_n_states(value::RenewableEnergyConverterTypeA) = value.n_states

"""Set [`RenewableEnergyConverterTypeA`](@ref) `T_g`."""
set_T_g!(value::RenewableEnergyConverterTypeA, val) = value.T_g = val
"""Set [`RenewableEnergyConverterTypeA`](@ref) `Rrpwr`."""
set_Rrpwr!(value::RenewableEnergyConverterTypeA, val) = value.Rrpwr = val
"""Set [`RenewableEnergyConverterTypeA`](@ref) `Brkpt`."""
set_Brkpt!(value::RenewableEnergyConverterTypeA, val) = value.Brkpt = val
"""Set [`RenewableEnergyConverterTypeA`](@ref) `Zerox`."""
set_Zerox!(value::RenewableEnergyConverterTypeA, val) = value.Zerox = val
"""Set [`RenewableEnergyConverterTypeA`](@ref) `Lvpl1`."""
set_Lvpl1!(value::RenewableEnergyConverterTypeA, val) = value.Lvpl1 = val
"""Set [`RenewableEnergyConverterTypeA`](@ref) `Vo_lim`."""
set_Vo_lim!(value::RenewableEnergyConverterTypeA, val) = value.Vo_lim = val
"""Set [`RenewableEnergyConverterTypeA`](@ref) `Lv_pnts`."""
set_Lv_pnts!(value::RenewableEnergyConverterTypeA, val) = value.Lv_pnts = val
"""Set [`RenewableEnergyConverterTypeA`](@ref) `Io_lim`."""
set_Io_lim!(value::RenewableEnergyConverterTypeA, val) = value.Io_lim = val
"""Set [`RenewableEnergyConverterTypeA`](@ref) `T_fltr`."""
set_T_fltr!(value::RenewableEnergyConverterTypeA, val) = value.T_fltr = val
"""Set [`RenewableEnergyConverterTypeA`](@ref) `K_hv`."""
set_K_hv!(value::RenewableEnergyConverterTypeA, val) = value.K_hv = val
"""Set [`RenewableEnergyConverterTypeA`](@ref) `Iqr_lims`."""
set_Iqr_lims!(value::RenewableEnergyConverterTypeA, val) = value.Iqr_lims = val
"""Set [`RenewableEnergyConverterTypeA`](@ref) `Accel`."""
set_Accel!(value::RenewableEnergyConverterTypeA, val) = value.Accel = val
"""Set [`RenewableEnergyConverterTypeA`](@ref) `Lvpl_sw`."""
set_Lvpl_sw!(value::RenewableEnergyConverterTypeA, val) = value.Lvpl_sw = val
"""Set [`RenewableEnergyConverterTypeA`](@ref) `Q_ref`."""
set_Q_ref!(value::RenewableEnergyConverterTypeA, val) = value.Q_ref = val
"""Set [`RenewableEnergyConverterTypeA`](@ref) `R_source`."""
set_R_source!(value::RenewableEnergyConverterTypeA, val) = value.R_source = val
"""Set [`RenewableEnergyConverterTypeA`](@ref) `X_source`."""
set_X_source!(value::RenewableEnergyConverterTypeA, val) = value.X_source = val
"""Set [`RenewableEnergyConverterTypeA`](@ref) `ext`."""
set_ext!(value::RenewableEnergyConverterTypeA, val) = value.ext = val


text:
"""
Reaction struct.

# Fields
````
id :: String
name :: String
metabolites :: Dict{Metabolite, Float64}
lb :: Float64
ub :: Float64
grr :: Vector{Vector{Gene}}
subsystem :: String
notes :: Dict{String, Vector{String}}
annotation :: Dict{String, Union{Vector{String}, String}}
objective_coefficient :: Float64
````
"""
mutable struct Reaction
    id::String
    name::Maybe{String}
    metabolites::Dict{String,Float64}
    lb::Float64
    ub::Float64
    grr::Maybe{GeneAssociation}
    subsystem::Maybe{String}
    notes::Notes
    annotations::Annotations
    objective_coefficient::Float64

    Reaction(
        id = "";
        name = nothing,
        metabolites = Dict{String,Float64}(),
        lb = -_constants.default_reaction_bound,
        ub = _constants.default_reaction_bound,
        grr = nothing,
        subsystem = nothing,
        notes = Notes(),
        annotations = Annotations(),
        objective_coefficient = 0.0,
    ) = new(
        id,
        name,
        metabolites,
        lb,
        ub,
        grr,
        subsystem,
        notes,
        annotations,
        objective_coefficient,
    )

    # this constructor is only used internally
    function Reaction(
        id::String,
        metabolites::Dict{String,Float64},
        dir = :bidirectional;
        default_bound = _constants.default_reaction_bound,
    )
        if dir == :forward
            lb = 0.0
            ub = default_bound
        elseif dir == :reverse
            lb = -default_bound
            ub = 0.0
        else
            lb = -default_bound
            ub = default_bound
        end
        Reaction(id; metabolites = metabolites, lb = lb, ub = ub)
    end
end


text:
precompile(Tuple{typeof(Base.print), Base.GenericIOBuffer{Array{UInt8, 1}}, UInt16})
precompile(Tuple{typeof(Base.Libc.Libdl.dlopen), String, UInt32})
precompile(Tuple{typeof(Base.invokelatest), Any})
precompile(Tuple{typeof(ColorTypes.__init__)})
precompile(Tuple{typeof(Base.getindex), Base.EnvDict, String})
precompile(Tuple{typeof(Gtk.GLib.GClosureMarshal), Ptr{Nothing}, Ptr{Gtk.GLib.GValue}, UInt32, Ptr{Gtk.GLib.GValue}, Ptr{Nothing}, Ptr{Nothing}})
precompile(Tuple{typeof(Gtk.GLib.uv_prepare), Ptr{Nothing}, Ptr{Int32}})
precompile(Tuple{typeof(Gtk.GLib.uv_check), Ptr{Nothing}})
precompile(Tuple{typeof(Gtk.GLib.uv_dispatch), Ptr{Nothing}, Ptr{Nothing}, Int64})
precompile(Tuple{typeof(Gtk.GLib.g_yield), UInt64})
precompile(Tuple{typeof(Gtk.GLib.__init__)})
precompile(Tuple{typeof(Base.cconvert), Type, Int64})
precompile(Tuple{typeof(Gtk.__init__)})
precompile(Tuple{typeof(Base.filter), Gtk.var"#311#317", Array{String, 1}})
precompile(Tuple{Type{NamedTuple{(:validate_strict, :libc, :call_abi, :libgfortran_version, :cxxstring_abi, :libstdcxx_version, :os_version), T} where T<:Tuple}, Tuple{Bool, String, Nothing, Base.VersionNumber, String, Base.VersionNumber, Nothing}})
precompile(Tuple{typeof(Base.merge), NamedTuple{(:validate_strict, :libc, :call_abi, :libgfortran_version, :cxxstring_abi, :libstdcxx_version, :os_version), Tuple{Bool, String, Nothing, Base.VersionNumber, String, Base.VersionNumber, Nothing}}, Array{Base.Pair{Symbol, Base.SubString{String}}, 1}})
precompile(Tuple{Base.var"##s79#169", Any, Any, Any, Any, Any})
precompile(Tuple{typeof(Base.merge_types), Tuple{Vararg{Symbol, N} where N}, Type{var"#s79"} where var"#s79"<:(NamedTuple{names, T} where T<:Tuple where names), Type{var"#s78"} where var"#s78"<:(NamedTuple{names, T} where T<:Tuple where names)})
precompile(Tuple{typeof(Base.merge), NamedTuple{(:validate_strict, :libc, :call_abi, :libgfortran_version, :cxxstring_abi, :libstdcxx_version, :os_version), Tuple{Bool, String, Nothing, Base.VersionNumber, String, Base.VersionNumber, Nothing}}, NamedTuple{(:julia_version,), Tuple{Base.SubString{String}}}})
precompile(Tuple{Core.var"#Type##kw", NamedTuple{(:validate_strict, :libc, :call_abi, :libgfortran_version, :cxxstring_abi, :libstdcxx_version, :os_version, :julia_version), Tuple{Bool, String, Nothing, Base.VersionNumber, String, Base.VersionNumber, Nothing, Base.SubString{String}}}, Type{Base.BinaryPlatforms.Platform}, String, String})
precompile(Tuple{Artifacts.var"##artifact_hash#13", Base.BinaryPlatforms.Platform, Nothing, typeof(Artifacts.artifact_hash), String, String})
precompile(Tuple{Gtk.GLib.var"#100#101"})

max_stars_repo_name: Lecrapouille/BacASable

max_issues_repo_name: Lecrapouille/BacASable

max_forks_repo_name: Lecrapouille/BacASable


text:
"""
Utilities to add cost functions together
"""

"""
A struct for holding differentiable cost functions. Let d::DiffCost.
d(p) returns cost at p
d(p,g) returns cost and **mutates** gradient at p
"""
struct DiffCost{F,F2} <: Function
    cost_function::F
    cost_function2::F2
end
(f::DiffCost)(x) = f.cost_function(x)
(f::DiffCost)(x,y) = f.cost_function2(x, y)

"""
given a cost function C(p), makes a new cost d::DiffCost, which has a method for returning the finite-difference gradient.
"""
function make_fd_differentiable(cost)

    function cost2(p, g)
        FiniteDiff.finite_difference_gradient!(g, cost, p)
        return cost(p)
    end
    return DiffCost(cost, cost2)
end

"""
    sum_losses(lArray::Array{T,1}, p0) where T <: Function
- Given array of cost functions c1...cn, makes new cost function D. Multithreads evaluation of D(p) and D(p,g) if threads are available.
- c1...cn must be differentiable: ci(p,g) returns cost and mutates gradient g

D(p) = sum_i c_i(p)
D(p,g) mutates g to give gradient of D.

"""
function sum_losses(lArray::Array{T,1}, p0) where T <: Function
    (Threads.nthreads() == 1) && (@info "Note that restarting julia with multiple threads will increase performance of the generated loss function from sum_losses()")
    dummy_gs = [deepcopy(p0) for el in lArray]::Vector{typeof(p0)}
    cs = [0. for el in lArray]
    function cost1(p)
        return sum(lArray) do loss
            loss(p)
        end
    end
    n = length(lArray)
    cs = Vector{Float64}(undef, n)
    pure_costs = map(lArray) do el
        (p, g) -> (el[i](p, g), g)
    end

    function cost2(p, g)
        ThreadsX.foreach(enumerate(lArray)) do (i, el)
            cs[i] = el(p, dummy_gs[i])
        end
        g[:] = sum(dummy_gs)
        return sum(cs)
    end
    return DiffCost(cost1, cost2)
end


text:
create(filename, A::AbstractRaster{T}; kw...) where T = create(filename, T, A; kw...)
function create(filename, T, A::AbstractRaster;
    name=name(A), metadata=metadata(A), missingval=missingval(A), kw...
)
    create(filename, T, dims(A); parent=parent(A), name, metadata, missingval, kw...)
end
function create(filename::AbstractString, T::Type, dims::Tuple; parent=nothing, suffix=nothing, kw...)
    filename = _maybe_add_suffix(filename, suffix)
    create(filename, _sourcetype(filename), T, dims; kw...)
end
function create(filename::Nothing, T::Type, dims::Tuple; parent=nothing, suffix=nothing, missingval, kw...)
    T = isnothing(missingval) ? T : promote_type(T, typeof(missingval))
    data = isnothing(parent) ? Array{T}(undef, dims) : similar(parent, T, size(dims))
    Raster(data, dims; missingval, kw...)
end

_maybe_add_suffix(filename::Nothing, suffix) = nothing
_maybe_add_suffix(filename::Nothing, suffix::Nothing) = nothing
_maybe_add_suffix(filename, suffix::Nothing) = filename
function _maybe_add_suffix(filename, suffix)
    base, ext = splitext(filename)
    if string(suffix) == ""
        filename
    else
        return string(base, "_", suffix, ext)
    end
end


text:
using Test
using SafeTestsets

include("generate_example_tests.jl")

# Note: comment outer @testset to stop after first @safetestset failure
@time @testset verbose=true "Krotov.jl Package" begin

    print("\n* Example 1 (examples/simple_state_to_state.jl):")
    @time @safetestset "Example 1" begin include(joinpath("examples", "simple_state_to_state.jl")) end

    print("\n* Example 2 (examples/state_to_state_rwa.jl):")
    @time @safetestset "Example 2" begin include(joinpath("examples", "state_to_state_rwa.jl")) end

    print("\n* Example 3 (examples/rho_3states.jl):")
    if Sys.isapple()
        println("\nSkipped (macOS)")
    else
        @time @safetestset "Example 3" begin include(joinpath("examples", "rho_3states.jl")) end
    end

    print("\n* Example 4 (examples/state_to_state_parametrizations.jl):")
    @time @safetestset "Example 4" begin include(joinpath("examples", "state_to_state_parametrizations.jl")) end

    println("")

end


text:
# This example demonstrated how to calculate one-loop diagram of free electrons using the Monte Carlo module
# Observable is normalized: Γ₄*N_F where N_F is the free electron density of states

using Distributed
using QuantumStatistics, LinearAlgebra, Random, Printf, StaticArrays, BenchmarkTools, InteractiveUtils, Parameters

const Ncpu = 4 # number of workers (CPU)
const totalStep = 1e8 # MC steps of each worker

addprocs(Ncpu)

@everywhere using QuantumStatistics, Parameters, StaticArrays, Random, LinearAlgebra
@everywhere include("../application/electron_gas/RPA.jl")

# parameters
@everywhere @with_kw struct Para
    kF::Float64 = (9π/4.0)^(1.0/3)
    m::Float64 = 0.5
    β::Float64 = 1000.0 / kF^2 *2m
    spin::Int = 2
    n::Int = 0 # external Matsubara frequency
end

################ construct RPA interaction ################################
@everywhere const para = Para()
@everywhere const e2tors = 2.0/(9π/4.0)^(1.0/3)
@everywhere const kF, m, e, spin, AngSize = para.kF, para.m, sqrt(1.0*e2tors), para.spin, 32
@everywhere const mass2 = 0.0#0.000000001
@everywhere const β, EF = para.β, kF^2 / (2m)

@everywhere const qgrid = Grid.boseKUL(kF, 6kF, 0.000001*sqrt(m^2/β/kF^2), 15,4)
@everywhere const τgrid = Grid.tauUL(β, 0.0001, 11,4)
@everywhere const vqinv = [(q^2 + mass2) / (4π * e^2) for q in qgrid.grid]
println(qgrid.grid)
println(τgrid.grid)

@everywhere const dW0 = dWRPA(vqinv, qgrid.grid, τgrid.grid, kF, β, spin, m) # dynamic part of the effective interaction

@everywhere function interaction(q, τIn, τOut)
    dτ = abs(τOut - τIn)

    kQ = sqrt(dot(q, q))
    v = 4π * e^2 / (kQ^2 + mass2)
    if kQ <= qgrid.grid[1]
        w = v * Grid.linear2D(dW0, qgrid, τgrid, qgrid.grid[1] + 1.0e-14, dτ) # the current interpolation vanishes at q=0, which needs to be corrected!
    else
        w = v * Grid.linear2D(dW0, qgrid, τgrid, kQ, dτ) # dynamic interaction, don't forget the singular factor vq
    end

    return -v / β, -w
end

@everywhere function integrand(config)
    if config.curr == 1
        return 1.0 # return a weight!
    elseif config.curr == 2
        return eval2(config)
    else
        error("Not implemented!")
    end
end

@everywhere function eval2(config)
    para = config.para
    β, kF, m, spin = para.β, para.kF, para.m, para.spin

    T, K = config.var[1], config.var[2]
    k1= K[1]
    t1, t2 = T[1], T[2]
    KIN = @SVector[kF,0,0]

    v,w = interaction(k1, t1, t2)

    kq=KIN+k1
    τ = (t2-t1) / β

    ϵ1 = (dot(kq, kq) - kF^2) / (2m)
    # possible green's functions on the top
    g1 = Spectral.kernelFermiT(τ, ϵ1)
    g0 = Spectral.kernelFermiT(t1 - t1, ϵ1, β)

    phase = 1.0 / (2π)^3
    return spin * g1 * w * phase * (-τ * β) * sin(π * (2*para.n+1) * abs(τ)) #+ g0 * v * phase
end

@everywhere function measure(config)
    obs = config.observable
    curr = config.curr
    factor = 1.0 / config.reweight[curr]
    if curr == 1
        obs[1][1] += factor
    elseif curr == 2
        weight = integrand(config)
        obs[2][1] += weight / abs(weight) * factor
    else
        return
    end
end

@everywhere normalize(config) = config.observable[2] / sum(config.observable[1])

function run(totalStep)

    para = Para()
    @unpack kF, β = para

    K = MonteCarlo.FermiK(3, kF, 0.2 * kF, 10.0 * kF)
    T = MonteCarlo.Tau(β, β / 2.0)

    dof = ([1, 0], [2, 1]) # degrees of freedom of the normalization diagram and the bubble

    obs = (zeros(Float64, 1), zeros(Float64, 1)) # observable for the normalization diagram and the bubble

    avg, std = MonteCarlo.sample(totalStep, (T, K), dof, obs, integrand, measure, normalize; para=para, print=10)

    @unpack kF, β, m, n = Para()

    @printf("%10.6f ± %10.6f\n", avg[1], std[1])
    @printf("%10.6f ± %10.6f\n", 1.0/(1.0 - avg[1]),std[1]/(1.0-avg[1])^2)
end

run(totalStep)
# @time run(totalStep)


text:
#=

Create ODE numerical integration function

Ignacio Quintero Mächler

t(-_-t)

September 19 2017

=#

"""
    make_solver(odef, p0::Array{Float64,1})
Make **ODE** solver for `odef`.
"""
function make_solver(odef, p0::Array{Float64,1}, u0::Array{Float64,1})

  prob = ODEProblem(odef, u0, (0.0,1.0), p0)

  int = init(prob,
             Tsit5(),
             save_everystep  = false,
             calck           = false,
             force_dtmin     = true,
             save_start      = false,
             initialize_save = false,
             maxiters        = 100_000_000,
             verbose         = false)

  function f(u ::Array{Float64,1},
             p ::Array{Float64,1},
             ti::Float64,
             tf::Float64)
    int.p = p
    reinit!(int, u, t0 = ti, tf = tf)
    solve!(int).u[1]
  end
  return f
end

"""
    solvef(int::DEIntegrator,
           u  ::Array{Float64,1},
           ti ::Float64,
           tf ::Float64)

Solve an IDE integrator for new `u` and and `ti` and `tf`.
"""
function solvef(int::DiffEqBase.DEIntegrator,
                u  ::Array{Float64,1},
                ti ::Float64,
                tf ::Float64)
  @inbounds begin
    reinit!(int, u, t0 = ti, tf = tf)
    return solve!(int).u[1]
  end
end

max_stars_repo_name: hmorlon/JPANDA

max_issues_repo_name: hmorlon/JPANDA

max_forks_repo_name: hmorlon/JPANDA


text:
using DataFrames
using Test

m = test_page_model()
include("../src/components/AdaptationCosts.jl")

adaptationcosts_noneconomic = MimiPAGE2009.addadaptationcosts_noneconomic(m)
adaptationcosts_economic = MimiPAGE2009.addadaptationcosts_economic(m)
adaptationcosts_sealevel = MimiPAGE2009.addadaptationcosts_sealevel(m)

set_param!(m, :y_year_0, 2008)
set_param!(m, :y_year, Mimi.dim_keys(m.md, :time))
set_param!(m, :gdp, readpagedata(m, "test/validationdata/gdp.csv"))
set_param!(m, :automult_autonomouschange, 0.65)

p = load_parameters(m)
set_leftover_params!(m, p)

run(m)

autofac = m[:AdaptiveCostsEconomic, :autofac_autonomouschangefraction]
atl_economic = m[:AdaptiveCostsEconomic, :atl_adjustedtolerablelevel]
imp_economic = m[:AdaptiveCostsEconomic, :imp_adaptedimpacts]
acp_economic = m[:AdaptiveCostsEconomic, :acp_adaptivecostplateau]
aci_economic = m[:AdaptiveCostsEconomic, :aci_adaptivecostimpact]

autofac_compare = readpagedata(m, "test/validationdata/autofac_autonomouschangefraction.csv")
atl_economic_compare = readpagedata(m, "test/validationdata/atl_1.csv")
imp_economic_compare = readpagedata(m, "test/validationdata/imp_1.csv")
acp_economic_compare = readpagedata(m, "test/validationdata/acp_adaptivecostplateau_economic.csv")
aci_economic_compare = readpagedata(m, "test/validationdata/aci_adaptivecostimpact_economic.csv")

@test autofac ≈ autofac_compare rtol=1e-8
@test atl_economic ≈ atl_economic_compare rtol=1e-8
@test imp_economic ≈ imp_economic_compare rtol=1e-8
@test acp_economic ≈ acp_economic_compare rtol=1e-3
@test aci_economic ≈ aci_economic_compare rtol=1e-3

ac_noneconomic = m[:AdaptiveCostsNonEconomic, :ac_adaptivecosts]
ac_economic = m[:AdaptiveCostsEconomic, :ac_adaptivecosts]
ac_sealevel = m[:AdaptiveCostsSeaLevel, :ac_adaptivecosts]

ac_noneconomic_compare = readpagedata(m, "test/validationdata/ac_adaptationcosts_noneconomic.csv")
ac_economic_compare = readpagedata(m, "test/validationdata/ac_adaptationcosts_economic.csv")
ac_sealevel_compare = readpagedata(m, "test/validationdata/ac_adaptationcosts_sealevelrise.csv")

@test ac_noneconomic ≈ ac_noneconomic_compare rtol=1e-2
@test ac_economic ≈ ac_economic_compare rtol=1e-3
@test ac_sealevel ≈ ac_sealevel_compare rtol=1e-2


text:
using StatsBase, LinearAlgebra
include("neuralNet.jl")
include("propagation.jl")
include("pulsesGandD.jl")

Eα_r = LinRange(0.0,1.25,1200)
δE = Eα_r[2] - Eα_r[1]
Eαr = h5read("Energy_V_alpha_beta.h5", "Energy_Levels");

data2 = h5open("dataSynthst10kfv2.h5", "r")
coeffsNoi = read(data2["coeffsNoi"])
xdata = reshape(coeffsNoi, (60,10*10000));

coeffsRefn = read(data2["coeffsRef"]);
bestTs_ref = read(data2["bestTs_ref"]);
ydata = coeffsRefn;

#duplicate to pair up with all references double pulse
function get_ydata(ref_ydata,i,j; nc= 40,sampref = 10)
    # i,j = 1,1, Guassian pulses
    outdata = zeros(Float64, nc, size(ref_ydata)[4]*sampref)
    ijdata = ref_ydata[:, i,j,:]
    increment = 1
    for i in 1:size(ydata)[4]
        outdata[:, increment:increment+sampref-1] = repeat(ijdata[:,i],1,sampref)
        increment += sampref
    end
    outdata
end
function splitdata(xdata, ydata)
    #μnoi = mean(xdata)
    #stdnoi = std(xdata)
    #μref = mean(ydata_1)
    #stdref = std(ydata_1)
    #xdata = (xdata .- μnoi)./stdnoi
    #ydata = (ydata .- μref)./stdref

    xtrain = xdata[:,1:70000]
    ytrain = ydata[:,1:70000]
    xval = xdata[:,70001:80000]
    yval = ydata[:,70001:80000]
    xtest = xdata[:,80001:end]
    ytest = ydata[:,80001:end]
    xtrain, ytrain, xval, yval, xtest, ytest #, [μnoi, stdnoi, μref, stdref]
end

i = 3
j = 3
ydataDij = get_ydata(ydata,i,j)
xtrain, ytrain, xval, yval, xtest, ytest = splitdata(xdata, ydataDij)
dict_models, losses = ensembleMachines(xtrain, ytrain, xval, yval;
                                                    nmodels = 10, nepoch = 100)

max_stars_repo_name: lazarusA/LightMatterInteraction

max_issues_repo_name: lazarusA/LightMatterInteraction

max_forks_repo_name: lazarusA/LightMatterInteraction


text:
using AWS.EC2
using AWS

const ami_ubuntu_13_04_64_bit_US_east = "ami-64a4de0d" # Julia nightly 25th july
const ami_ubuntu_13_04_32_bit_US_east = "ami-5b8be132" # With Julia

const ami_ubuntu_13_04_32bit_USeast_dev = "ami-8bdda8e2"

env = AWSEnv()
# env.dbg=true
# env.dry_run=true

CHK_ERR(resp::EC2Response) = (typeof(resp.obj) == EC2Error) ? error(ec2_error_str(resp.obj)) : resp.obj

function desc_avail_zones()
    resp = CHK_ERR(DescribeAvailabilityZones(env))
    resp
end

function terminate_instances (instances)
    req = TerminateInstancesType(instancesSet=instances)
    resp = CHK_ERR(TerminateInstances(env, req))
    instances
end

function terminate_instances_by_owner (owner::String)
    instances = get_instances_by_owner (owner)
    println("$owner has the following [$instances] instances")
    terminate_instances (instances)
end

function get_running_instances_by_owner (owner::String)
    instances = ASCIIString[]
    req = DescribeInstancesType(filterSet=[FilterType(name="tag:Owner", valueSet=[owner])])
    resp = CHK_ERR(DescribeInstances(env, req))
    reservs = resp.reservationSet
    for reserv in reservs
        instancesSet = reserv.instancesSet
        for i in instancesSet
            if i.instanceState.code == 16
                push!(instances, i.instanceId)
            end
        end
    end
    instances
end

function check_running(chk_instances)
    new_chk_instances = ASCIIString[]
    req = DescribeInstanceStatusType(instancesSet=chk_instances, includeAllInstances=true)
    resp = CHK_ERR(DescribeInstanceStatus(env, req))
    statuses = resp.instanceStatusSet
    for status in statuses
        println("Status of $(status.instanceId) is $(status.instanceState.code):$(status.instanceState.name)")
        if status.instanceState.code != 16
            push!(new_chk_instances, status.instanceId)
        end
    end
    new_chk_instances
end

function get_hostnames(instances)
    names = NTuple[]
    req = DescribeInstancesType(instancesSet=instances)
    resp = CHK_ERR(DescribeInstances(env, req))
    reservs = resp.reservationSet
    for reserv in reservs
        instancesSet = reserv.instancesSet
        for i in instancesSet
            push!(names, (i.instanceId, i.dnsName))
        end
    end
    names
end

function launch_n_ec2(n::Int)
# "c1.medium"
# "m1.small"

    resp = CHK_ERR(RunInstances(env, RunInstancesType(imageId=ami_ubuntu_13_04_32bit_USeast_dev, instanceType="m1.small", minCount=n, maxCount=n, keyName="jublr")))
    instances = ASCIIString[]
    for inst in resp.instancesSet
        push!(instances, inst.instanceId)
    end
    println("Launched instances : $instances" )

    # Tag the instances..also test ec2_basic API..
    resp = ec2_basic(env, "CreateTags", {"ResourceId"=>instances, "Tag"=>[{"Key"=>"Name","Value"=>"AWSTest"}, {"Key"=>"Owner","Value"=>"amitm"}]})
    if (typeof(resp.obj) == EC2Error)
        error(ec2_error_str(resp.obj))
    else
        resp.obj = CreateTagsResponseType(resp.pd)
        if resp.obj._return == true
            println("Successfully added tags!")
        else
            error("error adding tags!")
        end
    end
    println("Tagged instances : $instances" )

    # Wait for the instances to come to a running state....
    start = time()
    chk_instances = check_running(instances)
    while (((time() - start) < 60.0) && (length(chk_instances) > 0))
        println("All instances not yet in a running state. Waiting and trying again....")
        sleep(5.0)

        chk_instances = check_running(chk_instances)
    end

    if (length(chk_instances) > 0)
        println("All instances not yet in a running state. Please check after some time.")
    end
    instances
end


text:
using MPI
using DoubleFloats

struct MPIData
    comm::Any
    nb::Any
    ind::Any
    function MPIData()
        if !MPI.Initialized()
            MPI.Init()
        end
        comm = MPI.COMM_WORLD
        nb = MPI.Comm_size(comm)
        ind = MPI.Comm_rank(comm) + 1
        return new(comm, nb, ind)
    end
end

function mpibroadcast(mpid, t_split, data::Array{T,N}) where {T,N}

    MPI.Barrier(mpid.comm)
        if T == Float64
            # for Float64 or Double64 ... for example
            for i = 1:mpid.nb
                vbcast = view(data, t_split[i])
                MPI.Bcast!(vbcast, i - 1, mpid.comm)
            end
        else
            # for BigFloat ... for example
            for i = 1:mpid.nb
                vbcast = view(data, t_split[i])
                bufr = MPI.bcast(vbcast, i - 1, mpid.comm)
                if i != mpid.ind
                    copy!(vbcast, bufr)
                end
            end
        end

    MPI.Barrier(mpid.comm)
end
function splititr(nb, lgtot)
    lg, r = divrem(lgtot, nb)
    return vcat(
        map(x -> ((x-1)*(lg+1)+1):x*(lg+1), 1:r),
        map(x -> ((x-1)*lg+r+1):(x*lg+r), (r+1):nb),
    )
end

function splitvec(nb, v)
    return map(x -> v[x], splititr(nb, length(v)))
end

function fct(T, sz::NTuple{2,Int})
    mpid = MPIData()
    nbsplit = mpid.nb

    t_itr = (splitvec(nbsplit, CartesianIndices(sz)),)

    @show typeof(t_itr)

    li = LinearIndices(sz)
    it = t_itr[1]
    tt_split = ntuple(y -> (li[it[y][1]]:li[it[y][end]]), nbsplit)

    @show typeof(tt_split)

    data = zeros(T,sz)
    cpt = 1
    loc_t = ref_t = time_ns()
    while true
        dec1 = rand()-0.5
        dec2 = rand()-0.5
        for i in tt_split[mpid.ind]
            a = div(i,sz[1])
            b = mod(i,sz[2])
            data[i] = (data[i] + dec1*a + dec2*b)/(a+b)
        end
        mpibroadcast(mpid, tt_split, data)
        if mpid.ind == 1
                if cpt % 80 == 0
                    n_t = time_ns()
                    loc = (n_t - loc_t)*1e-9
                    gen = (n_t - ref_t)*1e-9
                    moy = gen *80 / cpt
                    loc_t = n_t
                    println(". $(sum(data)) $cpt $moy $loc $gen")
                else
                    if cpt % 1 == 0
                        print(".")
                    end
                end
        end
        flush(stdout)
        cpt += 1
    end

end
function run_main(T,sz)
    fct(T,sz)
end

run_main(Float64, (512,512) )


text:
## ExoplanetsSysSim/src/tess_catalog.jl
## (c) 2020 Aditya Sengupta?

using DataFrames
using CSV
using FileIO

mutable struct TESSPhysicalCatalog
  target::Array{TESSTarget,1}
end

mutable struct TESSObsCatalog
  target::Array{TESSTargetObs,1}
end
TESSObsCatalog() = TESSObsCatalog(TESSTargetObs[])

"""
    generate_tess_physical_catalog(sim_param)

Wrapper function to create catalog of simulated TESS targets.

# Arguments:
- `sim_param::SimParam`: Simulation parameter object; must have the following parameters set:
  - num_targets_sim_pass_one = Number of TESS targets in simulated catalog
  - generate_tess_target = Function which generates TESS targets
  - (Optional) stellar_catalog = Stellar catalog filename
  - (Optional) star_table_setup = Function that loads stellar catalog into DataFrame

# Returns:
TESS physical catalog object containing all simulated TESS targets and associated planetary systems.
"""
function generate_tess_physical_catalog(sim_param::SimParam)
   if haskey(sim_param,"stellar_catalog")
      star_tab_func = get_function(sim_param, "star_table_setup")
      star_tab_func(sim_param)
   end
   num_sys = get_int(sim_param,"num_targets_sim_pass_one")
   generate_tess_target = get_function(sim_param,"generate_tess_target")
   target_list = Array{TESSTarget}(undef,num_sys)
   map!(x->generate_tess_target(sim_param), target_list, 1:num_sys )
   return TESSPhysicalCatalog(target_list)
end

"""
    observe_tess_targets_sky_avg(input, sim_param)

Wrapper function to create catalog of simulated observations of TESS targets using the sky averaging observation scheme (i.e. each planet's detection probability is the average detection probability over all view-angles).

# Arguments:
- `input::TESSPhysicalCatalog`: Catalog object of simulated TESS targets and associated planetary systems to be observed
- `sim_param::SimParam`: Simulation parameter object; requires the following simulation parameters to be set:
  - calc_target_obs_sky_ave: Function name for sky averaging simulated observations

# Returns:
TESS observations catalog object containing all properties observed from the TESS targets and associated planetary systems that were detected during the simulation.
"""
function observe_tess_targets_sky_avg(input::TESSPhysicalCatalog, sim_param::SimParam )
  calc_target_obs = get_function(sim_param,"calc_target_obs_sky_ave")
  return observe_tess_targets(calc_target_obs, input, sim_param)
end

"""
    observe_tess_targets_single_obs(input, sim_param)

Wrapper function to create catalog of simulated observations of TESS targets using the single observer observation scheme (i.e. each planet's detection probability is the detection probability from the Earth).

# Arguments:
- `input::TESSPhysicalCatalog`: Catalog object of simulated TESS targets and associated planetary systems to be observed
- `sim_param::SimParam`: Simulation parameter object; requires the following simulation parameters to be set:
  - calc_target_obs_single_obs: Function name for single observer simulated observations

# Returns:
TESS observations catalog object containing all properties observed from the TESS targets and associated planetary systems that were detected during the simulation.
"""
function observe_tess_targets_single_obs(input::TESSPhysicalCatalog, sim_param::SimParam )
  calc_target_obs = get_function(sim_param,"calc_target_obs_single_obs")
  return observe_tess_targets(calc_target_obs, input, sim_param)
end

"""
    observe_tess_targets(calc_target_obs, input, sim_param)

Wrapper function to create catalog of simulated observations of TESS targets.

# Arguments:
- `calc_target_obs::Function`: Function to use in simulating observations of TESS targets (sky averaging vs. single observer schemes)
- `input::TESSPhysicalCatalog`: Catalog object of simulated TESS targets and associated planetary systems to be observed
- `sim_param::SimParam`: Simulation parameter object

# Returns:
TESS observations catalog object containing all properties observed from the TESS targets and associated planetary systems that were detected during the simulation.
"""
function observe_tess_targets(calc_target_obs::Function, input::TESSPhysicalCatalog, sim_param::SimParam )
  #calc_target_obs = get_function(sim_param,"calc_target_obs_sky_ave")
  #calc_target_obs = get_function(sim_param,"calc_target_obs_single_obs")
  output = TESSObsCatalog()
  if haskey(sim_param,"mem_tess_target_obs")
     output.target = get(sim_param,"mem_tess_target_obs",Array{TESSTargetObs}(0) )
  end
  num_targets_sim_pass_one = get_int(sim_param,"num_targets_sim_pass_one")
  if length(output.target) < num_targets_sim_pass_one
     output.target = Array{TESSTargetObs}(undef,num_targets_sim_pass_one)
  end
  map!(x::TESSTarget->calc_target_obs(x,sim_param)::TESSTargetObs, output.target, input.target)
  resize!(output.target,length(input.target))
  return output
end

# Remove undetected planets from physical catalog
# TODO: OPT: Maybe create array of bools for which planets to keep, rather than splicing out non-detections?
function generate_obs_targets(cat_phys::TESSPhysicalCatalog, sim_param::SimParam )
  for t in 1:length(cat_phys.target)
    kep_targ = cat_phys.target[t]
    for ps in 1:length(cat_phys.target[t].sys)
      sys = kep_targ.sys[ps]
      for pl in length(sys.orbit):-1:1    # Going in reverse since removing planets from end of list first is cheaper than starting at beginning
        ecc::Float64 = sys.orbit[pl].ecc
	incl::Float64 = sys.orbit[pl].incl
   	a::Float64 = semimajor_axis(sys,pl)
   	Rstar::Float64 = rsol_in_au*sys.star.radius

        does_it_transit = does_planet_transit(sys, pl)
        pdet_if_tr = does_it_transit ? calc_prob_detect_if_transit_with_actual_b(kep_targ, ps, pl, sim_param) : 0.
        if !does_it_transit || (rand()>pdet_if_tr)
    	  splice!(cat_phys.target[t].sys[ps].orbit, pl)
	  splice!(cat_phys.target[t].sys[ps].planet, pl)
     	end
      end
    end
  end
  return cat_phys
end

# The following function is primarily left for debugging.
# Create a catalog of observations of simulated TESS targets.
function simulated_read_tess_observations(sim_param::SimParam )
   println("# WARNING: Using simulated_read_tess_observations.")
   # if haskey(sim_param,"stellar_catalog")
   #    star_tab_func = get_function(sim_param, "star_table_setup")
   #    star_tab_func(sim_param)
   # end
   num_sys = get_int(sim_param,"num_tess_targets")
   generate_tess_target = get_function(sim_param,"generate_tess_target")
   target_list = Array{TESSTarget}(undef,num_sys)
   map!(x->generate_tess_target(sim_param), target_list, 1:num_sys )

   cat_phys_cut = generate_obs_targets(TESSPhysicalCatalog(target_list), sim_param)
   calc_target_obs = get_function(sim_param,"calc_target_obs_single_obs")
   output = TESSObsCatalog()
   output.target = map(x::TESSTarget->calc_target_obs(x,sim_param)::TESSTargetObs, cat_phys_cut.target)
   return output
end

"""
    read_toi_catalog(sim_param, force_reread=false)

Wrapper function to read TESS Object of Interest (TOI) catalog given SimParam

# Arguments:
- `sim_param::SimParam`: Simulation parameter object; this function uses the following parameters from the SimParam object:
  - toi_catalog: String filename of TESS Object of Interest catalog
- `force_reread::Bool`: Should the file be read in even if a DataFrame of the TOIs already exists in workspace?

# Returns:
- Dataframe of TOI objects (and their respective properties).
- Vector of booleans indicating which TOIs  were designated as planet candidates by the TESS pipeline and have a valid observed radius ratio and period (necessary for detection probability calculation).
"""
function read_toi_catalog(sim_param::SimParam, force_reread::Bool = false)
    filename = convert(String,joinpath(dirname(pathof(ExoplanetsSysSim)),"..", "data", convert(String,get(sim_param,"toi_catalog", "tesstargets/toi_catalog.csv")) ) )
    return read_toi_catalog(filename, force_reread)
end

"""
    read_toi_catalog(filename, force_reread=false)

Function to read TESS Object of Interest (TOI) catalog given filename string.

# Arguments:
- `filename::String`: String filename of TESS Object of Interest catalog
- `force_reread::Bool`: Should the file be read in even if a DataFrame of the TOIs already exists in workspace?

# Returns:
- Dataframe of TOI objects (and their respective properties).
- Vector of booleans indicating which TOIs were designated as planet candidates by the TESS pipeline and have a valid observed radius ratio and period (necessary for detection probability calculation).
"""
function read_toi_catalog(filename::String, force_reread::Bool = false)
    local df, usable

    if occursin(r".jld2$",filename) && !force_reread
        try
            data = load(filename)
            df = data["toi_catalog"]
            usable = data["toi_catalog_usable"]
            Core.typeassert(df,DataFrame)
            Core.typeassert(usable,Array{Int64,1})
        catch
            error(string("# Failed to read toi catalog >",filename,"< in jld2 format."))
        end
    else
       try
            df = CSV.read(filename,comment="#")

            # Choose which TOIs to keep
            #is_cand = (csv_data[!,:,toi_disposition_idx] .== "CONFIRMED") | (csv_data[!,:,toi_disposition_idx] .== "CANDIDATE")
            is_cand = df[!,:toi_pdisposition] .== "PC"
             # TODO: figure out a way to replace missing errors without zeroing them out (add in a baseline Price/Rogers variance?)
            has_radius = .!ismissing.(df[!,:toi_prad])
            has_period = .!(ismissing.(df[!,:toi_period]) .| ismissing.(df[!,:toi_period_err]))
            has_epoch = .!(ismissing.(df[!,:epoch]) .| ismissing.(df[!,:epoch_err]))
            has_depth = .!(ismissing.(df[!,:toi_transit_depth]) .| ismissing.(df[!,:toi_transit_depth_err]))
            is_usable = .&(is_cand, has_radius, has_period, has_depth, has_epoch)
            usable = findall(is_usable)
           #  symbols_to_keep = [:TIC, :kepoi_name, :toi_pdisposition, :toi_score, :toi_ror, :toi_period, :toi_period_err1, :toi_period_err2, :toi_time0bk, :toi_time0bk_err1, :toi_time0bk_err2, :toi_depth, :toi_depth_err1, :toi_depth_err2, :toi_duration, :toi_duration_err1, :toi_duration_err2]
           # df = df[usable, symbols_to_keep]
           # tmp_df = DataFrame()
           # for col in names(df)
           #     tmp_df[col] = collect(skipmissing(df[col]))
           # end
           # df = tmp_df
           # usable = collect(1:length(df[!,:TIC]))
        catch
            error(string("# Failed to read toi catalog >",filename,"< in ascii format."))
        end
    end
    df[!,:toi_prad] *= (earth_radius_eq_in_m_IAU2015 / sun_radius_in_m_IAU2015)
    return df, usable
end

"""
    setup_actual_planet_candidate_catalog(df_star, df_toi, usable_toi, sim_param)

Create (true) catalog of TESS observations of TESS targets

# Arguments:
- `df_star::DataFrame`: DataFrame containing all TESS target stars in catalog
   NOTE: df_star is assumed to have fields TIC, mass and radius for all targets in the survey)
- `df_toi::DataFrame`: DataFrame containing all TESS Object of Interests (TOIs)
- `usable_toi::Array{Integer}`: Array of TOI dataframe row indices corresponding to TOIs to use
- `sim_param::SimParam`: Simulation parameter object

# Returns:
- TESS observations catalog containing TESS targets and associated TOIs (to be used as true catalog in comparison with simulated observations)
"""
function setup_actual_pc_catalog_tess(df_star::DataFrame, df_toi::DataFrame, usable_toi::Array{Int64}, sim_param::SimParam)
    local target_obs, num_pl
    df_toi = df_toi[usable_toi,:]
    output = TESSObsCatalog()
    sort!(df_star, (:ticid))
    df_obs = innerjoin(select(df_star, Not(:sectors)), df_toi, on=:ticid => :TIC, makeunique=false, validate=(false, false))

    df_obs = sort!(df_obs, (:ticid))

    # if haskey(sim_param, "toi_subset_csv")
    #     tot_plan -= length(df_obs[!,:kepoi_name])
    #     println("# Number of planet candidates in subset file with no matching star in table: ", tot_plan)
    # end

    # Add each TOI planet candidate to TESS target object
    plid = 0
    for i in 1:length(df_obs[!,:toi_id])
        if plid == 0
            plid = 1
            while i+plid < length(df_obs[!,:toi_id]) && df_obs[i+plid,:ticid] == df_obs[i,:ticid]
                plid += 1
            end
            num_pl = plid
            target_obs = TESSTargetObs(num_pl)
            star_idx = searchsortedfirst(df_star[!,:ticid],df_obs[i,:ticid])
            if star_idx > length(df_star[!,:ticid])
                @warn "# Couldn't find TIC " * df_star[i,:ticid] * " in df_obs."
                star_idx = rand(1:length(df_star[!,:ticid]))
            end
            target_obs.star = ExoplanetsSysSim.StarObs(df_obs[i,:toi_prad],df_obs[i,:mass],star_idx)

        end

        target_obs.obs[plid] = ExoplanetsSysSim.TransitPlanetObs(df_obs[i,:toi_period],df_obs[i,:epoch],df_obs[i,:toi_transit_depth]/1.0e6,df_obs[i,:toi_transit_dur])
        target_obs.sigma[plid] = ExoplanetsSysSim.TransitPlanetObs(abs(df_obs[i,:toi_period_err]), abs(df_obs[i,:epoch_err]),abs(df_obs[i,:toi_transit_depth_err]/1.0e6),abs(df_obs[i,:toi_transit_dur_err]))
	#target_obs.prob_detect = ExoplanetsSysSim.SimulatedSystemDetectionProbs{OneObserver}( ones(num_pl), ones(num_pl,num_pl), ones(num_pl), fill(Array{Int64}(undef,0), 1) )  # Made line below to simplify calling
        target_obs.prob_detect = ExoplanetsSysSim.OneObserverSystemDetectionProbs(num_pl)
        plid -= 1
        if plid == 0
            push!(output.target,target_obs)
        end
    end
    return output
end

# Two functions below were just for debugging purposes
# Calculate SNR of every planet in simulated catalog
function calc_snr_list(cat::TESSPhysicalCatalog, sim_param::SimParam)
  snrlist = Array{Float64}(undef,0)
  for t in 1:length(cat.target)
    for p in 1:length(cat.target[t].sys[1].planet)
      snr = calc_snr_if_transit(cat.target[t],1,p,sim_param)
      if snr>0.0
        push!(snrlist,snr)
      end
    end
  end
  snrlist[findall(x->x>7.1,snrlist)]
end

# Calculate detection probability (assuming planet transits) for every planet in simulated catalog
function calc_prob_detect_list(cat::TESSPhysicalCatalog, sim_param::SimParam)
  pdetectlist = Array{Float64}(undef,0)
  for t in 1:length(cat.target)
    for p in 1:length(cat.target[t].sys[1].planet)
      #pdet = calc_prob_detect_if_transit(cat.target[t],1,p,sim_param)
      pdet = calc_prob_detect_if_transit_with_actual_b(cat.target[t],1,p,sim_param)

      if pdet>0.0
        push!(pdetectlist,pdet)
      end
    end
  end
  idx = findall(x->x>0.0,pdetectlist)
  pdetectlist[idx]
end