text:
using CImGui
using CImGui: ImVec2, ImVec4, IM_COL32, ImU32
using CSyntax
using CSyntax.CFor
using CSyntax.CStatic

"""
    ShowExampleAppCustomRendering(p_open::Ref{Bool})
Demonstrate using the low-level ImDrawList to draw custom shapes.
"""
function ShowExampleAppCustomRendering(p_open::Ref{Bool})
    CImGui.SetNextWindowSize((350, 560), CImGui.ImGuiCond_FirstUseEver)
    CImGui.Begin("Example: Custom rendering", p_open) || (CImGui.End(); return)

    draw_list = CImGui.GetWindowDrawList()

    # primitives
    CImGui.Text("Primitives")
    sz, thickness, col = @cstatic sz=Cfloat(36.0) thickness=Cfloat(4.0) col=Cfloat[1.0,1.0,0.4,1.0] begin
        @c CImGui.DragFloat("Size", &sz, 0.2, 2.0, 72.0, "%.0f")
        @c CImGui.DragFloat("Thickness", &thickness, 0.05, 1.0, 8.0, "%.02f")
        CImGui.ColorEdit4("Color", col)
    end

    p = CImGui.GetCursorScreenPos()
    col32 = CImGui.ColorConvertFloat4ToU32(ImVec4(col...))
    begin
        x::Cfloat = p.x + 4.0
        y::Cfloat = p.y + 4.0
        spacing = 8.0
        for n = 0:2-1
            th::Cfloat = (n == 0) ? 1.0 : thickness
            CImGui.AddCircle(draw_list, ImVec2(x+sz*0.5, y+sz*0.5), sz*0.5, col32, 6, th); x += sz + spacing; # hexagon
            CImGui.AddCircle(draw_list, ImVec2(x+sz*0.5, y+sz*0.5), sz*0.5, col32, 20, th); x += sz + spacing; # circle
            CImGui.AddRect(draw_list, ImVec2(x, y), ImVec2(x+sz, y+sz), col32, 0.0, CImGui.ImDrawCornerFlags_All, th); x += sz + spacing;
            CImGui.AddRect(draw_list, ImVec2(x, y), ImVec2(x+sz, y+sz), col32, 10.0, CImGui.ImDrawCornerFlags_All, th); x += sz + spacing;
            CImGui.AddRect(draw_list, ImVec2(x, y), ImVec2(x+sz, y+sz), col32, 10.0, CImGui.ImDrawCornerFlags_TopLeft|CImGui.ImDrawCornerFlags_BotRight, th); x += sz + spacing;
            CImGui.AddTriangle(draw_list, ImVec2(x+sz*0.5, y), ImVec2(x+sz,y+sz-0.5), ImVec2(x,y+sz-0.5), col32, th); x += sz + spacing;
            CImGui.AddLine(draw_list, ImVec2(x, y), ImVec2(x+sz,    y), col32, th); x += sz + spacing;  # horizontal line (note: drawing a filled rectangle will be faster!)
            CImGui.AddLine(draw_list, ImVec2(x, y), ImVec2(x,    y+sz), col32, th); x += spacing;       # vertical line (note: drawing a filled rectangle will be faster!)
            CImGui.AddLine(draw_list, ImVec2(x, y), ImVec2(x+sz, y+sz), col32, th); x += sz +spacing;   # diagonal line
            CImGui.AddBezierCurve(draw_list, ImVec2(x, y), ImVec2(x+sz*1.3,y+sz*0.3), (x+sz-sz*1.3,y+sz-sz*0.3), ImVec2(x+sz, y+sz), col32, th);
            x = p.x + 4
            y += sz + spacing
        end
        CImGui.AddCircleFilled(draw_list, ImVec2(x+sz*0.5, y+sz*0.5), sz*0.5, col32, 6); x += sz+spacing; # hexagon
        CImGui.AddCircleFilled(draw_list, ImVec2(x+sz*0.5, y+sz*0.5), sz*0.5, col32, 32); x += sz+spacing; # circle
        CImGui.AddRectFilled(draw_list, ImVec2(x, y), ImVec2(x+sz, y+sz), col32); x += sz+spacing;
        CImGui.AddRectFilled(draw_list, ImVec2(x, y), ImVec2(x+sz, y+sz), col32, 10.0); x += sz+spacing;
        CImGui.AddRectFilled(draw_list, ImVec2(x, y), ImVec2(x+sz, y+sz), col32, 10.0, CImGui.ImDrawCornerFlags_TopLeft|CImGui.ImDrawCornerFlags_BotRight); x += sz+spacing;
        CImGui.AddTriangleFilled(draw_list, ImVec2(x+sz*0.5, y), ImVec2(x+sz,y+sz-0.5), ImVec2(x,y+sz-0.5), col32); x += sz+spacing;
        CImGui.AddRectFilled(draw_list, ImVec2(x, y), ImVec2(x+sz, y+thickness), col32); x += sz+spacing;          # horizontal line (faster than AddLine, but only handle integer thickness)
        CImGui.AddRectFilled(draw_list, ImVec2(x, y), ImVec2(x+thickness, y+sz), col32); x += spacing+spacing;     # vertical line (faster than AddLine, but only handle integer thickness)
        CImGui.AddRectFilled(draw_list, ImVec2(x, y), ImVec2(x+1, y+1), col32);          x += sz;                  # pixel (faster than AddLine)
        CImGui.AddRectFilledMultiColor(draw_list, ImVec2(x, y), ImVec2(x+sz, y+sz), IM_COL32(0,0,0,255), IM_COL32(255,0,0,255), IM_COL32(255,255,0,255), IM_COL32(0,255,0,255))
        CImGui.Dummy(ImVec2((sz+spacing)*8, (sz+spacing)*3))
    end
    CImGui.Separator()
    @cstatic adding_line=false points=ImVec2[] begin
        CImGui.Text("Canvas example")
        CImGui.Button("Clear") && empty!(points)
        if length(points) ≥ 2
            CImGui.SameLine()
            CImGui.Button("Undo") && (pop!(points); pop!(points);)
        end
        CImGui.Text("Left-click and drag to add lines,\nRight-click to undo")

        # here we are using InvisibleButton() as a convenience to 1) advance the cursor and 2) allows us to use IsItemHovered()
        # but you can also draw directly and poll mouse/keyboard by yourself. You can manipulate the cursor using GetCursorPos() and SetCursorPos().
        # if you only use the ImDrawList API, you can notify the owner window of its extends by using SetCursorPos(max).
        canvas_pos = CImGui.GetCursorScreenPos()            # ImDrawList API uses screen coordinates!
        canvas_size = CImGui.GetContentRegionAvail()        # resize canvas to what's available

        cx, cy = canvas_size.x, canvas_size.y
        cx < 50.0 && (cx = 50.0)
        cy < 50.0 && (cy = 50.0)
        canvas_size = ImVec2(cx, cy)
        CImGui.AddRectFilledMultiColor(draw_list, canvas_pos, ImVec2(canvas_pos.x + canvas_size.x, canvas_pos.y + canvas_size.y), IM_COL32(50, 50, 50, 255), IM_COL32(50, 50, 60, 255), IM_COL32(60, 60, 70, 255), IM_COL32(50, 50, 60, 255))
        CImGui.AddRect(draw_list, canvas_pos, ImVec2(canvas_pos.x + canvas_size.x, canvas_pos.y + canvas_size.y), IM_COL32(255, 255, 255, 255))

        adding_preview = false
        CImGui.InvisibleButton("canvas", canvas_size)
        mouse_pos_in_canvas = ImVec2(CImGui.GetIO().MousePos.x - canvas_pos.x, CImGui.GetIO().MousePos.y - canvas_pos.y)
        if adding_line
            adding_preview = true
            push!(points, mouse_pos_in_canvas)
            !CImGui.IsMouseDown(0) && (adding_line = adding_preview = false;)
        end
        if CImGui.IsItemHovered()
            if !adding_line && CImGui.IsMouseClicked(0)
                push!(points, mouse_pos_in_canvas)
                adding_line = true
            end
            if CImGui.IsMouseClicked(1) && !isempty(points)
                adding_line = adding_preview = false
                pop!(points)
                pop!(points)
            end
        end
        CImGui.PushClipRect(draw_list, canvas_pos, ImVec2(canvas_pos.x + canvas_size.x, canvas_pos.y + canvas_size.y), true) # clip lines within the canvas (if we resize it, etc.)
        @cfor i=1 i<length(points) i+=2 begin
            CImGui.AddLine(draw_list, ImVec2(canvas_pos.x + points[i].x, canvas_pos.y + points[i].y), ImVec2(canvas_pos.x + points[i + 1].x, canvas_pos.y + points[i + 1].y), IM_COL32(255, 255, 0, 255), 2.0)
        end
        CImGui.PopClipRect(draw_list)
        adding_preview && pop!(points)
    end
    CImGui.End()
end


text:
bshow(i) = bitstring(i)[end-7:end]

function large_enough_unsigned(bit_cnt)
    unsigned_types = [UInt8, UInt16, UInt32, UInt64, UInt128]
    atype = nothing
    for xtype in unsigned_types
        if sizeof(xtype) * 8 >= bit_cnt
            atype = xtype
            break
        end
    end
    atype
end

"""
    is_power_of_two(v::Base.BitInteger)

This integer is a power of two.
"""
function is_power_of_two(v::Base.BitInteger)
    v != 0 && ((v & (v - 1)) == 0)
end

"""
    log_base2(v::Integer)

Count the number of bits in an integer, rounding down.
"""
function log_base2(v::Integer)::Integer
    r = zero(v)
    while (v = v >> one(v)) > 0
      r += one(v)
    end
    r
end

"""
    count_set_bits(v::Integer)

A naive algorithm to count bits that are set.
Look here to improve: https://graphics.stanford.edu/~seander/bithacks.html.
"""
function count_set_bits(v::Integer)
    c = zero(v)
    while v != zero(v)
        c += v & one(v)
        v = v>>one(v)
    end
    c
end

"""
    msb(v::Integer)

Most-significant bit, zero-based count. So `0b1` is 0,
`0b1010` is 3.
"""
function msb(v::Integer)
    r = 0
    while (v >>= 1) != 0
        r += 1
    end
    r
end

"""
    trailing_zero_bits(v::Integer)

The number of zero bits after the last one-bit.
"""
function trailing_zero_bits(v::Integer)
    c = zero(v)
    if v != zero(v)
        # Set v's trailing 0s to 1s and zero rest
        v = (v ⊻ (v - one(v))) >> one(v)
        while v != zero(v)
            v = v >> 1
            c += one(v)
        end
    else
        c = 8 * sizeof(v)
    end
    c
end

"""
    trailing_set_bits(v::Integer)

The number of one-bits after the last zero-bit.
"""
trailing_set_bits(v) = trailing_zero_bits(~v)

"""
Treat `x` as an `n`-bit unsigned integer. Rotate the bits
`k` places to the left, wrapping them around the right side.
"""
function rotateleft_naive(x::Base.BitInteger, k::Integer, n::Integer)
    y = zero(x)
    for i = 0:(n - 1)
        ind = (i - k) % n
        ind = (ind >= 0) ? ind : ind + n
        if (x & (one(x) << ind)) > 0
            y |= (one(x) << i)
        end  # else leave unset.
    end
    y
end

function rotateleft(x::T, k::Integer, n::Integer) where {T <: Base.BitInteger}
    @assert k >= 0
    @assert n > 0
    @assert k < n
    x &= fbvn1s(T, n)
    x = (x << k) | (x >> (n - k))
    x &= fbvn1s(T, n)
    x
end

function rotateright_naive(x::Base.BitInteger, k::Integer, n::Integer)
    y = zero(x)
    for i = 0:(n - 1)
        ind = (i + k) % n
        ind = (ind >= 0) ? ind : ind + n
        if (x & (one(x) << ind)) > 0
            y |= (one(x) << i)
        end  # else leave unset.
    end
    y
end

function rotateright(x::T, k::Integer, n::Integer) where {T <: Base.BitInteger}
    x &= fbvn1s(T, n)
    x = (x >> k) | (x << (n - k))
    x &= fbvn1s(T, n)
    x
end

"""
2^n - 1, computed for an unsigned integer.
"""
function fbvn1s(T::DataType, n)
    if n == T(sizeof(T) * 8)
        ~zero(T)
    else
        (one(T) << n) - one(T)
    end
end

"""
2^n - 1, computed for an unsigned integer.
"""
function fbvn1s(_::T, n) where {T}
    if n == T(sizeof(T) * 8)
        ~zero(T)
    else
        (one(T) << n) - one(T)
    end
end

function fbvmod(i, m)
    if i >= m
        i -= m * i / m
    else
        i
    end
end

function reverse_bits(w, n)
    r = zero(w)
    for i in 0:(n - 1)
        r |= ((w & (one(w)<<i)) >> i) << (n - 1 - i)
    end
    r
end

function set_bits(h, n, i, w)
    h |= w << (i * n)
    h
end

function set_bits_naive(h, n, i, w)
    for j in 0:(n - 1)
        b = (w & (1 << j)) >> j
        h |= (b << (n * i + j))
    end
    h
end

"""
Get b bits from v, starting at place i.
"""
function get_bits(v::T, b, i) where {T}
    n = 8 * sizeof(T)
    (v << (n - b - i)) >> (n - b)
end

function trailing_set_bits_hamilton(i::UInt64)
    T = UInt64
    c = zero(Int)
    if i == ~zero(T)
        return T(8) * sizeof(T)
    end
    if (i & fbvn1s(T, 32)) == fbvn1s(T, 32)
        i >>= T(32)
        c ⊻= T(32)
    end
    if (i & fbvn1s(T, 16)) == fbvn1s(T, 16)
        i >>= T(16)
        c ⊻= T(16)
    end
    if (i & fbvn1s(T, 8)) == fbvn1s(T, 8)
        i >>= T(8)
        c ⊻= T(8)
    end
    if (i & fbvn1s(T, 4)) == fbvn1s(T, 4)
        i >>= T(4)
        c ⊻= T(4)
    end
    if (i & fbvn1s(T, 2)) == fbvn1s(T, 2)
        i >>= T(2)
        c ⊻= T(2)
    end
    if (i & fbvn1s(T, 1)) == fbvn1s(T, 1)
        i >>= T(1)
        c ⊻= T(1)
    end
    c
end

function trailing_set_bits_hamilton(i::UInt32)
    T = UInt32
    c = zero(Int)
    if i == ~zero(T)
        return T(8) * sizeof(T)
    end
    if (i & fbvn1s(T, 16)) == fbvn1s(T, 16)
        i >>= T(16)
        c ⊻= T(16)
    end
    if (i & fbvn1s(T, 8)) == fbvn1s(T, 8)
        i >>= T(8)
        c ⊻= T(8)
    end
    if (i & fbvn1s(T, 4)) == fbvn1s(T, 4)
        i >>= T(4)
        c ⊻= T(4)
    end
    if (i & fbvn1s(T, 2)) == fbvn1s(T, 2)
        i >>= T(2)
        c ⊻= T(2)
    end
    if (i & fbvn1s(T, 1)) == fbvn1s(T, 1)
        i >>= T(1)
        c ⊻= T(1)
    end
    c
end

function trailing_set_bits_hamilton(i::UInt16)
    T = UInt16
    c = zero(Int)
    if i == ~zero(T)
        return T(8) * sizeof(T)
    end
    if (i & fbvn1s(T, 8)) == fbvn1s(T, 8)
        i >>= T(8)
        c ⊻= T(8)
    end
    if (i & fbvn1s(T, 4)) == fbvn1s(T, 4)
        i >>= T(4)
        c ⊻= T(4)
    end
    if (i & fbvn1s(T, 2)) == fbvn1s(T, 2)
        i >>= T(2)
        c ⊻= T(2)
    end
    if (i & fbvn1s(T, 1)) == fbvn1s(T, 1)
        i >>= T(1)
        c ⊻= T(1)
    end
    c
end

function trailing_set_bits_hamilton(i::UInt8)
    T = UInt8
    c = zero(Int)
    if i == ~zero(T)
        return T(8) * sizeof(T)
    end
    if (i & fbvn1s(T, 4)) == fbvn1s(T, 4)
        i >>= T(4)
        c ⊻= T(4)
    end
    if (i & fbvn1s(T, 2)) == fbvn1s(T, 2)
        i >>= T(2)
        c ⊻= T(2)
    end
    if (i & fbvn1s(T, 1)) == fbvn1s(T, 1)
        i >>= T(1)
        c ⊻= T(1)
    end
    c
end

function first_set_bit(i::T)::Int where {T <: Integer}
    if i == 0
        return 0
    end
    for j = 0:(8 * sizeof(T)  - 1)
        if (i & (T(1) << j)) != 0
            return j + 1
        end
    end
end

function first_set_bit(i::UInt128)::Int
    T = UInt128
    c = zero(Int)
    if i == zero(T)
        return c
    end
    if (i & fbvn1s(T, 64)) == zero(T)
        i >>= T(64)
        c ⊻= 64
    end
    if (i & fbvn1s(T, 32)) == zero(T)
        i >>= T(32)
        c ⊻= 32
    end
    if (i & fbvn1s(T, 16)) == zero(T)
        i >>= T(16)
        c ⊻= 16
    end
    if (i & fbvn1s(T, 8)) == zero(T)
        i >>= T(8)
        c ⊻= 8
    end
    if (i & fbvn1s(T, 4)) == zero(T)
        i >>= T(4)
        c ⊻= 4
    end
    if (i & fbvn1s(T, 2)) == zero(T)
        i >>= T(2)
        c ⊻= 2
    end
    if (i & fbvn1s(T, 1)) == zero(T)
        i >>= T(1)
        c ⊻= 1
    end
    c + one(Int)
end

function first_set_bit(i::UInt64)::Int
    T = UInt64
    c = zero(Int)
    if i == zero(T)
        return c
    end
    if (i & fbvn1s(T, 32)) == zero(T)
        i >>= T(32)
        c ⊻= 32
    end
    if (i & fbvn1s(T, 16)) == zero(T)
        i >>= T(16)
        c ⊻= 16
    end
    if (i & fbvn1s(T, 8)) == zero(T)
        i >>= T(8)
        c ⊻= 8
    end
    if (i & fbvn1s(T, 4)) == zero(T)
        i >>= T(4)
        c ⊻= 4
    end
    if (i & fbvn1s(T, 2)) == zero(T)
        i >>= T(2)
        c ⊻= 2
    end
    if (i & fbvn1s(T, 1)) == zero(T)
        i >>= T(1)
        c ⊻= 1
    end
    c + one(Int)
end

function first_set_bit(i::UInt32)::Int
    T = UInt32
    c = zero(Int)
    if i == zero(T)
        return c
    end
    if (i & fbvn1s(T, 16)) == zero(T)
        i >>= T(16)
        c ⊻= 16
    end
    if (i & fbvn1s(T, 8)) == zero(T)
        i >>= T(8)
        c ⊻= 8
    end
    if (i & fbvn1s(T, 4)) == zero(T)
        i >>= T(4)
        c ⊻= 4
    end
    if (i & fbvn1s(T, 2)) == zero(T)
        i >>= T(2)
        c ⊻= 2
    end
    if (i & fbvn1s(T, 1)) == zero(T)
        i >>= T(1)
        c ⊻= 1
    end
    c + one(Int)
end

function first_set_bit(i::UInt16)::Int
    T = UInt16
    c = zero(Int)
    if i == zero(T)
        return c
    end
    if (i & fbvn1s(T, 8)) == zero(T)
        i >>= T(8)
        c ⊻= 8
    end
    if (i & fbvn1s(T, 4)) == zero(T)
        i >>= T(4)
        c ⊻= 4
    end
    if (i & fbvn1s(T, 2)) == zero(T)
        i >>= T(2)
        c ⊻= 2
    end
    if (i & fbvn1s(T, 1)) == zero(T)
        i >>= T(1)
        c ⊻= 1
    end
    c + one(Int)
end

function first_set_bit(i::UInt8)::Int
    T = UInt8
    c = zero(Int)
    if i == zero(T)
        return c
    end
    if (i & fbvn1s(T, 4)) == zero(T)
        i >>= T(4)
        c ⊻= 4
    end
    if (i & fbvn1s(T, 2)) == zero(T)
        i >>= T(2)
        c ⊻= 2
    end
    if (i & fbvn1s(T, 1)) == zero(T)
        i >>= T(1)
        c ⊻= 1
    end
    c + one(Int)
end

function getloc_case(p, n, i, l, case)
    if case == 1
        if (n <= i + 1) && (p[i] & im != 0)
            l | (1 << i)
        else
            l
        end
    else
        l = getloc_case(p, n, i + case >> 1, l, case >> 1)
        getloc_case(p, n, i, l, case >> 1)
    end
end

function get_location(p::Vector{T}, n, i) where {T}
    im = one(T) << i
    # p, jo=0, jn=n, ir=0, im=1<<i, l is return
    l = 0
    getloc_case(p, n, 0, l, 8 * sizeof(T))
end

ext: jl


text:
function CC.OptimizationState(linfo::MethodInstance, params::OptimizationParams, interp::JETInterpreter)
    return OptimizationState(linfo, params, interp.native)
end

function CC.optimize(interp::JETInterpreter, opt::OptimizationState, params::OptimizationParams, result)
    # NOTE: don't recurse with `@invoke`
    # `JETInterpreter` shouldn't recur into the optimization step via `@invoke` macro,
    # but rather it should just go through `NativeInterpreter`'s optimization pass
    # this is necessary because `CC.get(wvc::WorldView{JETGlobalCache}, mi::MethodInstance, default)`
    # can be called from optimization pass otherwise, and it may cause errors because our report
    # construction is only valid on inference frames before the optimizer runs on it
    return optimize(interp.native, opt, params, result)
end


text:
export append, append!

"""
    append!(data::PopData, data2::PopData)
Add the rows of `data2` to the end of `data`. This will add the samples present
in the second `PopData` object to the first `PopData` object (mutating it).
**Note** that this is a simple appending, and you risk corrupting your `PopData` if
the two `PopData` objects do not have identical loci.

**Example**
```
julia> cats = @nancycats
PopData{Diploid, 9 Microsatellite Loci}
  Samples: 237
  Populations: 17

julia> purrfect_pairs = cross(cats, "N200", "N7", generation = "F1")
PopData{Diploid, 9 Microsatellite Loci}
  Samples: 100
  Populations: 1

julia> append!(cats, purrfect_pairs);

julia> cats
PopData{Diploid, 9 Microsatellite Loci}
  Samples: 337
  Populations: 18
```
"""
function Base.append!(data::PopData, data2::PopData)
  n1 = data.metadata.samples
  pl1 = data.metadata.ploidy
  if "parents" ∉ names(data.sampleinfo) && "parents" ∈ names(data2.sampleinfo)
      len = data.metadata.samples
      insertcols!(
          data.sampleinfo,
          :parents => Vector{Union{Missing, Tuple{String,String}}}(undef, len)
          )
  elseif "parents" ∉ names(data2.sampleinfo) && "parents" ∈ names(data.sampleinfo)
      len = length(data2.sampleinfo.name)
      insertcols!(
          data2.sampleinfo,
          :parents => Vector{Union{Missing, Tuple{String,String}}}(undef, len)
          )
  end

  append!(data.sampleinfo, data2.sampleinfo)

  append!(data.genodata, data2.genodata)
  # update metadata
  data.metadata.samples = n1 + data2.metadata.samples
  if pl1 != data2.metadata.ploidy
    data.metadata.ploidy = Int8[pl1, data2.metadata.ploidy]
  end
  data.metadata.populations = length(unique(data.sampleinfo.population))
  return data
end

"""
    append(data::PopData, data2::PopData)
Add the rows of `data2` to the end of `data`. This will combine the samples present
in both `PopData` objects and return a new `PopData` object. **Note** that this is
a simple appending, and you risk corrupting your `PopData` if the two `PopData`
objects do not have identical loci.

**Example**
```
julia> cats = @nanycats
PopData{Diploid, 9 Microsatellite Loci}
  Samples: 237
  Populations: 17

julia> purrfect_pairs = cross(cats, "N200", "N7", generation = "F1")
PopData{Diploid, 9 Microsatellite Loci}
  Samples: 100
  Populations: 1

julia> merged_cats = append(cats, purrfect_pairs)
PopData{Diploid, 9 Microsatellite Loci}
  Samples: 337
  Populations: 18
```
"""
function append(data::PopData, data2::PopData)
  #=
  new_pdi = deepcopy(data.info)
  append!(new_pdi.sampleinfo, data2.sampleinfo)
  new_pdi.samples = data.info.samples + data2.info.samples
  new_pdi.populations = length(unique(new_pdi.sampleinfo.population))
  ploid = vcat(data.metadata.ploidy, data2.metadata.ploidy)
  new_pdi.ploidy = length(unique(ploid)) == 1 ? ploid[1] : Int8[unique(ploidy)]
  PopData(new_pdi, vcat(data.genodata, data2.genodata))
  =#
  tmp  = copy(data)
  append!(tmp, data2)
  return tmp
end

function allele_pool(locus::T) where T <: GenoArray
  Tuple(Base.Iterators.flatten(skipmissing(locus)))
end

function allele_pool(data::PopData)
  # index dataframe by locus
  idx_df = groupby(data.genodata, [:locus])
  # instantiate dict to store alleles
  # pull out loci names
  loc = getindex.(keys(idx_df), :locus)
  allele_dict = Dict(i => allele_pool(idx_df[(;locus = i)].genotype) for i in loc)
  return string.(loc), allele_dict
end

"""
```
simulate_sample(alleles::Dict{String,NTuple}, loc::Vector{String}; ploidy::Int)
```
Using a global allele pool given by a Dict{loci,alleles} and a list of loci (`loc`), simulate
an individual with a given `ploidy`. Returns a Vector of genotypes.

**Example**
```
julia> cats = @nanycats ;
julia> loc, alleles = allele_pool(cats) ;
julia> simulate_sample(alleles, loc, ploidy = 2)
9-element Array{Array{Int16,1},1}:
 [139, 129]
 [146, 146]
 [145, 141]
 [126, 126]
 [150, 148]
 [148, 140]
 [185, 199]
 [91, 113]
 [208, 208]
```
"""
function simulate_sample(alleles::Dict{String,<:Tuple}, loc::Vector{String}; ploidy::Signed)
  map(i -> rand(alleles[i], ploidy), loc)
end

feature_req() = "\nPlease open an Issue or Pull Request on https://www.github.com/pdimens/PopGenSims.jl if you would like this feature implemented"


text:
using Test
using StaticArrays, IntervalSets
import ClimaCore.DataLayouts: IJFH
import ClimaCore:
    Fields, Domains, Meshes, Topologies, Spaces, Operators, Geometry
using StaticArrays, IntervalSets, LinearAlgebra

FT = Float64

radius = FT(1)
domain = Domains.SphereDomain(radius)

grad = Operators.Gradient()
wgrad = Operators.WeakGradient()

@testset "gradient" begin
    Ne = 5
    Nq = 6

    mesh = Meshes.EquiangularCubedSphere(domain, Ne)
    grid_topology = Topologies.Topology2D(mesh)

    quad = Spaces.Quadratures.GLL{Nq}()
    space = Spaces.SpectralElementSpace2D(grid_topology, quad)
    coords = Fields.coordinate_field(space)

    η = @. sind(coords.long) * cosd(coords.lat)
    ∇η = grad.(η)
    ∇η_w = wgrad.(η)
    ∇η_local = Geometry.transform.(Ref(Geometry.UVAxis()), ∇η)
    ∇ηw_local = Geometry.transform.(Ref(Geometry.UVAxis()), ∇η_w)
    Spaces.weighted_dss!(∇η_local)
    Spaces.weighted_dss!(∇ηw_local)

    ∇η_exact = @. Geometry.UVVector(
        cosd(coords.long) / radius,
        -sind(coords.long) * sind(coords.lat) / radius,
    )

    @test ∇η_local ≈ ∇η_exact rtol = 1e-2
    @test ∇ηw_local ≈ ∇η_exact rtol = 1e-2

end

"""
    convergence_rate(err, Δh)
Estimate convergence rate given vectors `err` and `Δh`
    err = C Δh^p+ H.O.T
    err_k ≈ C Δh_k^p
    err_k/err_m ≈ Δh_k^p/Δh_m^p
    log(err_k/err_m) ≈ log((Δh_k/Δh_m)^p)
    log(err_k/err_m) ≈ p*log(Δh_k/Δh_m)
    log(err_k/err_m)/log(Δh_k/Δh_m) ≈ p
"""
convergence_rate(err, Δh) =
    [log(err[i] / err[i - 1]) / log(Δh[i] / Δh[i - 1]) for i in 2:length(Δh)]

@testset "convergence tests for the gradient operator on the sphere" begin
    Nes = [3, 9, 27]
    Nqs = [4, 6]

    for (Iq, Nq) in enumerate(Nqs)
        err_sgrad = Array{FT}(undef, length(Nes))
        err_wgrad = Array{FT}(undef, length(Nes))
        Δh = Array{FT}(undef, length(Nes))

        for (Ie, Ne) in enumerate(Nes)
            mesh = Meshes.EquiangularCubedSphere(domain, Ne)
            grid_topology = Topologies.Topology2D(mesh)

            quad = Spaces.Quadratures.GLL{Nq}()
            space = Spaces.SpectralElementSpace2D(grid_topology, quad)
            coords = Fields.coordinate_field(space)

            η = @. sind(coords.long) * cosd(coords.lat)
            ∇η = grad.(η)
            ∇η_w = wgrad.(η)
            ∇η_local = Geometry.transform.(Ref(Geometry.UVAxis()), ∇η)
            ∇ηw_local = Geometry.transform.(Ref(Geometry.UVAxis()), ∇η_w)
            Spaces.weighted_dss!(∇η_local)
            Spaces.weighted_dss!(∇ηw_local)

            ∇η_exact = @. Geometry.UVVector(
                cosd(coords.long) / radius,
                -sind(coords.long) * sind(coords.lat) / radius,
            )

            err_sgrad[Ie] = norm(∇η_local .- ∇η_exact)
            err_wgrad[Ie] = norm(∇ηw_local .- ∇η_exact)
            Δh[Ie] = 1 / Ne

        end

        convergence_rate_sgrad = convergence_rate(err_sgrad, Δh)
        convergence_rate_wgrad = convergence_rate(err_wgrad, Δh)

        for Ie in range(1, length = length(Nes) - 1)
            @test convergence_rate_sgrad[Ie] ≈ (Nq - 1) atol = 0.1
            @test convergence_rate_wgrad[Ie] ≈ (Nq - 1) atol = 0.1
        end

    end

end


text:
importall ImmutableArrays

typealias GLboolean  Bool
typealias GLbyte     Int8
typealias GLubyte    Uint8
typealias GLchar     Uint8
typealias GLshort    Int16
typealias GLushort   Uint16
typealias GLint      Int32
typealias GLuint     Uint32
typealias GLint64    Int64
typealias GLuint64   Uint64
typealias GLsizei    Uint32
typealias GLenum     Uint32
typealias GLintptr   Ptr
typealias GLsizeiptr Int
typealias GLfloat    Float32
typealias GLdouble   Float64
typealias GLvoid     Void

immutable Object
	handle::GLuint
end
Object(x) = Object(convert(GLuint,x))

immutable Shader
	handle::Object
end
Shader(x) = Shader(Object(x))

immutable Program
	handle::Object
end
Program(x) = Program(Object(x))

typealias Buffer      Object
typealias Texture     Object
typealias VertexArray Object

typealias Attribute GLint

immutable Uniform
	location::GLint
end
Uniform(location) = Uniform(convert(GLint, location))
+(u::Uniform, i::Integer) = Uniform(u.location+i)

typealias Vec2 Vector2{Float32}
typealias Vec3 Vector3{Float32}
typealias Vec4 Vector4{Float32}


text:
module Nodal
using Dates, LinearAlgebra, Mmap, SeisIO, SeisIO.FastIO
path = Base.source_dir()

# Imports
include("Nodal/imports.jl")

# Constants
include("Nodal/Types/TDMSbuf.jl")
include("Nodal/constants.jl")

# Types
include("Nodal/Types/NodalData.jl")
include("Nodal/Types/NodalChannel.jl")

# Formats
for i in ls(path*"/Nodal/Formats/")
  if endswith(i, ".jl")
    include(i)
  end
end

# Utils
for i in ls(path*"/Nodal/Utils/")
  if endswith(i, ".jl")
    include(i)
  end
end

# Wrappers
include("Nodal/Wrappers/read_nodal.jl")

# Exports
include("Nodal/exports.jl")

end


text:
module LinearDynamicsModels

using LinearAlgebra
using StaticArrays
using DifferentialDynamicsModels
using ForwardDiff
using Requires
using MacroTools

import DifferentialDynamicsModels: SteeringBVP
import DifferentialDynamicsModels: state_dim, control_dim, duration, propagate, instantaneous_control
export LinearDynamics, ZeroOrderHoldLinearization, FirstOrderHoldLinearization, linearize
export NIntegratorDynamics, DoubleIntegratorDynamics, TripleIntegratorDynamics
export LinearQuadraticSteering, NIntegratorSteering, DoubleIntegratorSteering, TripleIntegratorSteering

include("utils.jl")

# Continous-Time Linear Time-Invariant Systems
struct LinearDynamics{Dx,Du,TA<:StaticMatrix{Dx,Dx},TB<:StaticMatrix{Dx,Du},Tc<:StaticVector{Dx}} <: DifferentialDynamics
    A::TA
    B::TB
    c::Tc
end
Base.zero(::Type{LinearDynamics{Dx,Du,TA,TB,Tc}}) where {Dx,Du,TA,TB,Tc} = LinearDynamics(zero(TA), zero(TB), zero(Tc))

state_dim(::LinearDynamics{Dx,Du}) where {Dx,Du} = Dx
control_dim(::LinearDynamics{Dx,Du}) where {Dx,Du} = Du

(f::LinearDynamics{Dx,Du})(x::StaticVector{Dx}, u::StaticVector{Du}) where {Dx,Du} = f.A*x + f.B*u + f.c
function propagate(f::LinearDynamics{Dx,Du}, x::StaticVector{Dx}, SC::StepControl{Du}) where {Dx,Du}
    y = f.B*SC.u + f.c
    eᴬᵗ, ∫eᴬᵗy = integrate_expAt_B(f.A, y, SC.t)
    eᴬᵗ*x + ∫eᴬᵗy
end
function propagate(f::LinearDynamics{Dx,Du}, x::StaticVector{Dx}, RC::RampControl{Du}) where {Dx,Du}
    y = f.B*RC.uf + f.c
    eᴬᵗ, ∫eᴬᵗy = integrate_expAt_B(f.A, y, RC.t)
    z = f.B*(RC.u0 - RC.uf)
    _, _, ∫eᴬᵗztdt⁻¹ = integrate_expAt_Bt_dtinv(f.A, z, RC.t)
    eᴬᵗ*x + ∫eᴬᵗy + ∫eᴬᵗztdt⁻¹
end

# Discrete-Time Linear Time-Invariant Systems
include("linearization.jl")

# NIntegrators (DoubleIntegrator, TripleIntegrator, etc.)
function NIntegratorDynamics(::Val{N}, ::Val{D}, ::Type{T} = Rational{Int}) where {N,D,T}
    A = diagm(Val(D) => ones(SVector{(N-1)*D,T}))
    B = [zeros(SMatrix{(N-1)*D,D,T}); SMatrix{D,D,T}(I)]
    c = zeros(SVector{N*D,T})
    LinearDynamics(A, B, c)
end
NIntegratorDynamics(N::Int, D::Int, ::Type{T} = Rational{Int}) where {T} = NIntegratorDynamics(Val(N), Val(D), T)
DoubleIntegratorDynamics(D::Int, ::Type{T} = Rational{Int}) where {T} = NIntegratorDynamics(2, D, T)
TripleIntegratorDynamics(D::Int, ::Type{T} = Rational{Int}) where {T} = NIntegratorDynamics(3, D, T)

# TimePlusQuadraticControl BVPs
function SteeringBVP(f::LinearDynamics{Dx,Du}, j::TimePlusQuadraticControl{Du};
                     compile::Union{Val{false},Val{true}}=Val(false)) where {Dx,Du}
    compile === Val(true) ? error("Run `import SymPy` to enable SteeringBVP compilation.") :
                            SteeringBVP(f, j, EmptySteeringConstraints(), EmptySteeringCache())
end
const LinearQuadraticSteering{Dx,Du,Cache} = SteeringBVP{<:LinearDynamics{Dx,Du},<:TimePlusQuadraticControl{Du},EmptySteeringConstraints,Cache}
function LinearQuadraticSteering(A, B, c, R; compile::Union{Val{false},Val{true}}=Val(false))
    SteeringBVP(LinearDynamics(A, B, c), TimePlusQuadraticControl(R), compile=compile)
end
function NIntegratorSteering(N::Int, D::Int, R=SMatrix{D,D,Rational{Int}}(I); compile::Union{Val{false},Val{true}}=Val(false))
    SteeringBVP(NIntegratorDynamics(N, D), TimePlusQuadraticControl(R), compile=compile)
end
function DoubleIntegratorSteering(D::Int, R=SMatrix{D,D,Rational{Int}}(I); compile::Union{Val{false},Val{true}}=Val(false))
    SteeringBVP(DoubleIntegratorDynamics(D), TimePlusQuadraticControl(R), compile=compile)
end
function TripleIntegratorSteering(D::Int, R=SMatrix{D,D,Rational{Int}}(I); compile::Union{Val{false},Val{true}}=Val(false))
    SteeringBVP(TripleIntegratorDynamics(D), TimePlusQuadraticControl(R), compile=compile)
end

## Ad Hoc Steering
struct LinearQuadraticSteeringControl{Dx,Du,T,
                                      Tx0<:StaticVector{Dx},
                                      Txf<:StaticVector{Dx},
                                      Tf<:LinearDynamics{Dx,Du},
                                      Tj<:TimePlusQuadraticControl{Du},
                                      Tz<:StaticVector{Dx}} <: ControlInterval
    t::T
    x0::Tx0
    xf::Txf
    dynamics::Tf
    cost::Tj
    z::Tz
end
duration(lqsc::LinearQuadraticSteeringControl) = lqsc.t
function Base.zero(::Type{LinearQuadraticSteeringControl{Dx,Du,T,Tx0,Txf,Tf,Tj,Tz}}) where {Dx,Du,T,Tx0,Txf,Tf,Tj,Tz}
    LinearQuadraticSteeringControl(zero(T), zero(Tx0), zero(Txf), zero(Tf), zero(Tj), zero(Tz))
end
propagate(f::LinearDynamics, x::State, lqsc::LinearQuadraticSteeringControl) = (x - lqsc.x0) + lqsc.xf
function propagate(f::LinearDynamics, x::State, lqsc::LinearQuadraticSteeringControl, s::Number)
    f, j = lqsc.dynamics, lqsc.cost
    x0, A, B, c, R, z = lqsc.x0, f.A, f.B, f.c, j.R, lqsc.z
    eᴬˢ, ∫eᴬˢc = integrate_expAt_B(A, c, s)
    Gs = integrate_expAt_B_expATt(A, B*(R\B'), s)
    (x - x0) + eᴬˢ*x0 + ∫eᴬˢc + Gs*(eᴬˢ'\z)
end
function instantaneous_control(lqsc::LinearQuadraticSteeringControl, s::Number)
    A, B, R, z = lqsc.dynamics.A, lqsc.dynamics.B, lqsc.cost.R, lqsc.z
    eᴬˢ = exp(A*s)
    (R\B')*(eᴬˢ'\z)
end
function (j::TimePlusQuadraticControl)(lqsc::LinearQuadraticSteeringControl)
    @assert j == lqsc.cost
    cost(lqsc.dynamics, lqsc.cost, lqsc.x0, lqsc.xf, lqsc.t)
end

function (bvp::LinearQuadraticSteering{Dx,Du,EmptySteeringCache})(x0::StaticVector{Dx}, xf::StaticVector{Dx},
                                                                  c_max::T=eltype(x0)(1e6)) where {Dx,Du,T<:Number}    # TODO: handle c_max == Inf
    f, j = bvp.dynamics, bvp.cost
    A, B, c, R = f.A, f.B, f.c, j.R
    x0 == xf && return (cost=T(0), controls=LinearQuadraticSteeringControl(T(0), x0, xf, f, j, zeros(similar_type(typeof(c), T))))
    t = optimal_time(bvp, x0, xf, c_max)
    Q = B*(R\B')
    G = integrate_expAt_B_expATt(A, Q, t)
    eᴬᵗ, ∫eᴬᵗc = integrate_expAt_B(A, c, t)
    x̄ = eᴬᵗ*x0 + ∫eᴬᵗc
    z = eᴬᵗ'*(G\(xf - x̄))
    (cost=cost(f, j, x0, xf, t), controls=LinearQuadraticSteeringControl(t, x0, xf, f, j, z))
end

function cost(f::LinearDynamics{Dx,Du}, j::TimePlusQuadraticControl{Du},
              x0::StaticVector{Dx}, xf::StaticVector{Dx}, t) where {Dx,Du}
    A, B, c, R = f.A, f.B, f.c, j.R
    Q = B*(R\B')
    G = integrate_expAt_B_expATt(A, Q, t)
    eᴬᵗ, ∫eᴬᵗc = integrate_expAt_B(A, c, t)
    x̄ = eᴬᵗ*x0 + ∫eᴬᵗc
    t + (xf - x̄)'*(G\(xf - x̄))
end

function dcost(f::LinearDynamics{Dx,Du}, j::TimePlusQuadraticControl{Du},
               x0::StaticVector{Dx}, xf::StaticVector{Dx}, t) where {Dx,Du}
    A, B, c, R = f.A, f.B, f.c, j.R
    Q = B*(R\B')
    G = integrate_expAt_B_expATt(A, Q, t)
    eᴬᵗ, ∫eᴬᵗc = integrate_expAt_B(A, c, t)
    x̄ = eᴬᵗ*x0 + ∫eᴬᵗc
    z = eᴬᵗ'*(G\(xf - x̄))
    1 - 2*(A*x0 + c)'*z - z'*Q*z
end

function optimal_time(bvp::SteeringBVP{D,C,EmptySteeringConstraints,EmptySteeringCache},
                      x0::StaticVector{Dx},
                      xf::StaticVector{Dx},
                      t_max::T) where {Dx,Du,T<:Number,D<:LinearDynamics{Dx,Du},C<:TimePlusQuadraticControl{Du}}
    t = bisection(t -> dcost(bvp.dynamics, bvp.cost, x0, xf, t), t_max/100, t_max)
    t !== nothing ? t : golden_section(cost, t_max/100, t_max)
end

## Compiled Steering Functions (enabled by `import SymPy`; compiled `SteeringBVP`s return `BVPControl`s)
function __init__()
    @require SymPy="24249f21-da20-56a4-8eb1-6a02cf4ae2e6" include("sympy_bvp_compilation.jl")
end

end # module


text:
abstract type QueryStruct end
abstract type QWidget end
abstract type DataAttribute <: QueryStruct end

############ Providers

struct DataAttributeContext
    ctx::Set{Tuple{String,String,String}}
end

struct Sysstat <: DataAttribute
    property::String
    #context::DataAttributeContext
end
#Sysstat(name::String) = Sysstat(name, DataAttributeContext(Set()))

struct Aws <: DataAttribute
    property::String
end

struct Postfix <: DataAttribute
    property::String
end

##############

using Dates

mutable struct DataStore
    data::Dict{Tuple{String,String},Dict{String,Tuple{DataAttributeContext,Array{DateTime,1},Array{Float64,1}}}}
    points::Dict{Tuple{String,String},Array{DateTime,1}}
    descriptions::Dict{Tuple{String,String},String}
end
DataStore() = DataStore(Dict(), Dict(), Dict())

############### widgets

mutable struct QPlot <: QWidget
    title::String
    days::Int64
    data::Array{<:QueryStruct,1}
    size::Union{Tuple{Int64,Int64},Nothing}
end
# line plot of the given structure
QPlot(
    title::String;
    days::Int64 = 3,
    data::Array{<:QueryStruct,1} = [],
    size::Union{Tuple{Int64,Int64},Nothing} = nothing,
) = QPlot(title, days, data, size)

mutable struct QStack <: QueryStruct
    stacked::Array{<:QueryStruct,1}
    ctxs::Union{Array{DataAttributeContext, 1},Nothing}
end
# stacks the select lines
QStack(data::Array{<:QueryStruct,1}; ctxs = nothing) = QStack(data, ctxs)
QStack(content::QueryStruct; ctxs = nothing) = QStack([content], ctxs)

# overrides the default styles for the given data
#   fillrange: nothing to prevent filling (for example when stacked) or a float to fill towards this number
#   color: line color
#   fillcolor: fill color
#   fillalpha: fill alpha
mutable struct QStyle <: QueryStruct
    data::QueryStruct
    overloads::Dict{String,Any}
end

mutable struct QGroup <: QueryStruct
    unit::String
    max::Union{Nothing,<:Number}
    min::Union{Nothing,<:Number}
    log::Bool
    data::Array{<:QueryStruct,1}
    ctxs::Union{Array{DataAttributeContext, 1},Nothing}
end
# assigns the unit to this group. Use this, to plot different units in the same plot
QGroup(
    unit::String;
    max::Union{Nothing,<:Number} = nothing,
    min::Union{Nothing,<:Number} = nothing,
    log::Bool = false,
    data::Array{<:QueryStruct,1} = [],
    ctxs = nothing
) = QGroup(unit, max, min, log, data, ctxs)

###########

mutable struct Palette
    palette
    index
end

mutable struct PlotContext
    plot
    palette::Palette
    styles::Dict{String,Any}
    today::DateTime
    store::DataStore
    maxDate::DateTime
    minDate::DateTime
    group::Int
end
PlotContext(;
    plot = nothing,
    palette = Palette(nothing, 1),
    styles::Dict{String,Any} = defaultStyles(),
    today::DateTime = Dates.now(),
    store::DataStore = DataStore(),
    maxDate::DateTime = Dates.now(),
    minDate::DateTime = Dates.now(),
    group::Int = 0,
) = PlotContext(plot, palette, styles, today, store, maxDate, minDate, group)

max_stars_repo_name: Agapanthus/server-diary

max_issues_repo_name: Agapanthus/server-diary

max_forks_repo_name: Agapanthus/server-diary


text:
using Test
import Pluto
import Pluto: update_run!, WorkspaceManager, ClientSession, ServerSession, Notebook, Cell
import Distributed

@testset "Bonds" begin

    🍭 = ServerSession()
    🍭.options.evaluation.workspace_use_distributed = false
    fakeclient = ClientSession(:fake, nothing)
    🍭.connected_clients[fakeclient.id] = fakeclient

    @testset "AbstractPlutoDingetjes.jl" begin
        🍭.options.evaluation.workspace_use_distributed = true
        notebook = Notebook([
                # 1
                Cell("""
                begin
                    import AbstractPlutoDingetjes
                    const APD = AbstractPlutoDingetjes
                    import AbstractPlutoDingetjes.Bonds
                end
                """),
                # 2
                Cell("""
                begin
                    struct HTMLShower
                        f::Function
                    end
                    Base.show(io::IO, m::MIME"text/html", hs::HTMLShower) = hs.f(io)
                end
                """),
                Cell("""
                APD.is_inside_pluto()
                """),
                # 4
                Cell("""
                HTMLShower() do io
                    write(io, string(APD.is_inside_pluto(io)))
                end
                """),
                Cell("""
                HTMLShower() do io
                    write(io, string(APD.is_supported_by_display(io, Bonds.initial_value)))
                end
                """),
                Cell("""
                HTMLShower() do io
                    write(io, string(APD.is_supported_by_display(io, Bonds.transform_value)))
                end
                """),
                Cell("""
                HTMLShower() do io
                    write(io, string(APD.is_supported_by_display(io, Bonds.validate_value)))
                end
                """),
                Cell("""
                HTMLShower() do io
                    write(io, string(APD.is_supported_by_display(io, sqrt)))
                end
                """),
                # 9
                Cell("""
                @bind x_simple html"<input type=range>"
                """),
                Cell("""
                x_simple
                """),
                # 11
                Cell("""
                begin
                    struct OldSlider end

                    Base.show(io::IO, m::MIME"text/html", os::OldSlider) = write(io, "<input type=range value=1>")

                    Base.get(os::OldSlider) = 1
                end
                """),
                Cell("""
                @bind x_old OldSlider()
                """),
                Cell("""
                x_old
                """),
                # 14
                Cell("""
                begin
                    struct NewSlider end

                    Base.show(io::IO, m::MIME"text/html", os::NewSlider) = write(io, "<input type=range value=1>")

                    Bonds.initial_value(os::NewSlider) = 1

                    Bonds.possible_values(s::NewSlider) = [1,2,3]
                end
                """),
                Cell("""
                @bind x_new NewSlider()
                """),
                Cell("""
                x_new
                """),
                # 17
                Cell("""
                begin
                    struct TransformSlider end

                    Base.show(io::IO, m::MIME"text/html", os::TransformSlider) = write(io, "<input type=range value=1>")

                    Bonds.initial_value(os::TransformSlider) = "x"
                    Bonds.possible_values(os::TransformSlider) = 1:10
                    Bonds.transform_value(os::TransformSlider, from_js) = repeat("x", from_js)
                end
                """),
                Cell("""
                @bind x_transform TransformSlider()
                """),
                Cell("""
                x_transform
                """),
                # 20
                Cell("""
                begin
                    struct BadTransformSlider end

                    Base.show(io::IO, m::MIME"text/html", os::BadTransformSlider) = write(io, "<input type=range value=1>")

                    Bonds.initial_value(os::BadTransformSlider) = "x"
                    Bonds.possible_values(os::BadTransformSlider) = 1:10
                    Bonds.transform_value(os::BadTransformSlider, from_js) = error("bad")
                end
                """),
                Cell("""
                @bind x_badtransform BadTransformSlider()
                """),
                Cell("""
                x_badtransform
                """),
                # 23
                Cell("""
                count = Ref(0)
                """),
                Cell("""
                @bind x_counter NewSlider() #or OldSlider(), same idea
                """),
                Cell("""
                let
                    x_counter
                    count[] += 1
                end
                """),
                # 26
                Cell("""
                @assert x_old == 1
                """),
                Cell("""
                @assert x_new == 1
                """),
                Cell("""
                @assert x_transform == "x"
                """),
                # 29
                Cell("""
                begin
                    struct PossibleValuesTest
                        possible_values::Any
                    end

                    Base.show(io::IO, m::MIME"text/html", ::PossibleValuesTest) = write(io, "hello")

                    Bonds.possible_values(pvt::PossibleValuesTest) = pvt.possible_values
                end
                """),
                Cell("@bind pv1 PossibleValuesTest(Bonds.NotGiven())"),
                Cell("@bind pv2 PossibleValuesTest(Bonds.InfinitePossibilities())"),
                Cell("@bind pv3 PossibleValuesTest([1,2,3])"),
                Cell("@bind pv4 PossibleValuesTest((x+1 for x in 1:10))"),
                # 34
                Cell("@bind pv5 PossibleValuesTest(1:10)"),
            ])
        fakeclient.connected_notebook = notebook

        function set_bond_value(name, value, is_first_value=false)
            notebook.bonds[name] = Dict("value" => value)
            Pluto.set_bond_values_reactive(; session=🍭, notebook, bound_sym_names=[name],
                is_first_values=[is_first_value],
                run_async=false,
            )
        end

        # before loading AbstractPlutoDingetjes, test the default behaviour:
        update_run!(🍭, notebook, notebook.cells[9:10])
        @test noerror(notebook.cells[9])
        @test noerror(notebook.cells[10])

        @test Pluto.possible_bond_values(🍭, notebook, :x_simple) == :NotGiven
        @test notebook.cells[10].output.body == "missing"
        set_bond_value(:x_simple, 1, true)
        @test notebook.cells[10].output.body == "1"

        update_run!(🍭, notebook, notebook.cells)

        @test noerror(notebook.cells[1])
        @test noerror(notebook.cells[2])
        @test noerror(notebook.cells[3])
        @test noerror(notebook.cells[4])
        @test noerror(notebook.cells[5])
        @test noerror(notebook.cells[6])
        @test noerror(notebook.cells[7])
        @test noerror(notebook.cells[8])

        @test notebook.cells[3].output.body == "true"
        @test notebook.cells[4].output.body == "true"
        @test notebook.cells[5].output.body == "true"
        @test notebook.cells[6].output.body == "true"
        @test notebook.cells[7].output.body == "false"
        @test notebook.cells[8].output.body == "false"

        @test noerror(notebook.cells[9])
        @test noerror(notebook.cells[10])
        @test noerror(notebook.cells[11])
        @test noerror(notebook.cells[12])
        @test noerror(notebook.cells[13])
        @test noerror(notebook.cells[14])
        @test noerror(notebook.cells[15])
        @test noerror(notebook.cells[16])
        @test noerror(notebook.cells[17])
        @test noerror(notebook.cells[18])
        @test noerror(notebook.cells[19])
        @test noerror(notebook.cells[20])
        @test noerror(notebook.cells[21])
        @test noerror(notebook.cells[22])
        @test noerror(notebook.cells[23])
        @test noerror(notebook.cells[24])
        @test noerror(notebook.cells[25])
        @test noerror(notebook.cells[26])
        @test noerror(notebook.cells[27])
        @test noerror(notebook.cells[28])
        @test noerror(notebook.cells[29])
        @test noerror(notebook.cells[30])
        @test noerror(notebook.cells[31])
        @test noerror(notebook.cells[32])
        @test noerror(notebook.cells[33])
        @test noerror(notebook.cells[34])
        @test length(notebook.cells) == 34

        @test Pluto.possible_bond_values(🍭, notebook, :x_new) == [1,2,3]
        @test_throws Exception Pluto.possible_bond_values(🍭, notebook, :asdfasdfx_new)
        @test Pluto.possible_bond_values(🍭, notebook, :pv1) == :NotGiven
        @test Pluto.possible_bond_values(🍭, notebook, :pv2) == :InfinitePossibilities
        @test Pluto.possible_bond_values(🍭, notebook, :pv3) == [1,2,3]
        @test Pluto.possible_bond_values(🍭, notebook, :pv4) == 2:11
        @test Pluto.possible_bond_values(🍭, notebook, :pv5) === 1:10
        @test Pluto.possible_bond_values_length(🍭, notebook, :pv1) == :NotGiven
        @test Pluto.possible_bond_values_length(🍭, notebook, :pv2) == :InfinitePossibilities
        @test Pluto.possible_bond_values_length(🍭, notebook, :pv3) == 3
        @test Pluto.possible_bond_values_length(🍭, notebook, :pv4) == 10
        @test Pluto.possible_bond_values_length(🍭, notebook, :pv5) == 10

        @test notebook.cells[10].output.body == "missing"
        set_bond_value(:x_simple, 1, true)
        @test notebook.cells[10].output.body == "1"

        @test notebook.cells[13].output.body == "1"
        set_bond_value(:x_old, 1, true)
        @test notebook.cells[13].output.body == "1"
        set_bond_value(:x_old, 99, false)
        @test notebook.cells[13].output.body == "99"

        @test notebook.cells[16].output.body == "1"
        set_bond_value(:x_new, 1, true)
        @test notebook.cells[16].output.body == "1"
        set_bond_value(:x_new, 99, false)
        @test notebook.cells[16].output.body == "99"

        @test notebook.cells[19].output.body == "\"x\""
        set_bond_value(:x_transform, 1, true)
        @test notebook.cells[19].output.body == "\"x\""
        set_bond_value(:x_transform, 3, false)
        @test notebook.cells[19].output.body == "\"xxx\""

        @test notebook.cells[22].output.body != "missing"
        @info "The following error is expected:"
        set_bond_value(:x_badtransform, 1, true)
        @test notebook.cells[22].output.body != "missing"

        @test notebook.cells[25].output.body == "1"
        set_bond_value(:x_counter, 1, true)
        @test notebook.cells[25].output.body == "1"
        set_bond_value(:x_counter, 7, false)
        @test notebook.cells[25].output.body == "2"

        WorkspaceManager.unmake_workspace((🍭, notebook))
        🍭.options.evaluation.workspace_use_distributed = false

        # test that the notebook file is runnable:

        test_proc = Distributed.addprocs(1)[1]

        Distributed.remotecall_eval(Main, test_proc, quote
            import Pkg
            try
                Pkg.UPDATED_REGISTRY_THIS_SESSION[] = true
            catch; end
            Pkg.activate(mktempdir())
            Pkg.add("AbstractPlutoDingetjes")
        end)
        @test Distributed.remotecall_eval(Main, test_proc, quote
            include($(notebook.path))
            true
        end)
        Distributed.rmprocs(test_proc)
    end

    @testset "Dependent Bound Variables" begin
        🍭 = ServerSession()
        🍭.options.evaluation.workspace_use_distributed = false
        fakeclient = ClientSession(:fake, nothing)
        🍭.connected_clients[fakeclient.id] = fakeclient
        🍭.options.evaluation.workspace_use_distributed = true
        notebook = Notebook([
            Cell(raw"""@bind x HTML("<input type=range min=1 max=10>")"""),
            Cell(raw"""@bind y HTML("<input type=range min=1 max=$(x)>")"""),
            Cell(raw"""x"""), #3
            Cell(raw"""y"""), #4
            Cell(raw"""
            begin
                struct TransformSlider
                    range::AbstractRange
                end

                Base.show(io::IO, m::MIME"text/html", os::TransformSlider) = write(io, "<input type=range value=$(minimum(os.range)) min=$(minimum(os.range)) max=$(maximum(os.range))>")

                Bonds.initial_value(os::TransformSlider) = Bonds.transform_value(os, minimum(os.range))
                Bonds.possible_values(os::TransformSlider) = os.range
                Bonds.transform_value(os::TransformSlider, from_js) = from_js * 2
            end
            """),
            Cell(raw"""begin
                hello1 = 123
                @bind a TransformSlider(1:10)
            end"""),
            Cell(raw"""begin
                hello2 = 234
                @bind b TransformSlider(1:a)
            end"""),
            Cell(raw"""a"""), #8
            Cell(raw"""b"""), #9
            Cell(raw"""hello1"""), #10
            Cell(raw"""hello2"""), #11
            Cell(raw"""using AbstractPlutoDingetjes"""),
        ])
        fakeclient.connected_notebook = notebook
        update_run!(🍭, notebook, notebook.cells)

        function set_bond_values!(notebook:: Notebook, bonds:: Dict; is_first_value=false)
            for (name, value) in bonds
                notebook.bonds[name] = Dict("value" => value)
            end
            Pluto.set_bond_values_reactive(; session=🍭, notebook, bound_sym_names=collect(keys(bonds)), run_async=false, is_first_values=fill(is_first_value, length(bonds)))
        end

        @test notebook.cells[3].output.body == "missing"
        @test notebook.cells[4].output.body == "missing" # no initial value defined for simple html slider (in contrast to TransformSlider)
        @test notebook.cells[8].output.body == "2"
        @test notebook.cells[9].output.body == "2"
        @test notebook.cells[10].output.body == "123"
        @test notebook.cells[11].output.body == "234"

        set_bond_values!(notebook, Dict(:x => 1, :a => 1); is_first_value=true)
        @test notebook.cells[3].output.body == "1"
        @test notebook.cells[4].output.body == "missing" # no initial value defined for simple html slider (in contrast to TransformSlider)
        @test notebook.cells[8].output.body == "2" # TransformSlider scales values *2
        @test notebook.cells[9].output.body == "2"
        @test notebook.cells[10].output.body == "123"
        @test notebook.cells[11].output.body == "234"

        set_bond_values!(notebook, Dict(:y => 1, :b => 1); is_first_value=true)
        @test notebook.cells[3].output.body == "1"
        @test notebook.cells[4].output.body == "1"
        @test notebook.cells[8].output.body == "2"
        @test notebook.cells[9].output.body == "2"
        @test notebook.cells[10].output.body == "123"
        @test notebook.cells[11].output.body == "234"

        set_bond_values!(notebook, Dict(:x => 5))
        @test notebook.cells[3].output.body == "5"
        @test notebook.cells[4].output.body == "missing" # the slider object is re-defined, therefore its value is the default one

        set_bond_values!(notebook, Dict(:y => 3))
        @test notebook.cells[3].output.body == "5"
        @test notebook.cells[4].output.body == "3"

        set_bond_values!(notebook, Dict(:x => 10, :y => 5))
        @test notebook.cells[3].output.body == "10"
        @test notebook.cells[4].output.body == "5" # this would fail without PR #2014 - previously `y` was reset to the default value `missing`

        set_bond_values!(notebook, Dict(:b => 2))
        @test notebook.cells[8].output.body == "2"
        @test notebook.cells[9].output.body == "4"
        @test notebook.cells[10].output.body == "123"
        @test notebook.cells[11].output.body == "234"

        set_bond_values!(notebook, Dict(:a => 8, :b => 12))
        @test notebook.cells[8].output.body == "16"
        @test notebook.cells[9].output.body == "24" # this would fail without PR #2014
        @test notebook.cells[10].output.body == "123"
        @test notebook.cells[11].output.body == "234"

        set_bond_values!(notebook, Dict(:a => 1, :b => 1))
        setcode!(notebook.cells[10], "a + hello1")
        setcode!(notebook.cells[11], "b + hello2")
        update_run!(🍭, notebook, notebook.cells[10:11])

        @test notebook.cells[10].output.body == "125"
        @test notebook.cells[11].output.body == "236"

        set_bond_values!(notebook, Dict(:a => 2, :b => 2))
        @test notebook.cells[10].output.body == "127"
        @test notebook.cells[11].output.body == "238"
        set_bond_values!(notebook, Dict(:b => 3))
        @test notebook.cells[10].output.body == "127"
        @test notebook.cells[11].output.body == "240"
        set_bond_values!(notebook, Dict(:a => 1))
        @test notebook.cells[10].output.body == "125"
        @test notebook.cells[11].output.body == "236" # changing a will reset b

        WorkspaceManager.unmake_workspace((🍭, notebook))

    end
end


text:
module TestIteration
    using Test, DataFrames

    dv = [1, 2, missing]
    dm = Union{Int, Missing}[1 2; 3 4]
    df = Array{Union{Int, Missing}}(zeros(2, 2, 2))

    df = DataFrame(A = Vector{Union{Int, Missing}}(1:2), B = Vector{Union{Int, Missing}}(2:3))

    @test size(eachrow(df)) == (size(df, 1),)
    @test eachrow(df)[1] == DataFrameRow(df, 1)
    @test collect(eachrow(df)) isa Vector{DataFrameRow{DataFrame}}
    @test eltype(eachrow(df)) == DataFrameRow{DataFrame}
    for row in eachrow(df)
        @test isa(row, DataFrameRow)
        @test (row[:B] - row[:A]) == 1
        # issue #683 (https://github.com/JuliaData/DataFrames.jl/pull/683)
        @test collect(row) isa Vector{Pair{Symbol, Int}}
    end

    @test size(eachcol(df)) == (size(df, 2),)
    @test length(eachcol(df)) == size(df, 2)
    @test eachcol(df)[1] == df[:, 1]
    @test collect(eachcol(df)) isa Vector{Tuple{Symbol, Any}}
    @test eltype(eachcol(df)) == Tuple{Symbol, Any}
    for col in eachcol(df)
        @test isa(col, Tuple{Symbol, AbstractVector})
    end

    @test map(x -> minimum(convert(Array, x)), eachrow(df)) == Any[1,2]
    @test map(minimum, eachcol(df)) == DataFrame(A = [1], B = [2])

    row = DataFrameRow(df, 1)

    row[:A] = 100
    @test df[1, :A] == 100

    row[1] = 101
    @test df[1, :A] == 101

    df = DataFrame(A = Vector{Union{Int, Missing}}(1:4), B = Union{String, Missing}["M", "F", "F", "M"])

    s1 = view(df, 1:3)
    s1[2,:A] = 4
    @test df[2, :A] == 4
    @test view(s1, 1:2) == view(df, 1:2)

    s2 = view(df, 1:2:3)
    s2[2, :B] = "M"
    @test df[3, :B] == "M"
    @test view(s2, 1:1:2) == view(df, [1,3])

    # @test_fail for x in df; end # Raises an error
end


text:
## Load the packages.
using PyPlot
using HDF5
using Statistics
using Interpolations

## Read the data.
x_01 = h5read("drycbl.01.00002400.h5", "x")
x_02 = h5read("drycbl.02.00002400.h5", "x")
y_01 = h5read("drycbl.01.00002400.h5", "y")
y_02 = h5read("drycbl.02.00002400.h5", "y")
z_01 = h5read("drycbl.01.00002400.h5", "z")
z_02 = h5read("drycbl.02.00002400.h5", "z")
zh_01 = h5read("drycbl.01.00002400.h5", "zh")
zh_02 = h5read("drycbl.02.00002400.h5", "zh")

s_bot_01 = h5read("drycbl.01.00002400.h5", "s_bot", (:, :))
s_bot_02 = h5read("drycbl.02.00002400.h5", "s_bot", (:, :))

s_bot_01 .-= mean(s_bot_01)
s_bot_02 .-= mean(s_bot_02)

println(zh_01[13], " ", zh_02[37])
w_xy_01 = h5read("drycbl.01.00002400.h5", "w", (:, :, 13)) # at 300 m
w_xy_02 = h5read("drycbl.02.00002400.h5", "w", (:, :, 37))

w_xz_01 = h5read("drycbl.01.00002400.h5", "w", (:, 1, :))
w_xz_02 = h5read("drycbl.02.00002400.h5", "w", (:, 1, :))

## Read for plotting.
xh_01 = h5read("drycbl.01.00002400.h5", "xh")
xh_02 = h5read("drycbl.02.00002400.h5", "xh")
yh_01 = h5read("drycbl.01.00002400.h5", "yh")
yh_02 = h5read("drycbl.02.00002400.h5", "yh")
push!(xh_01, 3200.)
push!(xh_02, 3200.)
push!(yh_01, 3200.)
push!(yh_02, 3200.)

## Surface temperature figure.
figure(figsize=(9,4.5))
subplot(121)
pcolormesh(xh_01, yh_01, s_bot_01', vmin=-0.6, vmax=0.6, cmap=plt.cm.cividis)
xlabel("x (m)")
ylabel("y (m)")

subplot(122)
pcolormesh(xh_02, yh_02, s_bot_02', vmin=-0.6, vmax=0.6, cmap=plt.cm.cividis)
xlabel("x (m)")
ylabel("y (m)")

tight_layout()
display(gcf())

println("variance in s_bot: ", var(s_bot_01), " ", var(s_bot_02))

## Surface temperature figure tweaked
s_01 = copy(s_bot_01)
s_02 = copy(s_bot_02)
@. s_01[s_01 < 0.1] = NaN
@. s_02[s_02 < 0.1] = NaN

figure(figsize=(5,4.5))
pcolormesh(xh_02, yh_02, s_02', vmin=0., vmax=0.6, cmap=PyPlot.cm.Reds, alpha=0.4)
pcolormesh(xh_01, yh_01, s_01', vmin=0., vmax=0.6, cmap=PyPlot.cm.Blues, alpha=0.4)
xlabel("x (m)")
ylabel("y (m)")

tight_layout()
display(gcf())

## Vertical velocity xy
figure(figsize=(9,4.5))
subplot(121)
pcolormesh(xh_01, yh_01, w_xy_01', vmin=-1.5, vmax=1.5, cmap=plt.cm.seismic)
xlabel("x (m)")
ylabel("y (m)")

subplot(122)
pcolormesh(xh_02, yh_02, w_xy_02', vmin=-1.5, vmax=1.5, cmap=plt.cm.seismic)
xlabel("x (m)")
ylabel("y (m)")

tight_layout()
display(gcf())

println("variance in w_xy: ", var(w_xy_01), " ", var(w_xy_02))

## Vertical velocity xy interpolated
interp = LinearInterpolation((x_02, y_02), w_xy_02, extrapolation_bc=Periodic())
w_xy_02_to_01 = interp(x_01, y_01)

"""
figure(figsize=(9,4.5))
subplot(121)
pcolormesh(xh_01, yh_01, w_xy_01', vmin=-1.5, vmax=1.5, cmap=plt.cm.seismic)
xlabel("x (m)")
ylabel("y (m)")

subplot(122)
pcolormesh(xh_01, yh_01, w_xy_02_to_01', vmin=-1.5, vmax=1.5, cmap=plt.cm.seismic)
xlabel("x (m)")
ylabel("y (m)")

tight_layout()
display(gcf())
"""

## Scatter xy
w_min = max(minimum(w_xy_01), minimum(w_xy_02_to_01))
w_max = min(maximum(w_xy_01), maximum(w_xy_02_to_01))

figure(figsize=(5.5, 4.5))
hexbin(vec(w_xy_01), vec(w_xy_02_to_01), gridsize=60, bins="log", cmap=plt.cm.cividis)
colorbar()
plot([-2, 5], [-2, 5], "w:")
xlim(w_min, w_max)
ylim(w_min, w_max)
xlabel(L"w_xy_01 (m s$^{-1}$)")
ylabel(L"w_xy_02 (m s$^{-1}$)")

tight_layout()
display(gcf())

println("correlation in w_xy: ", cor(vec(w_xy_01), vec(w_xy_02_to_01)))

## Vertical velocity xz
zp_01 = copy(z_01)
zp_02 = copy(z_02)
pushfirst!(zp_01, 0.)
pushfirst!(zp_02, 0.)
push!(zp_01, 3200.)
push!(zp_02, 3200.)

figure(figsize=(9,4.5))
subplot(121)
pcolormesh(xh_01, zp_01, w_xz_01', vmin=-1.5, vmax=1.5, cmap=plt.cm.seismic)
xlabel("x (m)")
ylabel("y (m)")

subplot(122)
pcolormesh(xh_02, zp_02, w_xz_02', vmin=-1.5, vmax=1.5, cmap=plt.cm.seismic)
xlabel("x (m)")
ylabel("y (m)")

tight_layout()
display(gcf())

show()


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    signarea(p₁, p₂, p₃)

Compute signed area of triangle formed
by points `p₁`, `p₂` and `p₃`.
"""
function signarea(p₁::Point{2}, p₂::Point{2}, p₃::Point{2})
  a = coordinates(p₁)
  b = coordinates(p₂)
  c = coordinates(p₃)
  ((b[1]-a[1])*(c[2]-a[2]) - (b[2]-a[2])*(c[1]-a[1])) / 2
end

"""
    signarea(triangle)

Compute the signed area of `triangle`.
"""
function signarea(t::Triangle{2})
  vs = vertices(t)
  signarea(vs[1], vs[2], vs[3])
end

"""
    sideof(point, segment)

Determines on which side of the oriented `segment`
the `point` lies. Possible results are `:LEFT`,
`:RIGHT` or `:ON` the segment.
"""
function sideof(p::Point{2,T}, s::Segment{2,T}) where {T}
  a, b = vertices(s)
  area = signarea(p, a, b)
  ifelse(area > atol(T), :LEFT, ifelse(area < -atol(T), :RIGHT, :ON))
end

"""
    sideof(point, chain)

Determines on which side of the closed `chain` the
`point` lies. Possible results are `:INSIDE` or
`:OUTSIDE` the chain.
"""
function sideof(p::Point{2,T}, c::Chain{2,T}) where {T}
  w = windingnumber(p, c)
  ifelse(isapprox(w, zero(T), atol=atol(T)), :OUTSIDE, :INSIDE)
end


text:
# Defines some function to find linear parameters
#
export solve_linear_parameters

@doc raw"""

Find optimal linear parameters
"""
function solve_linear_parameters(model::Abstract_Model2Fit,
                                 X::AbstractVector,
                                 Y::AbstractVector,
                                 θ::AbstractVector,
                                 linear_param_idx::AbstractVector{Int})
    @assert length(X) == length(Y)
    @assert length(θ) == parameter_size(model)

    # Build J of the least squares problem
    #
    # | J.θ_linear - Y |^2
    #
    θ_copy = copy(θ)
    θ_copy[linear_param_idx] .= 0

    J = zeros(promote_type(eltype(X),eltype(Y),eltype(θ)),
              length(X),
              length(linear_param_idx))

    for (col_idx,θ_idx) in enumerate(linear_param_idx)
        θ_copy[θ_idx] = 1
        accumulate_y!(model,(@view J[:,col_idx]),X,θ_copy)
        θ_copy[θ_idx] = 0
    end

    # Solve the problem
    #
    # H.θ_linear = rhs
    #
    # where H = J'J, rhs = J'Y
    #
    rhs = J'*Y
    H = J'*J

    θ_linear = H \ rhs

    θ_copy[linear_param_idx] .= θ_linear

    θ_copy
end


text:
type Account <: AbstractTimeSeries
    timestamp::Vector{Date}
    values::Matrix{Float64}
    colnames::Vector{ASCIIString}
    portfolios::Vector{Portfolio}

    function Account(timestamp::Vector{Date},
                          values::Matrix{Float64},
                          colnames::Vector{ASCIIString},
                          portfolios::Vector{Account})

                          nrow, ncol = size(values, 1), size(values, 2)
                          nrow != size(timestamp, 1) ? error("values must match length of timestamp"):
                          ncol != size(colnames,1) ? error("column names must match width of array"):
                          timestamp != unique(timestamp) ? error("there are duplicate dates"):
                          ~(flipud(timestamp) == sort(timestamp) || timestamp == sort(timestamp)) ? error("dates are mangled"):
                          flipud(timestamp) == sort(timestamp) ?
                          new(flipud(timestamp), flipud(values), colnames):
                          new(timestamp, values, colnames)
    end
end


text:
# [[file:~/GitHub/J4Org.jl/docs/main.org::*A%20documented%20Julia%20=Foo=%20module][A documented Julia =Foo= module:1]]
module Foo

export Point, norm, foo

import Base: norm

#+Point L:Point_struct
# This is my Point structure
#
# *Example:*
#
# Creates a point $p$ of coordinates $(x=1,y=2)$.
#
# #+BEGIN_SRC julia :eval never :exports code
# p=Point(1,2)
# #+END_SRC
#
# You can add any valid Org mode directive. If you want to use
# in-documentation link, use [[norm_link_example][]]
#
struct Point
    x::Float64
    y::Float64
end

#+Point
# Creates Point at origin $(0,0)$
Point() = Point(0,0)

#+Enum
# An enum
@enum Alpha A=1 B C # just for example

#+Point,Method L:norm_link_example
# A simple function that computes $\sqrt{x^2+y^2}$
#
# *Example:*
#!p=Point(1.0,2.0);
#!norm(p)
#
# See: [[Point_struct][]]
#
norm(p::Point)::Float64 = sqrt(p.x*p.x+p.y*p.y)

# +Method,Internal
# An internal function
#
# For symbol that are not exported, do not forget the "Foo." prefix:
# !p=Point(1.0,2.0)
# !Foo.foo(2.0,p)
foo(r::Float64,p::Point) = Point(r*p.x,r*p.y)

end
# A documented Julia =Foo= module:1 ends here


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "ROCmCompilerSupport"
version = v"4.0.0"

# Collection of sources required to build
sources = [
    ArchiveSource("https://github.com/RadeonOpenCompute/ROCm-CompilerSupport/archive/rocm-$(version).tar.gz",
                  "f389601fb70b2d9a60d0e2798919af9ddf7b8376a2e460141507fe50073dfb31"), # 4.0.0
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
ln -s ${bindir}/../tools/* ${bindir}/

cd ${WORKSPACE}/srcdir/ROCm-CompilerSupport*/lib/comgr

# first manually build bc2h utility
cd cmake && mkdir bc2h && cd bc2h
cp ../bc2h.cmake CMakeLists.txt
CC=$HOSTCC cmake -DCMAKE_C_COMPILER=$HOSTCC .
make
cd ../..

# then build everything else
atomic_patch -p1 $WORKSPACE/srcdir/patches/disable-1031.patch
atomic_patch -p1 $WORKSPACE/srcdir/patches/disable-tests.patch
atomic_patch -p1 $WORKSPACE/srcdir/patches/disable-bc2h.patch
atomic_patch -p1 $WORKSPACE/srcdir/patches/fix-objdump-clopt.patch
mkdir build && cd build
cp ../cmake/bc2h/bc2h .
# TODO: -DROCM_DIR=${prefix}
cmake -DCMAKE_PREFIX_PATH=${prefix} \
      -DCMAKE_INSTALL_PREFIX=${prefix} \
      -DCMAKE_BUILD_TYPE=Release \
      -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN%.*}_clang.cmake \
      -DLLVM_DIR="${prefix}/lib/cmake/llvm" \
      -DClang_DIR="${prefix}/lib/cmake/clang" \
      -DLLD_DIR="${prefix}/lib/cmake/ldd" \
      ..
PATH=../cmake/bc2h:$PATH make -j${nproc}
make install

find $WORKSPACE/destdir/bin -type l -delete
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    # TODO: cxx03
    Platform("x86_64", "linux"; libc="glibc", cxxstring_abi="cxx11"),
    Platform("x86_64", "linux"; libc="musl", cxxstring_abi="cxx11"),
]

# The products that we will ensure are always built
products = [
    LibraryProduct(["libamd_comgr"], :libamd_comgr),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("hsa_rocr_jll"),
    Dependency("ROCmDeviceLibs_jll"),
    BuildDependency(PackageSpec(; name="LLVM_full_jll", version=v"11.0.1")),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies,
               preferred_gcc_version=v"8", preferred_llvm_version=v"11")

max_stars_repo_name: jameskermode/Yggdrasil

max_issues_repo_name: jameskermode/Yggdrasil

max_forks_repo_name: jameskermode/Yggdrasil


text:
Base.iszero(t::Tuple) = iszero(t[1]) && iszero(t[2])
Base.occursin(re, char::Char) = occursin(re, string(char))

rev(char) = char == '(' ? ')' : '{'


text:
struct PolyakIterable{T, F, S, M}
    w0::T
    f::F
    stepsize::S
    momentum::M
end

PolyakIterable(w0, f; stepsize, momentum) = PolyakIterable(w0, f, to_iterator(stepsize), to_iterator(momentum))

Base.IteratorSize(::Type{<:PolyakIterable}) = Base.IsInfinite()

mutable struct PolyakState{T, R, G, S, M}
    w::T
    f_w::R
    grad_f_w::G
    stepsize_iterator::S
    momentum_iterator::M
    step::T
end

function Base.iterate(iter::PolyakIterable)
    w = copy(iter.w0)
    grad_f_w, f_w = gradient(iter.f, w)
    state = PolyakState(w, f_w, grad_f_w, Iterators.Stateful(iter.stepsize), Iterators.Stateful(iter.momentum), zero(w))
    return IterationOutput(state.w, state.f_w, state.grad_f_w), state
end

function Base.iterate(iter::PolyakIterable, state::PolyakState)
    stepsize = popfirst!(state.stepsize_iterator)
    momentum = popfirst!(state.momentum_iterator)
    state.step .= momentum .* state.step .- stepsize .* state.grad_f_w
    state.w .+= state.step
    state.grad_f_w, state.f_w = gradient(iter.f, state.w)
    return IterationOutput(state.w, state.f_w, state.grad_f_w), state
end

struct Polyak
    kwargs
    Polyak(; kwargs...) = new(kwargs)
end

(alg::Polyak)(args...) = PolyakIterable(args...; alg.kwargs...)


text:
using SnoopCompileBot

snoop_bench(BotConfig("SnoopCompileWithMocking"));


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

abstract AbstractSolver

type Solver{S<:AbstractSolver}
    name :: AbstractString       # some descriptive name for problem
    time :: Float64              # current time
    problems :: Vector{Problem}
    norms :: Vector{Tuple}       # solution norms for convergence studies
    ndofs :: Int                 # number of degrees of freedom in problem
    xdmf :: Nullable{Xdmf}       # input/output handle
    initialized :: Bool
    u :: Vector{Float64}
    la :: Vector{Float64}
    properties :: S
end

function Solver{S<:AbstractSolver}(::Type{S}, name="solver", properties...)
    variant = S(properties...)
    solver = Solver{S}(name, 0.0, [], [], 0, nothing, false, [], [], variant)
    return solver
end

function Solver{S<:AbstractSolver}(::Type{S}, problems::Problem...)
    solver = Solver(S, "$(S)Solver")
    push!(solver.problems, problems...)
    return solver
end

function get_problems(solver::Solver)
    return solver.problems
end

function push!(solver::Solver, problem)
    push!(solver.problems, problem)
end

function getindex(solver::Solver, problem_name)
    for problem in get_problems(solver)
        if problem.name == problem_name
            return problem
        end
    end
    throw(KeyError(problem_name))
end

# one-liner helpers to identify problem types

is_field_problem(problem) = false
is_field_problem{P<:FieldProblem}(problem::Problem{P}) = true
is_boundary_problem(problem) = false
is_boundary_problem{P<:BoundaryProblem}(problem::Problem{P}) = true
get_field_problems(solver::Solver) = filter(is_field_problem, get_problems(solver))
get_boundary_problems(solver::Solver) = filter(is_boundary_problem, get_problems(solver))

"""
Posthook for field assembly. By default, do nothing.
This can be used to make some modifications for assembly
after all elements are assembled.

Examples
--------
function field_assembly_posthook!(solver::Solver,
                                  K::SparseMatrixCSC,
                                  Kg::SparseMatrixCSC,
                                  f::SparseMatrixCSC,
                                  fg::SpareMatrixCSC)
    info("doing stuff, size(K) = ", size(K))
end
"""
function field_assembly_posthook!
end

"""Return one combined field assembly for a set of field problems.

Parameters
----------
solver :: Solver

Returns
-------
M, K, Kg, f, fg :: SparseMatrixCSC

Notes
-----
If several field problems exists, they are simply summed together, so
problems must have unique node ids.

"""
function get_field_assembly(solver::Solver; show_info=true)
    problems = get_field_problems(solver)

    M = SparseMatrixCOO()
    K = SparseMatrixCOO()
    Kg = SparseMatrixCOO()
    f = SparseMatrixCOO()
    fg = SparseMatrixCOO()

    for problem in problems
        append!(M, problem.assembly.M)
        append!(K, problem.assembly.K)
        append!(Kg, problem.assembly.Kg)
        append!(f, problem.assembly.f)
        append!(fg, problem.assembly.fg)
    end

    if solver.ndofs == 0
        solver.ndofs = size(K, 1)
        show_info && info("automatically determined problem dimension, ndofs = $(solver.ndofs)")
    end

    M = sparse(M, solver.ndofs, solver.ndofs)
    K = sparse(K, solver.ndofs, solver.ndofs)
    if nnz(K) == 0
        warn("Field assembly seems to be empty. Check that elements are pushed to problem and formulation is correct.")
    end
    Kg = sparse(Kg, solver.ndofs, solver.ndofs)
    f = sparse(f, solver.ndofs, 1)
    fg = sparse(fg, solver.ndofs, 1)

    # run any posthook for assembly if defined
    args = Tuple{Solver, SparseMatrixCSC, SparseMatrixCSC, SparseMatrixCSC, SparseMatrixCSC}
    if method_exists(field_assembly_posthook!, args)
        field_assembly_posthook!(solver, K, Kg, fg, fg)
    end

    return M, K, Kg, f, fg
end

""" Posthook for boundary assembly. By default, do nothing. """
function boundary_assembly_posthook!
end

""" Return one combined boundary assembly for a set of boundary problems.

Returns
-------
C1, C2, D, g :: SparseMatrixCSC

Notes
-----
When some dof is constrained by multiple boundary problems an algorithm is
launched what tries to do it's best to solve issue. It's far from perfect
but is able to handle some basic situations occurring in corner nodes and
crosspoints.

"""
function get_boundary_assembly(solver::Solver)
    ndofs = solver.ndofs
    @assert ndofs != 0
    K = spzeros(ndofs, ndofs)
    C1 = spzeros(ndofs, ndofs)
    C2 = spzeros(ndofs, ndofs)
    D = spzeros(ndofs, ndofs)
    f = spzeros(ndofs, 1)
    g = spzeros(ndofs, 1)
    for problem in get_boundary_problems(solver)
        assembly = problem.assembly
        K_ = sparse(assembly.K, ndofs, ndofs)
        C1_ = sparse(assembly.C1, ndofs, ndofs)
        C2_ = sparse(assembly.C2, ndofs, ndofs)
        D_ = sparse(assembly.D, ndofs, ndofs)
        f_ = sparse(assembly.f, ndofs, 1)
        g_ = sparse(assembly.g, ndofs, 1)
        # check for overconstraint situation and handle it if possible
        already_constrained = get_nonzero_rows(C2)
        new_constraints = get_nonzero_rows(C2_)
        overconstrained_dofs = intersect(already_constrained, new_constraints)
        if length(overconstrained_dofs) != 0
            overconstrained_dofs = sort(overconstrained_dofs)
            overconstrained_nodes = find_nodes_by_dofs(problem, overconstrained_dofs)
            handle_overconstraint_error!(problem, overconstrained_nodes,
                overconstrained_dofs, C1, C1_, C2, C2_, D, D_, g, g_)
        end
        K += K_
        C1 += C1_
        C2 += C2_
        D += D_
        f += f_
        g += g_
    end
    return K, C1, C2, D, f, g
end

function resize!(A::SparseMatrixCSC, m::Int64, n::Int64)
    (n == A.n) && (m == A.m) && return
    @assert n >= A.n
    @assert m >= A.m
    append!(A.colptr, A.colptr[end]*ones(Int, m-A.m))
    A.n = n
    A.m = m
end

"""
Given C and g, construct new basis such that v = P*u + g

Parameters
----------
S  set of linearly independent dofs.
"""
function create_projection(C::SparseMatrixCSC, g; S=nothing, tol=1.0e-12)
    n, m = size(C)
    @assert n == m
    if S == nothing
        S = get_nonzero_rows(C)
    end
    # FIXME: this creates dense matrices
    # efficiency / memory usage is a question
    M = get_nonzero_columns(C)
    F = qrfact(C[S,:])
    P = spzeros(n,m)
    P[:,M] = sparse(F \ full(C[S,M]))
    h = sparse(F \ full(g[S]))
    resize!(P, n, m)
    resize!(h, n, 1)
    P = speye(n) - P
    SparseMatrix.droptol!(P, tol)
    return P, h
end

""" Assume C is invertible. """
function create_projection(C, g, ::Type{Val{:invertible}})
    nz1, nz2 = get_nonzeros(C)
    P = spzeros(size(C)...)
    for j=1:size(C,1)
        j in nz1 && continue
        P[j,j] = 1.0
    end
    v = lufact(C[nz1,nz2]) \ full(g[nz1])
    return P, v
end

"""
Solve linear system using LDLt factorization (SuiteSparse). This version
requires that final system is symmetric and positive definite, so boundary
conditions are first eliminated before solution.
"""
function solve!(solver::Solver, K, C1, C2, D, f, g, u, la, ::Type{Val{1}}; debug=false)

    nnz(D) == 0 || return false
    nz = get_nonzero_rows(C2)
    B = get_nonzero_rows(C2')
    # C2^-1 exists or this doesn't work
    length(nz) == length(B) || return false

    A = get_nonzero_rows(K)
    I = setdiff(A, B)

    if debug
        info("# nz = $(length(nz))")
        info("# A = $(length(A))")
        info("# B = $(length(B))")
        info("# I = $(length(I))")
    end

    # solver boundary dofs
    try
        u[B] = lufact(C2[nz,B]) \ full(g[nz])
    catch
        info("solver #1 failed to solve boundary dofs (you should not see this message).")
        info("# nz = $(length(nz))")
        info("# A = $(length(A))")
        info("# B = $(length(B))")
        info("# I = $(length(I))")
        info("nz = $nz")
        info("B = $B")
        rethrow()
    end

    # solve interior domain using LDLt factorization
    F = ldltfact(K[I,I])
    u[I] = F \ (f[I] - K[I,B]*u[B])
    # solve lambda
    la[B] = lufact(C1[B,nz]) \ full(f[B] - K[B,I]*u[I] - K[B,B]*u[B])

    return true
end

"""
Solve linear system using LU factorization (UMFPACK). This version solves
directly the saddle point problem without elimination of boundary conditions.
"""
function solve!(solver::Solver, K, C1, C2, D, f, g, u, la, ::Type{Val{2}})
    # construct global system Ax = b and solve using lufact (UMFPACK)
    A = [K C1'; C2  D]
    b = [f; g]
    x = lufact(A) \ full(b)
    u[:] = x[1:solver.ndofs]
    la[:] = x[solver.ndofs+1:end]
    return true
end

""" Default linear system solver for solver. """
function solve!(solver::Solver; empty_assemblies_before_solution=true, show_info=true)
    show_info && info("Solving problems ...")
    t0 = Base.time()

    # assemble field & boundary problems
    # TODO: return same kind of set for both assembly types
    # M1, K1, Kg1, f1, fg1, C11, C21, D1, g1 = get_field_assembly(solver)
    # M2, K2, Kg2, f2, fg2, C12, C22, D2, g2 = get_boundary_assembly(solver)

    M, K, Kg, f, fg = get_field_assembly(solver)
    Kb, C1, C2, D, fb, g = get_boundary_assembly(solver)
    K = K + Kg + Kb
    f = f + fg + fb
    K = 1/2*(K + K')
    M = 1/2*(M + M')

    nz = ones(solver.ndofs)
    nz[get_nonzero_rows(C2)] = 0.0
    nz[get_nonzero_rows(D)] = 0.0
    D += spdiagm(nz)

    # free up some memory before solution
    for problem in get_problems(solver)
        if empty_assemblies_before_solution
            empty!(problem.assembly)
        else
            optimize!(problem.assembly)
        end
        gc()
    end

    ndofs = solver.ndofs
    u = zeros(ndofs)
    la = zeros(ndofs)
    status = false
    i = 0
    for i in [1, 2]
        status = solve!(solver, K, C1, C2, D, f, g, u, la, Val{i})
        status && break
    end
    status || error("Failed to solve linear system!")
    t1 = round(Base.time()-t0, 2)
    norms = (norm(u), norm(la))
    show_info && info("Solved problems in $t1 seconds using solver $i. Solution norms = $norms.")
    push!(solver.norms, norms)

    solver.u = u
    solver.la = la

    return
end

""" Default assembler for solver. """
function assemble!(solver::Solver; show_info=true, timing=true)
    show_info && info("Assembling problems ...")
    t0 = Base.time()
    assembly_times = Dict()
    nproblems = 0
    ndofs = 0
    for problem in solver.problems
        t00 = Base.time()
        empty!(problem.assembly)
        assemble!(problem, solver.time)
        nproblems += 1
        Ks = size(problem.assembly.K, 2)
        Cs = size(problem.assembly.C1, 2)
        ndofs = max(ndofs, Ks, Cs)
        t11 = Base.time()
        assembly_times[problem.name] = t11-t00
    end
    solver.ndofs = ndofs
    t1 = round(Base.time()-t0, 2)
    show_info && info("Assembled $nproblems problems in $t1 seconds. ndofs = $ndofs.")
    if timing
        info("Assembly times:")
        for (i, problem) in enumerate(solver.problems)
            pn = problem.name
            pt = round(assembly_times[pn], 2)
            info("$i $pn $pt")
        end
    end
end

function get_unknown_fields(solver::Solver)
    fields = Dict()
    for problem in get_field_problems(solver)
        field_name = get_unknown_field_name(problem)
        field_dim = get_unknown_field_dimension(problem)
        fields[field_name] = field_dim
    end
    return fields
end

function get_unknown_field_name(solver::Solver)
    fields = get_unknown_fields(solver)
    return join(sort(collect(keys(fields))), ", ")
end

function get_unknown_field_dimension(solver::Solver)
    fields = get_unknown_fields(solver)
    return sum(values(fields))
end

""" Default initializer for solver. """
function initialize!(solver::Solver; show_info=true)
    if solver.initialized
        show_info && info("initialize!(): solver already initialized")
        return
    end
    show_info && info("Initializing solver ...")
    problems = get_problems(solver)
    length(problems) != 0 || error("Empty solver, add problems to solver using push!")
    t0 = Base.time()
    field_problems = get_field_problems(solver)
    length(field_problems) != 0 || warn("No field problem found from solver, add some..?")
    field_name = get_unknown_field_name(solver)
    field_dim = get_unknown_field_dimension(solver)
    info("initialize!(): looks we are solving $field_name, $field_dim dofs/node")
    nodes = Set{Int64}()
    for problem in problems
        initialize!(problem, solver.time)
        for element in get_elements(problem)
            conn = get_connectivity(element)
            push!(nodes, conn...)
        end
    end
    nnodes = length(nodes)
    info("Total number of nodes in problems: $nnodes")
    maxdof = maximum(nnodes)*field_dim
    info("# of max dof (=size of solution vector) is $maxdof")
    solver.u = zeros(maxdof)
    solver.la = zeros(maxdof)
    # TODO: this could be used to initialize elements too...
    # TODO: cannot initialize to zero always, construct vector from elements.
    for problem in problems
        problem.assembly.u = zeros(maxdof)
        problem.assembly.la = zeros(maxdof)
        # initialize(problem, ....)
    end
    t1 = round(Base.time()-t0, 2)
    show_info && info("Initialized solver in $t1 seconds.")
    solver.initialized = true
end

function get_all_elements(solver::Solver)
    elements = [get_elements(problem) for problem in get_problems(solver)]
    return [elements...;]
end

function get_element_type{E}(element::Element{E})
    return E
end

function get_element_id{E}(element::Element{E})
    return element.id
end

function is_element_type{E}(element::Element{E}, element_type)
    return is(E, element_type)
end

function filter_by_element_type(element_type, elements)
    return filter(element -> is_element_type(element, element_type), elements)
end

function call(solver::Solver, field_name::AbstractString, time::Float64)
    fields = [problem(field_name, time) for problem in get_problems(solver)]
    return merge(fields...)
end

function get_temporal_collection(xdmf::Xdmf)
    domain = find_element(xdmf.xml, "Domain")
    grid = nothing
    if domain == nothing
        info("Xdmf: creating new temporal collection")
        domain = new_child(xdmf.xml, "Domain")
        grid = new_child(domain, "Grid")
        set_attribute(grid, "CollectionType", "Temporal")
        set_attribute(grid, "GridType", "Collection")
    end
    grid = find_element(domain, "Grid")
    return grid
end

""" Default update for solver. """
function update!{S}(solver::Solver{S}; show_info=true)
    u = solver.u
    la = solver.la

    show_info && info("Updating problems ...")
    t0 = Base.time()

    for problem in solver.problems
        assembly = get_assembly(problem)
        elements = get_elements(problem)
        # update solution, first for assembly (u,la) ...
        update!(problem, assembly, u, la)
        # .. and then from assembly (u,la) to elements
        update!(problem, assembly, elements, solver.time)
    end

    # if io is attached to solver, update hdf / xml also
    if !isnull(solver.xdmf)
        update_xdmf!(solver)
    end

    t1 = round(Base.time()-t0, 2)
    show_info && info("Updated problems in $t1 seconds.")
end

function update_xdmf!{S}(solver::Solver{S}; show_info=true)
    xdmf = get(solver.xdmf)
    temporal_collection = get_temporal_collection(xdmf)

    frame = new_element("Grid")
    new_child(frame, "Time", Dict("Value" => solver.time))

    # save geometry
    X = solver("geometry", solver.time)
    node_ids = sort(collect(keys(X)))
    geometry = hcat([X[nid] for nid in node_ids]...)
    ndim, nnodes = size(geometry)
    geom_type = ndim == 2 ? "XY" : "XYZ"
    dataitem = new_dataitem(xdmf, "/Node IDs", node_ids)
    geom = new_child(frame, "Geometry", Dict("Type" => geom_type))
    dataitem = new_dataitem(xdmf, "/Geometry", geometry)
    add_child(geom, dataitem)

    # save topology
    all_elements = get_all_elements(solver)
    nelements = length(all_elements)
    element_types = unique(map(get_element_type, all_elements))

    xdmf_element_mapping = Dict(
        "Seg2" => "Polyline",
        "Tri3" => "Triangle",
        "Quad4" => "Quadrilateral",
        "Tet4" => "Tetrahedron",
        "Pyramid5" => "Pyramid",
        "Wedge6" => "Wedge",
        "Hex8" => "Hexahedron",
        "Seg3" => "Edge_3",
        "Tri6" => "Tri_6",
        "Quad8" => "Quad_8",
        "Tet10" => "Tet_10",
        "Pyramid13" => "Pyramid_13",
        "Wedge15" => "Wedge_15",
        "Hex20" => "Hex_20")

    for element_type in element_types
        elements = filter_by_element_type(element_type, all_elements)
        sort!(elements, by=get_element_id)
        element_ids = map(get_element_id, elements)
        element_conn = map(get_connectivity, elements)
        element_conn = transpose(hcat(element_conn...)) - 1
        element_code = split(string(element_type), ".")[end]
        dataitem = new_dataitem(xdmf, "/Topology/$element_code/Element IDs", element_ids)
        dataitem = new_dataitem(xdmf, "/Topology/$element_code/Connectivity", element_conn)
        topology = new_child(frame, "Topology")
        set_attribute(topology, "TopologyType", xdmf_element_mapping[element_code])
        set_attribute(topology, "NumberOfElements", length(elements))
        add_child(topology, dataitem)
    end

    # save solved fields
    unknown_field_name = get_unknown_field_name(solver)
    U = solver(unknown_field_name, solver.time)
    node_ids2 = sort(collect(keys(U)))

    @assert node_ids == node_ids2
    ndim = length(U[first(node_ids)])
    field_type = ndim == 1 ? "Scalar" : "Vector"
    field_center = "Node"
    if ndim == 2
        for nid in node_ids
            U[nid] = [U[nid]; 0.0]
        end
        ndim = 3
    end
    U = hcat([U[nid] for nid in node_ids]...)
    unknown_field_name = ucfirst(unknown_field_name)
    time = solver.time
    path = ""
    if S == Nonlinear
        iteration = solver.properties.iteration
        path = "/Results/Time $time/Iteration $iteration/Nodal Fields/$unknown_field_name"
    elseif S == Linear
        path = "/Results/Time $time/Nodal Fields/$unknown_field_name"
    end
    dataitem = new_dataitem(xdmf, path, U)
    attribute = new_child(frame, "Attribute")
    set_attribute(attribute, "Name", unknown_field_name)
    set_attribute(attribute, "Center", field_center)
    set_attribute(attribute, "AttributeType", field_type)
    add_child(attribute, dataitem)
    if (S == Linear) || ((S == Nonlinear) && has_converged(solver))
        add_child(temporal_collection, frame)
    end
    save!(xdmf)
end

### Nonlinear quasistatic solver

type Nonlinear <: AbstractSolver
    iteration :: Int        # iteration counter
    min_iterations :: Int64 # minimum number of iterations
    max_iterations :: Int64 # maximum number of iterations
    convergence_tolerance :: Float64
    error_if_no_convergence :: Bool # throw error if no convergence
end

function Nonlinear()
    solver = Nonlinear(0, 1, 10, 5.0e-5, true)
    return solver
end

""" Check convergence of problems.

Notes
-----
Default convergence criteria is obtained by checking each sub-problem convergence.
"""
function has_converged(solver::Solver{Nonlinear}; show_info=false,
                       check_convergence_for_boundary_problems=false)
    properties = solver.properties
    converged = true
    eps = properties.convergence_tolerance
    for problem in solver.problems
        has_converged = true
        if is_field_problem(problem)
            has_converged = problem.assembly.u_norm_change < eps
            if isapprox(norm(problem.assembly.u), 0.0)
                # trivial solution
                has_converged = true
            end
            show_info && info("Details for problem $(problem.name)")
            show_info && info("Norm: $(norm(problem.assembly.u))")
            show_info && info("Norm change: $(problem.assembly.u_norm_change)")
            show_info && info("Has converged? $(has_converged)")
        end
        if is_boundary_problem(problem) && check_convergence_for_boundary_problems
            has_converged = problem.assembly.la_norm_change/norm(problem.assembly.la) < eps
            show_info && info("Details for problem $(problem.name)")
            show_info && info("Norm: $(norm(problem.assembly.la))")
            show_info && info("Norm change: $(problem.assembly.la_norm_change)")
            show_info && info("Has converged? $(has_converged)")
        end
        converged &= has_converged
    end
    return converged
end

type NonlinearConvergenceError <: Exception
    solver :: Solver
end

function Base.showerror(io::IO, exception::NonlinearConvergenceError)
    max_iters = exception.solver.properties.max_iterations
    print(io, "nonlinear iteration did not converge in $max_iters iterations!")
end

""" Default solver for quasistatic nonlinear problems. """
function call(solver::Solver{Nonlinear}; show_info=true)

    properties = solver.properties

    # 1. initialize each problem so that we can start nonlinear iterations
    initialize!(solver)

    # 2. start non-linear iterations
    for properties.iteration=1:properties.max_iterations
        show_info && info(repeat("-", 80))
        show_info && info("Starting nonlinear iteration #$(properties.iteration)")
        show_info && info("Increment time t=$(round(solver.time, 3))")
        show_info && info(repeat("-", 80))

        # 2.1 update linearized assemblies
        assemble!(solver)
        # 2.2 call solver for linearized system
        solve!(solver)
        # 2.3 update solution back to elements
        update!(solver)

        # 2.4 check convergence
        if has_converged(solver)
            info("Converged in $(properties.iteration) iterations.")
            properties.iteration >= properties.min_iterations && return true
            info("Convergence criteria met, but iteration < min_iterations, continuing...")
        end
    end

    # 3. did not converge
    properties.error_if_no_convergence && throw(NonlinearConvergenceError(solver))
end

""" Convenience function to call nonlinear solver. """
function NonlinearSolver(problems...)
    solver = Solver(Nonlinear, "default nonlinear solver")
    if length(problems) != 0
        push!(solver, problems...)
    end
    return solver
end
function NonlinearSolver(name::AbstractString, problems::Problem...)
    solver = NonlinearSolver(problems...)
    solver.name = name
    return solver
end

### Linear quasistatic solver

""" Quasistatic solver for linear problems.

Notes
-----
Main differences in this solver, compared to nonlinear solver are:
1. system of problems is assumed to converge in one step
2. reassembly of problem is done only if it's manually requested using empty!(problem.assembly)

"""
type Linear <: AbstractSolver
end

function assemble!(solver::Solver{Linear}; show_info=true)
    show_info && info("Assembling problems ...")
    tic()
    nproblems = 0
    ndofs = 0
    for problem in get_problems(solver)
        if isempty(problem.assembly)
            assemble!(problem, solver.time)
            nproblems += 1
        else
            show_info && info("$(problem.name) already assembled, skipping.")
        end
        ndofs = max(ndofs, size(problem.assembly.K, 2))
    end
    solver.ndofs = ndofs
    t1 = round(toq(), 2)
    show_info && info("Assembled $nproblems problems in $t1 seconds. ndofs = $ndofs.")
end

function call(solver::Solver{Linear}; show_info=true)
    t0 = Base.time()
    show_info && info(repeat("-", 80))
    show_info && info("Starting linear solver")
    show_info && info("Increment time t=$(round(solver.time, 3))")
    show_info && info(repeat("-", 80))
    initialize!(solver)
    assemble!(solver)
    solve!(solver)
    update!(solver)
    t1 = round(Base.time()-t0, 2)
    show_info && info("Linear solver ready in $t1 seconds.")
end

""" Convenience function to call linear solver. """
function LinearSolver(problems::Problem...)
    solver = Solver(Linear, "default linear solver")
    if length(problems) != 0
        push!(solver, problems...)
    end
    return solver
end
function LinearSolver(name::AbstractString, problems::Problem...)
    solver = LinearSolver(problems...)
    solver.name = name
    return solver
end

### End of linear quasistatic solver

### Postprocessor

type Postprocessor <: AbstractSolver
    assembly :: Assembly
    F :: Union{Factorization, Void}
end

function Postprocessor()
    Postprocessor(Assembly(), nothing)
end

function assemble!(solver::Solver{Postprocessor}; show_info=true)
    show_info && info("Assembling problems ...")
    tic()
    nproblems = 0
    ndofs = 0
    assembly = solver.properties.assembly
    empty!(assembly)
    for problem in get_problems(solver)
        for element in get_elements(problem)
            postprocess!(assembly, problem, element, solver.time)
        end
        nproblems += 1
        ndofs = max(ndofs, size(problem.assembly.K, 2))
    end
    solver.ndofs = ndofs
    t1 = round(toq(), 2)
    show_info && info("Assembled $nproblems problems in $t1 seconds. ndofs = $ndofs.")
end

function call(solver::Solver{Postprocessor}; show_info=true)
    t0 = Base.time()
    show_info && info(repeat("-", 80))
    show_info && info("Starting postprocessor")
    show_info && info("Increment time t=$(round(solver.time, 3))")
    show_info && info(repeat("-", 80))
    initialize!(solver)
    assemble!(solver)
    assembly = solver.properties.assembly
    M = sparse(assembly.M)
    f = sparse(assembly.f)
    F = cholfact(M)
    q = F \ f
    t1 = round(Base.time()-t0, 2)
    show_info && info("Postprocess of results ready in $t1 seconds.")
    return q
end

""" Convenience function to call postprocessor. """
function Postprocessor(problems::Problem...)
    solver = Solver(Postprocessor, "default postprocessor")
    if length(problems) != 0
        push!(solver, problems...)
    end
    return solver
end

function Postprocessor(name::AbstractString, problems::Problem...)
    solver = Postprocessor(problems...)
    solver.name = name
    return solver
end


text:
using Documenter, MonteCarloObservable

makedocs(
    # options

)

makedocs(
    modules = [MonteCarloObservable],
    format = :html,
    sitename = "MonteCarloObservable.jl",
    pages = [
        "Home" => "index.md",
        "Manual" => [
            "Getting started" => "manual/gettingstarted.md",
            "Type of the mean" => "manual/meantype.md",
            "Error estimation" => "manual/errorestimation.md",
            "Memory / disk storage" => "manual/memdisk.md"
        ],
        "Methods" => [
            # "All" => "methods/all.md",
            "General" => "methods/general.md",
            "Statistics" => "methods/statistics.md",
            "IO" => "methods/io.md",
            "Plotting" => "methods/plotting.md"
        ]
    ]
)

deploydocs(
    repo   = "github.com/crstnbr/MonteCarloObservable.jl.git",
    target = "build",
    deps   = nothing,
    make   = nothing,
    julia  = "release",
    osname = "linux"
)


text:
function solve()
    N, M, T = [parse(Int, x) for x in split(readline())]
    t = 0
    n = N
    for i=1:M
        A, B = [parse(Int, x) for x in split(readline())]
        n -= (A - t)
        if n ≤ 0
            return "No"
        end
        n += (B - A)
        n = min(n, N)
        t = B
    end
    n -= (T - t)
    n ≤ 0 ? "No" : "Yes"
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
module Lu_julia_dense
#
	using Lu
#
		ArgName = "julia_dense"
		Descr = "Julia's default algorithm on a dense matrix."
#
	type Dat{T}
		prob::Lu.Problem{T}
		A::Array{T,2}
	#
		Dat() = new()
	end
#
	function construct_dat(T::DataType)
		return Dat{T}()
	end
	function fill_dat{T}(prob::Lu.Problem, dat::Dat{T})
		dat.prob = prob
		dat.A = full(prob.A)
	end
	function solve_dat{T}(dat::Dat{T}, sol::Lu.Solution{T})
		sol.L, sol.U, sol.p = lu(dat.A)
		sol.solved = true
	end
end
#
#
module Lu_julia_sparse
#
	using Lu
#
	ArgName = "julia_sparse"
	Descr = "Julia's default algorithm on a sparse matrix (uses umfpack)."
#
	type Dat{T}
		prob::Lu.Problem{T}
		A::SparseMatrixCSC{T}
	#
		Dat() = new()
	end
#
	function construct_dat(T::DataType)
		return Dat{T}()
	end
	function fill_dat{T}(prob::Lu.Problem, dat::Dat{T})
		dat.prob = prob
		dat.A = sparse(prob.A)
	end
	function solve_dat{T}(dat::Dat{T}, sol::Lu.Solution{T})
		# Note, lufact with sparse matrices seems to only be possible with a specific type (e.g Float64)
		# try: methods(lufact)
		# https://github.com/JuliaLang/julia/issues/4439 says: "I do believe that UMFPACK does not have a Float32/Complex64 interface"
		luF = lufact(dat.A); L,U,p,q,Rs = luF[:(:)]
		sol.L = L; sol.U = U; sol.p = p; sol.q = q; sol.rs = Rs; sol.solved = true;
	end
end

max_stars_repo_name: jadnohra/mlcp-lib

max_issues_repo_name: jadnohra/mlcp-lib

max_forks_repo_name: jadnohra/mlcp-lib


text:
"""
    fields(tbl::AbstractTable)::Vector{Symbol}

Produce a list of column names of `tbl` in canonical order.
"""
function fields end

"""
    eltypes(tbl::AbstractTable[, fields...])::Vector{DataType}

Produce a list of eltypes of columns of `tbl` in canonical order. If `fields`
are specified, lists only eltypes of respective columns.
"""
function eltypes end

function eltypes(tbl::AbstractTable, fields...)
    idx = index(tbl)
    col_indices = [ idx[field] for field in fields ]
    return eltypes(tbl)[col_indices]
end

"""
    nrow(tbl::AbstractTable)::Int

Return the number of rows in `tbl` if known, -1 if unknown.
"""
function nrow end

"""
    index(tbl::AbstractTable)::Dict{Symbol, Int}

Produce a mapping from fields to numerical indices, where the latter respect
the canonical order of columns of `tbl`.
"""
function index end