text:
# Workaround for libz loading confusion.
@static if Sys.islinux()
    using ImageMagick
end

using GtkReactive, Gtk.ShortNames, IntervalSets, Graphics, Colors,
      TestImages, FileIO, FixedPointNumbers, RoundingIntegers, Dates, Cairo,
      IdentityRanges
using Test

rtask = Reactive.runner_task # starting with Reactive 0.7.0, this became a Ref
if isa(rtask, Base.RefValue)
    rtask = rtask[]
end
if !istaskdone(rtask)
    Reactive.stop()
    fetch(rtask)
end

include("tools.jl")

@testset "Widgets" begin
    ## label
    l = label("Hello")
    @test signal(l) == l.signal
    @test signal(signal(l)) == l.signal
    @test get_gtk_property(l, :label, String) == "Hello"
    push!(signal(l), "world")
    rr()
    @test get_gtk_property(l, :label, String) == "world"
    @test string(l) == string("Gtk.GtkLabelLeaf with ", string(signal(l)))
    # map with keywords
    lsig0 = map(l) do lbl  # "regular" map runs the function
        lbl
    end
    rr()
    @test value(lsig0) == "world"
    lsig = map(l; init="foo") do lbl   # with "init", you avoid running
        lbl
    end
    rr()
    @test value(lsig) == "foo"
    push!(l, "bar")
    rr()
    @test value(lsig) == "bar"

    ## checkbox
    w = Window("Checkbox")
    check = checkbox(label="click me")
    push!(w, check)
    Gtk.showall(w)
    @test value(check) == false
    @test Gtk.G_.active(check.widget) == false
    push!(check, true)
    rr()
    @test value(check)
    @test Gtk.G_.active(check.widget)
    destroy(w)

    ## togglebutton
    w = Window("Togglebutton")
    tgl = togglebutton(label="click me")
    push!(w, tgl)
    Gtk.showall(w)
    @test value(tgl) == false
    @test Gtk.G_.active(tgl.widget) == false
    push!(tgl, true)
    rr()
    @test value(tgl)
    @test Gtk.G_.active(tgl.widget)
    destroy(w)

    ## textbox (aka Entry)
    txt = textbox("Type something")
    num = textbox(5, range=1:10)
    win = Window("Textboxes") |> (bx = Box(:h))
    push!(bx, txt)
    push!(bx, num)
    Gtk.showall(win)
    @test get_gtk_property(txt, :text, String) == "Type something"
    push!(txt, "ok")
    rr()
    @test get_gtk_property(txt, :text, String) == "ok"
    set_gtk_property!(txt, :text, "other direction")
    signal_emit(widget(txt), :activate, Nothing)
    rr()
    @test value(txt) == "other direction"
    @test get_gtk_property(num, :text, String) == "5"
    push!(signal(num), 11, (sig, val, osig, capex) -> throw(capex.ex))
    @test_throws ArgumentError rr()
    push!(num, 8)
    rr()
    @test get_gtk_property(num, :text, String) == "8"
    meld = map(txt, num) do t, n
        join((t, n), 'X')
    end
    rr()
    @test value(meld) == "other directionX8"
    push!(num, 4)
    rr()
    @test value(meld) == "other directionX4"
    push!(txt, "4")
    rr()
    @test value(meld) == "4X4"
    destroy(win)

    ## textarea (aka TextView)
    v = textarea("Type something longer")
    win = Window(v)
    Gtk.showall(win)
    @test value(v) == "Type something longer"
    push!(v, "ok")
    rr()
    @test get_gtk_property(Gtk.G_.buffer(v.widget), :text, String) == "ok"
    destroy(win)

    ## slider
    s = slider(1:15)
    sleep(0.01)    # For the Gtk eventloop
    @test value(s) == 8
    push!(s, 3)
    rr()
    @test value(s) == 3
    s3 = slider(IdentityRange(-3:3))
    sleep(0.01)
    @test value(s3) == 0
    push!(s3, -3)
    rr()
    @test value(s3) == -3

    # Use a single signal for two widgets
    s2 = slider(1:15, signal=signal(s), orientation='v')
    @test value(s2) == 3
    push!(s2, 11)
    rr()
    @test value(s) == 11
    destroy(s2)
    destroy(s)

    # Updating the limits of the slider
    s = slider(1:15)
    sleep(0.01)    # For the Gtk eventloop
    @test value(s) == 8
    push!(s, 1:7, 5)
    sleep(0.01)
    rr()
    @test value(s) == 5

    ## dropdown
    dd = dropdown(("Strawberry", "Vanilla", "Chocolate"))
    @test value(dd) == "Strawberry"
    push!(dd, "Chocolate")
    rr()
    @test get_gtk_property(dd, :active, Int) == 2
    destroy(dd.widget)

    r = Ref(0)
    dd = dropdown(["Five"=>x->x[]=5,
                   "Seven"=>x->x[]=7])
    ddsig = map(f->f(r), dd.mappedsignal)
    rr()
    @test value(dd) == "Five"
    @test r[] == 5
    push!(dd, "Seven")
    run_till_empty()
    @test value(dd) == "Seven"
    @test r[] == 7
    push!(dd, "Five")
    rr()
    @test r[] == 5
    destroy(dd.widget)

    ## spinbutton
    s = spinbutton(1:15)
    sleep(0.01)    # For the Gtk eventloop
    @test value(s) == 1
    push!(s, 3)
    rr()
    @test value(s) == 3
    destroy(s)

    s = spinbutton(0:59, orientation="vertical")
    @test G_.orientation(Orientable(widget(s))) == Gtk.GConstants.GtkOrientation.VERTICAL
    destroy(s)

    # Updating the limits of the spinbutton
    s = spinbutton(1:15)
    sleep(0.01)    # For the Gtk eventloop
    @test value(s) == 1
    push!(s, 1:7, 5)
    rr()
    @test value(s) == 5

    ## cyclicspinbutton
    a = spinbutton(1:10, value = 5)
    carry_up = Signal(false)
    foreach(carry_up; init=nothing) do up
        push!(a, value(a) - (-1)^up)
    end
    b = cyclicspinbutton(1:3, carry_up)
    run_till_empty()
    @test value(a) == 5
    @test value(b) == 1
    push!(b, 2)
    run_till_empty()
    @test value(a) == 5
    @test value(b) == 2
    push!(b, 0)
    run_till_empty()
    @test value(a) == 4
    @test value(b) == 3
    push!(b, 4)
    run_till_empty()
    @test value(a) == 5
    @test value(b) == 1
    destroy(a)

    s = cyclicspinbutton(0:59, carry_up, orientation="vertical")
    @test G_.orientation(Orientable(widget(s))) == Gtk.GConstants.GtkOrientation.VERTICAL
    destroy(s)

    # timewidget
    t = Dates.Time(1,1,1)
    s = Signal(t)
    tw = timewidget(t, signal=s)
    run_till_empty()
    @test value(tw) == value(s) == t
    t = Dates.Time(2,2,2)
    push!(tw, t)
    run_till_empty()
    @test value(tw) == value(s) == t
    t = Dates.Time(3,3,3)
    push!(s, t)
    run_till_empty()
    @test value(tw) == value(s) == t

    # datetimewidget
    t = DateTime(1,1,1,1,1,1)
    s = Signal(t)
    tw = datetimewidget(t, signal=s)
    run_till_empty()
    @test value(tw) == value(s) == t
    t = DateTime(2,2,2,2,2,2)
    push!(tw, t)
    run_till_empty()
    @test value(tw) == value(s) == t
    t = DateTime(3,3,3,3,3,3)
    push!(s, t)
    run_till_empty()
    @test value(tw) == value(s) == t

    # progressbar
    pb = progressbar(1..10)
    @test value(pb) == 1
    push!(pb, 5)
    run_till_empty()
    @test value(pb) == 5
    pb = progressbar(2:8)
    @test value(pb) == 2

end

## button
# For reasons I don't understand, this often fails if it's inside a @testset
counter = 0

w = Window("Widgets")
b = button("Click me")
push!(w, b)
action = map(b) do val
    global counter
    counter::Int += 1
end
Gtk.showall(w)
rr()
cc = counter  # map seems to fire it once, so record the "new" initial value
click(b::GtkReactive.Button) = ccall((:gtk_button_clicked,Gtk.libgtk),Cvoid,(Ptr{Gtk.GObject},),b.widget)
GC.gc(true)
click(b)
rr()
@test counter == cc+1
destroy(w)

# Make sure we can also put a ToolButton in a Button
button(; widget=ToolButton("Save as..."))

if Gtk.libgtk_version >= v"3.10"
    # To support GtkBuilder, we need this as the minimum libgtk version
    @testset "Compound widgets" begin
        ## player widget
        s = Signal(1)
        p = player(s, 1:8)
        win = Window() |> (g = Grid())
        g[1,1] = p
        Gtk.showall(win)
        rr()
        btn_fwd = p.widget.step_forward
        @test value(s) == 1
        push!(btn_fwd, nothing)
        run_till_empty()
        @test value(s) == 2
        push!(p.widget.play_forward, nothing)
        for i = 1:7
            run_till_empty()
            sleep(0.1)
        end
        @test value(s) == 8
        destroy(win)

        p = player(1:1000)
        win = Window(frame(p))
        Gtk.showall(win)
        push!(widget(p).direction, 1)
        destroy(win)  # this should not generate a lot of output
    end
end

@testset "CairoUnits" begin
    x = UserUnit(0.2)
    @test UserUnit(x) === x
    @test convert(UserUnit, x) === x
    @test x+x === UserUnit(0.2+0.2)
    @test x-x === UserUnit(0.0)
    @test Float64(x) === 0.2
    @test convert(Float64, x) === 0.2
    y = UserUnit(-0.3)
    @test x > y
    @test y < x
    @test abs(x) === x
    @test abs(y) === UserUnit(0.3)
    @test min(x, y) === y
    @test max(x, y) === x
    z = DeviceUnit(2.0)
    @test_throws ErrorException x+z
    @test Bool(DeviceUnit(1.0)) === true
    @test Integer(DeviceUnit(3.0)) === 3
end

@testset "Canvas" begin
    @test XY(5, 5) === XY{Int}(5, 5)
    @test XY(5, 5.0) === XY{Float64}(5.0, 5.0)
    @test XY{UserUnit}(5, 5.0) === XY{UserUnit}(5.0, 5.0) === XY{UserUnit}(UserUnit(5), UserUnit(5))
    @test XY(5.0, 5)+XY(4, 4.1) === XY(9, 9.1)
    @test XY(5, 5)-XY(4, 4) === XY(1, 1)

    @test isa(MouseButton{UserUnit}(), MouseButton{UserUnit})
    @test isa(MouseButton{DeviceUnit}(), MouseButton{DeviceUnit})
    @test isa(MouseScroll{UserUnit}(), MouseScroll{UserUnit})
    @test isa(MouseScroll{DeviceUnit}(), MouseScroll{DeviceUnit})

    @test BoundingBox(XY(2..4, -15..15)) === BoundingBox(2, 4, -15, 15)

    c = canvas(208, 207)
    win = Window(c)
    Gtk.showall(win)
    sleep(0.1)
    @test Graphics.width(c) == 208
    @test Graphics.height(c) == 207
    @test isa(c, GtkReactive.Canvas{DeviceUnit})
    destroy(win)
    c = canvas(UserUnit, 208, 207)
    win = Window(c)
    Gtk.showall(win)
    reveal(c, true)
    sleep(0.3)
    @test isa(c, GtkReactive.Canvas{UserUnit})
    corner_dev = (DeviceUnit(208), DeviceUnit(207))
    for (coords, corner_usr) in ((BoundingBox(0, 1, 0, 1), (UserUnit(1), UserUnit(1))),
                                 (ZoomRegion((5:10, 3:5)), (UserUnit(5), UserUnit(10))),
                                 ((-1:1, 101:110), (UserUnit(110), UserUnit(1))))
        set_coordinates(c, coords)
        @test GtkReactive.convertunits(UserUnit, c, corner_dev...) == corner_usr
        @test GtkReactive.convertunits(DeviceUnit, c, corner_dev...) == corner_dev
        @test GtkReactive.convertunits(UserUnit, c, corner_usr...) == corner_usr
        @test GtkReactive.convertunits(DeviceUnit, c, corner_usr...) == corner_dev
    end

    destroy(win)

    c = canvas()
    f = Frame(c)
    @test isa(f, Gtk.GtkFrameLeaf)
    destroy(f)
    c = canvas()
    f = AspectFrame(c, "Some title", 0.5, 0.5, 3.0)
    @test isa(f, Gtk.GtkAspectFrameLeaf)
    @test get_gtk_property(f, :ratio, Float64) == 3.0
    destroy(f)
end

# @testset "Canvas events" begin
    win = Window() |> (c = canvas(UserUnit))
    Gtk.showall(win)
    sleep(0.2)
    lastevent = Ref("nothing")
    press   = map(btn->lastevent[] = "press",   c.mouse.buttonpress)
    release = map(btn->lastevent[] = "release", c.mouse.buttonrelease)
    motion  = map(btn->lastevent[] = string("motion to ", btn.position.x, ", ", btn.position.y),
                  c.mouse.motion)
    scroll  = map(btn->lastevent[] = "scroll", c.mouse.scroll)
    rr()
    lastevent[] = "nothing"
    @test lastevent[] == "nothing"
    signal_emit(widget(c), "button-press-event", Bool, eventbutton(c, BUTTON_PRESS, 1))
    sleep(0.1)
    rr()
    @test lastevent[] == "press"
    signal_emit(widget(c), "button-release-event", Bool, eventbutton(c, GtkReactive.BUTTON_RELEASE, 1))
    sleep(0.1)
    rr()
    sleep(0.1)
    rr()
    @test lastevent[] == "release"
    signal_emit(widget(c), "scroll-event", Bool, eventscroll(c, UP))
    sleep(0.1)
    rr()
    sleep(0.1)
    rr()
    @test lastevent[] == "scroll"
    signal_emit(widget(c), "motion-notify-event", Bool, eventmotion(c, 0, UserUnit(20), UserUnit(15)))
    sleep(0.1)
    rr()
    sleep(0.1)
    rr()
    @test lastevent[] == "motion to UserUnit(20.0), UserUnit(15.0)"
    destroy(win)
# end

@testset "Popup" begin
    popupmenu = Menu()
    popupitem = MenuItem("Popup menu...")
    push!(popupmenu, popupitem)
    Gtk.showall(popupmenu)
    win = Window() |> (c = canvas())
    popuptriggered = Ref(false)
    push!(c.preserved, map(c.mouse.buttonpress) do btn
        if btn.button == 3 && btn.clicktype == BUTTON_PRESS
            popup(popupmenu, btn.gtkevent)  # use the raw Gtk event
            popuptriggered[] = true
            nothing
        end
    end)
    yield()
    @test !popuptriggered[]
    evt = eventbutton(c, BUTTON_PRESS, 1)
    signal_emit(widget(c), "button-press-event", Bool, evt)
    yield()
    @test !popuptriggered[]
    evt = eventbutton(c, BUTTON_PRESS, 3)
    signal_emit(widget(c), "button-press-event", Bool, evt)
    run_till_empty()
    @test popuptriggered[]
    destroy(win)
    destroy(popupmenu)
end

@testset "Drawing" begin
    img = testimage("lighthouse")
    c = canvas(UserUnit, size(img, 2), size(img, 1))
    win = Window(c)
    xsig, ysig = Signal(20), Signal(20)
    draw(c, xsig, ysig) do cnvs, x, y
        copy!(c, img)
        ctx = getgc(cnvs)
        set_source(ctx, colorant"red")
        set_line_width(ctx, 2)
        circle(ctx, x, y, 5)
        stroke(ctx)
    end
    Gtk.showall(win)
    rr()
    push!(xsig, 100)
    rr()
    sleep(1)
    # Check that we get the right answer
    fn = tempname()
    Cairo.write_to_png(getgc(c).surface, fn)
    imgout = load(fn)
    rm(fn)
    @test imgout[25,100] == imgout[16,100] == imgout[20,105] == colorant"red"
    @test imgout[20,100] == img[20,100]
    destroy(win)
end

# For testing ZoomRegion support for non-AbstractArray objects
struct Foo end
Base.axes(::Foo) = (Base.OneTo(7), Base.OneTo(9))

@testset "Zoom/pan" begin
    @test string(UserUnit(3)) == "UserUnit(3.0)"
    @test string(DeviceUnit(3)) == "DeviceUnit(3.0)"

    xy = @inferred(XY(1, 3))
    @test isa(xy, XY{Int})
    @test xy.x == 1
    @test xy.y == 3
    @test string(xy) == "XY(1, 3)"
    xy = @inferred(XY{Float64}(1, 3))
    @test isa(xy, XY{Float64})
    @test xy.x == 1
    @test xy.y == 3
    @test string(xy) == "XY(1.0, 3.0)"
    @test isa(convert(XY{Int}, xy), XY{Int})
    xy = XY{Float64}(3.2, 4.8)
    xyr = convert(XY{RInt}, xy)
    @test isa(xyr, XY{RInt}) && xyr.x == 3 && xyr.y == 5
    xy = XY(UserUnit(3), UserUnit(5))
    @test string(xy) == "XY{UserUnit}(3.0, 5.0)"
    @test @inferred(XY{UserUnit}(3, 5)) == xy

    zr = ZoomRegion((1:80, 1:100))  # y, x order
    zrz = GtkReactive.zoom(zr, 0.5)
    @test zrz.currentview.x == 26..75
    @test zrz.currentview.y == 21..60
    zrp = GtkReactive.pan_x(zrz, 0.2)
    @test zrp.currentview.x == 36..85
    @test zrp.currentview.y == 21..60
    zrp = GtkReactive.pan_x(zrz, -0.2)
    @test zrp.currentview.x == 16..65
    @test zrp.currentview.y == 21..60
    zrp = GtkReactive.pan_y(zrz, -0.2)
    @test zrp.currentview.x == 26..75
    @test zrp.currentview.y == 13..52
    zrp = GtkReactive.pan_y(zrz, 0.2)
    @test zrp.currentview.x == 26..75
    @test zrp.currentview.y == 29..68
    zrp = GtkReactive.pan_x(zrz, 1.0)
    @test zrp.currentview.x == 51..100
    @test zrp.currentview.y == 21..60
    zrp = GtkReactive.pan_y(zrz, -1.0)
    @test zrp.currentview.x == 26..75
    @test zrp.currentview.y == 1..40
    zrz2 = GtkReactive.zoom(zrz, 2.0001)
    @test zrz2 == zr
    zrz2 = GtkReactive.zoom(zrz, 3)
    @test zrz2 == zr
    zrz2 = GtkReactive.zoom(zrz, 1.9)
    @test zrz2.currentview.x == 4..97
    @test zrz2.currentview.y == 3..78
    zrz = GtkReactive.zoom(zr, 0.5, GtkReactive.XY{DeviceUnit}(50.5, 40.5))
    @test zrz.currentview.x == 26..75
    @test zrz.currentview.y == 21..60
    zrz = GtkReactive.zoom(zr, 0.5, GtkReactive.XY{DeviceUnit}(60.5, 30.5))
    @test zrz.currentview.x == 31..80
    @test zrz.currentview.y == 16..55
    zrr = GtkReactive.reset(zrz)
    @test zrr == zr

    zrbb = ZoomRegion(zr.fullview, BoundingBox(5, 15, 35, 75))
    @test zrbb.fullview === zr.fullview
    @test zrbb.currentview.x == 5..15
    @test zrbb.currentview.y == 35..75
    @test typeof(zrbb.currentview) == typeof(zr.currentview)

    zrsig = Signal(zr)
    push!(zrsig, (3:5, 4:7))
    rr()
    zr = value(zrsig)
    @test zr.fullview.y == 1..80
    @test zr.fullview.x == 1..100
    @test zr.currentview.y == 3..5
    @test zr.currentview.x == 4..7
    push!(zrsig, XY(1..2, 3..4))
    rr()
    zr = value(zrsig)
    @test zr.fullview.y == 1..80
    @test zr.fullview.x == 1..100
    @test zr.currentview.y == 3..4
    @test zr.currentview.x == 1..2

    zr = ZoomRegion(Foo())
    @test zr.fullview.y == 1..7
    @test zr.fullview.x == 1..9

    zr = ZoomRegion((1:100, 1:80), (11:20, 8:12))
    @test zr.fullview.x == 1..80
    @test zr.fullview.y == 1..100
    @test zr.currentview.x == 8..12
    @test zr.currentview.y == 11..20
    @test axes(zr) == (11:20, 8:12)
end

### Simulate the mouse clicks, etc. to trigger zoom/pan
# Again, this doesn't seem to work inside a @testset
win = Window() |> (c = canvas(UserUnit))
zr = Signal(ZoomRegion((1:11, 1:20)))
zoomrb = init_zoom_rubberband(c, zr)
zooms = init_zoom_scroll(c, zr)
pans = init_pan_scroll(c, zr)
pand = init_pan_drag(c, zr)
draw(c) do cnvs
    set_coordinates(c, value(zr))
    fill!(c, colorant"blue")
end
Gtk.showall(win)
sleep(0.1)

# Zoom by rubber band
signal_emit(widget(c), "button-press-event", Bool,
            eventbutton(c, BUTTON_PRESS, 1, UserUnit(5), UserUnit(3), CONTROL))
sleep(0.1)
rr()
signal_emit(widget(c), "motion-notify-event", Bool,
            eventmotion(c, mask(1), UserUnit(10), UserUnit(4)))
sleep(0.1)
rr()
signal_emit(widget(c), "button-release-event", Bool,
            eventbutton(c, GtkReactive.BUTTON_RELEASE, 1, UserUnit(10), UserUnit(4)))
sleep(0.1)
rr()
sleep(0.1)
rr()
@test value(zr).currentview.x == 5..10
@test value(zr).currentview.y == 3..4
# Ensure that the rubber band damage has been repaired
fn = tempname()
Cairo.write_to_png(getgc(c).surface, fn)
imgout = load(fn)
rm(fn)
@test all(x->x==colorant"blue", imgout)

# Pan-drag
signal_emit(widget(c), "button-press-event", Bool,
            eventbutton(c, BUTTON_PRESS, 1, UserUnit(6), UserUnit(3), 0))
sleep(0.1)
rr()
signal_emit(widget(c), "motion-notify-event", Bool,
            eventmotion(c, mask(1), UserUnit(7), UserUnit(2)))
sleep(0.1)
rr()
sleep(0.1)
rr()
@test value(zr).currentview.x == 4..9
@test value(zr).currentview.y == 4..5

# Reset
signal_emit(widget(c), "button-press-event", Bool,
            eventbutton(c, DOUBLE_BUTTON_PRESS, 1, UserUnit(5), UserUnit(4.5), CONTROL))
sleep(0.1)
rr()
sleep(0.1)
rr()
@test value(zr).currentview.x == 1..20
@test value(zr).currentview.y == 1..11

# Zoom-scroll
signal_emit(widget(c), "scroll-event", Bool,
            eventscroll(c, UP, UserUnit(8), UserUnit(4), CONTROL))
sleep(0.1)
rr()
sleep(0.1)
rr()
@test value(zr).currentview.x == 4..14
@test value(zr).currentview.y == 2..8

# Pan-scroll
signal_emit(widget(c), "scroll-event", Bool,
            eventscroll(c, RIGHT, UserUnit(8), UserUnit(4), 0))
sleep(0.1)
rr()
sleep(0.1)
rr()
@test value(zr).currentview.x == 5..15
@test value(zr).currentview.y == 2..8

signal_emit(widget(c), "scroll-event", Bool,
            eventscroll(c, DOWN, UserUnit(8), UserUnit(4), 0))
sleep(0.1)
rr()
sleep(0.1)
rr()
@test value(zr).currentview.x == 5..15
@test value(zr).currentview.y == 3..9

destroy(win)

@testset "Surfaces" begin
    for (val, cmp) in ((0.2, Gray24(0.2)),
                       (Gray(N0f8(0.5)), Gray24(0.5)),
                       (RGB(0, 1, 0), RGB24(0, 1, 0)),
                       (RGBA(1, 0, 0.5, 0.8), ARGB32(1, 0, 0.5, 0.8)))
        surf = GtkReactive.image_surface(fill(val, 3, 5))
        @test surf.height == 3 && surf.width == 5
        @test all(x->x == reinterpret(UInt32, cmp), surf.data)
        destroy(surf)
    end
end

@testset "Layout" begin
    g = Grid()
    g[1,1] = textbox("hello")
end

# Ensure that the examples run (but the Reactive queue is stopped, so
# they won't work unless one calls `@async Reactive.run()` manually)
examplepath = joinpath(dirname(dirname(@__FILE__)), "examples")
include(joinpath(examplepath, "imageviewer.jl"))
include(joinpath(examplepath, "widgets.jl"))
include(joinpath(examplepath, "drawing.jl"))


text:
# custom types to represent some special parsing states a column might be in

# not just that a column type is Missing, but that we'll end up dropping the column
# eventually (i.e. implies col.willdrop) or the user specifically set the column
# type to Missing (i.e. col.userprovided && col.type === Missing)
# a column type may be Missing and we still want to try and detect a proper column type
# but in HardMissing case, we don't; we know the column will *always* be Missing
struct HardMissing end

function willdrop!(columns, i)
    @inbounds col = columns[i]
    col.willdrop = true
    col.type = HardMissing
    return
end

struct NeedsTypeDetection end

nonmissingtypeunlessmissingtype(::Type{T}) where {T} = ifelse(T === Missing, Missing, nonmissingtype(T))
pooltype(col) = nonmissingtype(keytype(col.refpool.refs))
finaltype(T) = T
finaltype(::Type{HardMissing}) = Missing
finaltype(::Type{NeedsTypeDetection}) = Missing
coltype(col) = ifelse(col.anymissing, Union{finaltype(col.type), Missing}, finaltype(col.type))

pooled(col) = col.pool == 1.0
maybepooled(col) = 0.0 < col.pool

function getpool(x::Real)::Float64
    if x isa Bool
        return x ? 1.0 : 0.0
    else
        y = Float64(x)
        (isnan(y) || 0.0 <= y <= 1.0) || throw(ArgumentError("pool argument must be in the range: 0.0 <= x <= 1.0"))
        return y
    end
end

tupcat(::Type{Tuple{}}, S) = Tuple{S}
tupcat(::Type{Tuple{T}}, S) where {T} = Tuple{T, S}
tupcat(::Type{Tuple{T, T2}}, S) where {T, T2} = Tuple{T, T2, S}
tupcat(::Type{Tuple{T, T2, T3}}, S) where {T, T2, T3} = Tuple{T, T2, T3, S}
tupcat(::Type{Tuple{T, T2, T3, T4}}, S) where {T, T2, T3, T4} = Tuple{T, T2, T3, T4, S}
tupcat(::Type{T}, S) where {T <: Tuple} = Tuple{Any[(fieldtype(T, i) for i = 1:fieldcount(T))..., S]...}

const StringTypes = Union{Type{String}, Type{PosLenString}, Type{<:InlineString}}
pickstringtype(T, maxstringsize) = T === InlineString ? (maxstringsize < DEFAULT_MAX_INLINE_STRING_LENGTH ? InlineStringType(maxstringsize) : String) : T

# we define our own bit flag on a Parsers.ReturnCode to signal if a column needs to promote to string
const PROMOTE_TO_STRING = 0b0100000000000000 % Int16
promote_to_string(code) = code & PROMOTE_TO_STRING > 0

isinttype(T) = T === Int8 || T === Int16 || T === Int32 || T === Int64 || T === Int128

@inline function promote_types(@nospecialize(T), @nospecialize(S))
    T === S && return T
    T === NeedsTypeDetection && return S
    S === NeedsTypeDetection && return T
    T === Missing && return S
    S === Missing && return T
    isinttype(T) && isinttype(S) && return promote_type(T, S)
    isinttype(T) && S === Float64 && return S
    T === Float64 && isinttype(S) && return T
    T <: InlineString && S <: InlineString && return promote_type(T, S)
    (T === String || S === String) && return String
    T <: InlineString && return T
    S <: InlineString && return S
    return nothing
end

# when users pass non-standard types, we need to keep track of them in a Tuple{...} to generate efficient custom parsing kernel codes
@inline function nonstandardtype(T)
    T = nonmissingtype(T)
    if T === Union{} ||
       T === NeedsTypeDetection ||
       T isa StringTypes ||
       isinttype(T) ||
       T === Float16 ||
       T === Float32 ||
       T === Float64 ||
       T === Bool ||
       T === Date ||
       T === DateTime ||
       T === Time
        return Union{}
    end
    return T
end

## column array allocating
# we don't want to use SentinelVector for small integer types due to the higher risk of
# sentinel value collision, so we just use Vector{Union{T, Missing}} and convert to Vector{T} if no missings were found
const SmallIntegers = Union{Int8, UInt8, Int16, UInt16, Int32, UInt32}

const SVec{T} = SentinelVector{T, T, Missing, Vector{T}}
const SVec2{T} = SentinelVector{T, typeof(undef), Missing, Vector{T}}

struct Pooled end

promotevectype(::Type{T}) where {T <: Union{Bool, SmallIntegers}} = vectype(T)
promotevectype(::Type{T}) where {T} = SentinelVector{T}

# allocate columns for a full file
function allocate!(columns, rowsguess)
    for i = 1:length(columns)
        @inbounds col = columns[i]
        # if the type hasn't been detected yet, then column will get allocated
        # in the detect method while parsing
        if col.type !== NeedsTypeDetection
            if maybepooled(col) && (col.type isa StringTypes || col.columnspecificpool)
                col.column = allocate(Pooled, rowsguess)
                if !isdefined(col, :refpool)
                    col.refpool = RefPool(col.type)
                end
            else
                col.column = allocate(col.type, rowsguess)
            end
        end
    end
    return
end

setmissing!(col, i) = col[i] = missing
const POSLEN_MISSING = PosLen(0, 0, true)
setmissing!(col::Vector{PosLen}, i) = col[i] = POSLEN_MISSING

@inline function allocate(T, len)
    if T === NeedsTypeDetection || T === HardMissing || T === Missing
        # MissingVector is an efficient representation in SentinelArrays.jl package
        return MissingVector(len)
    elseif T === PosLenString
        A = Vector{PosLen}(undef, len)
        memset!(pointer(A), typemax(UInt8), sizeof(A))
        return A
    elseif T === String
        return SentinelVector{String}(undef, len)
    elseif T === Pooled
        return fill(UInt32(1), len) # initialize w/ all missing values
    elseif T === Bool
        return Vector{Union{Missing, Bool}}(undef, len)
    elseif T <: SmallIntegers
        return Vector{Union{Missing, T}}(undef, len)
    else
        return SentinelVector{T}(undef, len)
    end
end

function reallocate!(@nospecialize(A), len)
    if A isa Vector{UInt32}
        oldlen = length(A)
        resize!(A, len)
        # make sure new values are initialized to missing
        for i = (oldlen + 1):len
            @inbounds A[i] = UInt32(1)
        end
    elseif A isa Vector{PosLen}
        oldlen = length(A)
        resize!(A, len)
        # when reallocating, we just need to make sure the missing bit is set for lazy string PosLen
        memset!(pointer(A, oldlen + 1), typemax(UInt8), (len - oldlen) * 8)
    else
        resize!(A, len)
    end
    return
end

firstarray(x::ChainedVector) = x.arrays[1]

vectype(::Type{T}) where {T <: Union{Bool, SmallIntegers}} = Vector{Union{T, Missing}}
vectype(::Type{T}) where {T} = isbitstype(T) ? SVec{T} : SVec2{T}
vectype(::Type{PosLenString}) = PosLenStringVector{Union{PosLenString, Missing}}
vectype(::Type{Pooled}) = Vector{UInt32}
pooledvectype(::Type{T}) where {T} = PooledVector{Union{T, Missing}, UInt32, Vector{UInt32}}
pooledtype(::Type{T}) where {T} = PooledVector{T, UInt32, Vector{UInt32}}
# missingvectype(::PooledVector{T, R, AT}) where {T, R, AT} = PooledVector{Union{T, Missing}, R, AT}

_promote(::Type{A}, x::A) where {A} = x
_promote(::Type{A}, x::MissingVector) where {A <: AbstractVector{T}} where {T} = allocate(Base.nonmissingtype(T), length(x))
_promote(::Type{A}, x::Vector) where {A <: SentinelArray{T}} where {T} = convert(SentinelVector{T}, x)
_promote(::Type{Vector{T}}, x::Vector) where {T} = convert(Vector{T}, x)
_promote(::Type{A}, x::SentinelVector) where {A <: SentinelVector{T}} where {T} = convert(SentinelVector{T}, x)
_promote(::Type{PosLenStringVector{Union{PosLenString, Missing}}}, x::PosLenStringVector{PosLenString}) = PosLenStringVector{Union{PosLenString, Missing}}(x.data, x.poslens, x.e)
_promote(::Type{PosLenStringVector{Union{PosLenString, Missing}}}, x::MissingVector) = PosLenStringVector{Union{PosLenString, Missing}}(UInt8[], fill(POSLEN_MISSING, x.len), 0x00)
_promote(::Type{PooledVector{Union{T, Missing}, R, RA}}, x::PooledVector) where {T, R, RA} =
    PooledArray(PooledArrays.RefArray(x.refs), convert(Dict{Union{T, Missing}, UInt32}, x.invpool), convert(Vector{Union{T, Missing}}, x.pool))
_promote(::Type{PooledVector{Union{T, Missing}, R, RA}}, x::MissingVector) where {T, R, RA} =
    PooledArray(PooledArrays.RefArray(fill(UInt32(1), length(x))), Dict{Union{T, Missing}, UInt32}(missing => UInt32(1)), Union{T, Missing}[missing])
_promote(::Type{PooledVector{T, R, RA}}, x) where {T, R, RA} = PooledArray{T}(x)
_promote(::Type{PooledVector{T, R, RA}}, x::PooledVector{T, R, RA}) where {T, R, RA} = x  # avoid ambiguity
_promote(::Type{PooledVector{T, R, RA}}, x::PooledVector{T, R, RA}) where {T>:Missing, R, RA} = x  # avoid ambiguity

function chaincolumns!(@nospecialize(a), @nospecialize(b))
    if a isa PooledArray || b isa PooledArray
        # special-case PooledArrays apart from other container types
        # because we want the outermost array to be PooledArray instead of ChainedVector
        if eltype(a) == eltype(b)
            if a isa PooledArray
                return append!(a, b)
            else
                px = _promote(pooledtype(eltype(a)), a)
                return append!(px, b)
            end
        elseif a isa MissingVector
            P = pooledvectype(eltype(b))
        elseif b isa MissingVector
            P = pooledvectype(eltype(a))
        elseif eltype(a) >: Missing || eltype(b) >: Missing
            P = pooledvectype(promote_types(Base.nonmissingtype(eltype(a)), Base.nonmissingtype(eltype(b))))
        else
            # both arrays are non-missing, but not same eltype, just need to promote
            P = pooledtype(promote_types(eltype(a), eltype(b)))
        end
        px = _promote(P, a)
        py = _promote(P, b)
        return append!(px, py)
    end
    c = firstarray(a)
    if typeof(c) == typeof(b)
        # easiest case; vector types match, so just append
        return append!(a, b)
    elseif c isa MissingVector
        A = vectype(Base.nonmissingtype(eltype(b)))
    elseif b isa MissingVector
        A = vectype(Base.nonmissingtype(eltype(c)))
    elseif c isa Vector && b isa SentinelVector
        A = vectype(promote_types(eltype(c), Base.nonmissingtype(eltype(b))))
    elseif c isa SentinelVector && b isa Vector
        A = vectype(promote_types(Base.nonmissingtype(eltype(c)), eltype(b)))
    elseif c isa PosLenStringVector{PosLenString} && b isa PosLenStringVector{Union{PosLenString, Missing}}
        A = typeof(b)
    elseif c isa PosLenStringVector{Union{PosLenString, Missing}} && b isa PosLenStringVector{PosLenString}
        A = typeof(c)
    elseif c isa Vector{Bool} && b isa Vector{Union{Bool, Missing}}
        A = typeof(b)
    elseif c isa Vector{Union{Bool, Missing}} && b isa Vector{Bool}
        A = typeof(c)
    elseif c isa Vector{<:SmallIntegers} && b isa Vector{<:Union{SmallIntegers, Missing}}
        A = typeof(b)
    elseif c isa Vector{<:Union{SmallIntegers, Missing}} && b isa Vector{<:SmallIntegers}
        A = typeof(c)
    elseif c isa Vector && b isa Vector
        # two vectors, but we know eltype doesn't match, so try to promote
        A = Vector{promote_types(eltype(c), eltype(b))}
    elseif c isa SentinelVector && b isa SentinelVector
        A = vectype(promote_types(Base.nonmissingtype(eltype(c)), Base.nonmissingtype(eltype(b))))
    end
    x = ChainedVector([_promote(A, x) for x in a.arrays])
    y = _promote(A, b)
    return append!(x, y)
end

# one-liner suggested from ScottPJones
consumeBOM(buf, pos) = (length(buf) >= 3 && buf[pos] == 0xef && buf[pos + 1] == 0xbb && buf[pos + 2] == 0xbf) ? pos + 3 : pos

# whatever input is given, turn it into an AbstractVector{UInt8} we can parse with
@inline function getbytebuffer(x, buffer_in_memory)
    tfile = nothing
    if x isa Vector{UInt8}
        return x, 1, length(x), tfile
    elseif x isa SubArray{UInt8, 1, Vector{UInt8}}
        return parent(x), first(x.indices[1]), last(x.indices[1]), tfile
    elseif x isa StringCodeUnits
        return unsafe_wrap(Vector{UInt8}, x.s), 1, length(x), tfile
    elseif x isa IOBuffer
        if x.data isa Vector{UInt8}
            return x.data, x.ptr, x.size, tfile
        elseif x.data isa SubArray{UInt8}
            x = x.data
            return parent(x), first(x.indices), last(x.indices), tfile
        end
    elseif x isa Cmd || x isa IO
        if buffer_in_memory
            buf = Base.read(x isa Cmd ? open(x) : x)
        else
            buf, tfile = buffer_to_tempfile(CodecZlib.TranscodingStreams.Noop(), x isa Cmd ? open(x) : x)
        end
        return buf, 1, length(buf), tfile
    else
        try
            buf = Mmap.mmap(string(x))
            return buf, 1, length(buf), tfile
        catch e
            # if we can't mmap, try just buffering the whole thing into a tempfile byte vector
            if buffer_in_memory
                buf = Base.read(x)
            else
                buf, tfile = buffer_to_tempfile(CodecZlib.TranscodingStreams.Noop(), open(x))
            end
            return buf, 1, length(buf), tfile
        end
    end
end

function getsource(@nospecialize(x), buffer_in_memory)
    buf, pos, len, tfile = getbytebuffer(x, buffer_in_memory)::Tuple{Vector{UInt8},Int,Int,Union{Nothing,String}}
    if length(buf) >= 2 && buf[1] == 0x1f && buf[2] == 0x8b
        # gzipped source, gunzip it
        if buffer_in_memory
            buf = transcode(GzipDecompressor, buf)
        else
            # 917; if we already buffered input to tempfile, make sure the compressed tempfile is
            # cleaned up since we're only passing the *uncompressed* tempfile up for removal post-parsing
            tfile1 = tfile === nothing ? nothing : tfile
            buf, tfile = buffer_to_tempfile(GzipDecompressor(), IOBuffer(buf))
            if tfile1 !== nothing
                rm(tfile1; force=true)
            end
        end
        pos = 1
        len = length(buf)
    end
    return buf, pos, len, tfile
end

@inline function buffer_to_tempfile(codec, x)
    file, output = mktemp()
    stream = CodecZlib.TranscodingStream(codec, output)
    Base.write(stream, x)
    close(stream)
    return Mmap.mmap(file), file
end

@inline function getname(x)
    if x isa AbstractVector{UInt8}
        return "<raw byte buffer: $(hash(x))>"
    elseif x isa IO
        return string("<", typeof(x), ": $(hash(x))>")
    else
        return string(x)
    end
end

# normalizing column name utilities
const RESERVED = Set(["local", "global", "export", "let",
    "for", "struct", "while", "const", "continue", "import",
    "function", "if", "else", "try", "begin", "break", "catch",
    "return", "using", "baremodule", "macro", "finally",
    "module", "elseif", "end", "quote", "do"])

normalizename(name::Symbol) = name
function normalizename(name::String)::Symbol
    uname = strip(Unicode.normalize(name))
    id = Base.isidentifier(uname) ? uname : map(c->Base.is_id_char(c) ? c : '_', uname)
    cleansed = string((isempty(id) || !Base.is_id_start_char(id[1]) || id in RESERVED) ? "_" : "", id)
    return Symbol(replace(cleansed, r"(_)\1+"=>"_"))
end

function makeunique(names)
    set = Set(names)
    length(set) == length(names) && return Symbol[Symbol(x) for x in names]
    nms = Symbol[]
    for nm in names
        if nm in nms
            k = 1
            newnm = Symbol("$(nm)_$k")
            while newnm in set || newnm in nms
                k += 1
                newnm = Symbol("$(nm)_$k")
            end
            nm = newnm
        end
        push!(nms, nm)
    end
    @assert length(names) == length(nms)
    return nms
end

getordefault(x::AbstractDict{String}, nm, i, def) = haskey(x, string(nm)) ? x[string(nm)] : def
getordefault(x::AbstractDict{Symbol}, nm, i, def) = haskey(x, nm) ? x[nm] : def
getordefault(x::AbstractDict{Int}, nm, i, def) = haskey(x, i) ? x[i] : def
getordefault(x::AbstractDict, nm, i, def) = haskey(x, i) ? x[i] : haskey(x, nm) ? x[nm] : haskey(x, string(nm)) ? x[string(nm)] : def

# given a DateFormat, is it meant for parsing Date, DateTime, or Time?
function timetype(df::Parsers.Format)::Union{Type{Date}, Type{Time}, Type{DateTime}}
    date = false
    time = false
    for token in df.tokens
        T = typeof(token)
        if T in (Dates.DatePart{'H'}, Dates.DatePart{'I'}, Dates.DatePart{'M'}, Dates.DatePart{'S'}, Dates.DatePart{'s'})
            time = true
        elseif T in (Dates.DatePart{'y'}, Dates.DatePart{'Y'}, Dates.DatePart{'m'}, Dates.DatePart{'d'}, Dates.DatePart{'u'}, Dates.DatePart{'U'})
            date = true
        end
    end
    return ifelse(date & time, DateTime, ifelse(time, Time, Date))
end

# if a cell value of a csv file has escape characters, we need to unescape it
function unescape(s, e)
    n = ncodeunits(s)
    buf = Base.StringVector(n)
    len = 1
    i = 1
    @inbounds begin
        while i <= n
            b = codeunit(s, i)
            if b == e
                i += 1
                b = codeunit(s, i)
            end
            @inbounds buf[len] = b
            len += 1
            i += 1
        end
    end
    resize!(buf, len - 1)
    return String(buf)
end

"""
    CSV.detect(str::String)

Use the same logic used by `CSV.File` to detect column types, to parse a value from a plain string.
This can be useful in conjunction with the `CSV.Rows` type, which returns each cell of a file as a String.
The order of types attempted is: `Int`, `Float64`, `Date`, `DateTime`, `Bool`, and if all fail, the input String is returned.
No errors are thrown.
For advanced usage, you can pass your own `Parsers.Options` type as a keyword argument `option=ops` for sentinel value detection.
"""
function detect end

@inline pass(code, tlen, x, typecode) = (code, tlen, x, typecode)

function detect(str::String; opts=Parsers.OPTIONS)
    code, tlen, x, xT = detect(pass, codeunits(str), 1, sizeof(str), opts, true)
    return something(x, str)
end

const DetectTypes = Union{Missing, Int8, Int16, Int32, Int64, Float64, Date, DateTime, Time, Bool, Nothing, PosLen}

const NEEDSTYPEDETECTION = 0x01
const HARDMISSING = 0x02
const MISSING = 0x03
const INT8 = 0x04
const INT16 = 0x05
const INT32 = 0x06
const INT64 = 0x07
const FLOAT64 = 0x08
const DATE = 0x09
const DATETIME = 0x0a
const TIME = 0x0b
const BOOL = 0x0c
const STRING = 0x0d
const TYPES = Union{Nothing, Type}[NeedsTypeDetection, HardMissing, Missing, Int8, Int16, Int32, Int64, Float64, Date, DateTime, Time, Bool, nothing]
isinttypecode(T) = T === INT8 || T === INT16 || T === INT32 || T === INT64
promote_typecode(a, b) = max(a, b)
typecode(@nospecialize(T)) = T === Missing ? MISSING : T === Int8 ? INT8 : T === Int16 ? INT16 : T === Int32 ? INT32 : T === Int64 ? INT64 :
                           T === Float64 ? FLOAT64 : T === Date ? DATE : T === DateTime ? DATETIME : T === Time ? TIME : T === Bool ? BOOL :
                           T === NeedsTypeDetection ? NEEDSTYPEDETECTION : T === HardMissing ? HARDMISSING : STRING

concrete_or_concreteunion(T) = isconcretetype(T) ||
    (T isa Union && concrete_or_concreteunion(T.a) && concrete_or_concreteunion(T.b))

@inline smallestint(cb, code, tlen, x) = x < typemax(Int8) ? cb(code, tlen, unsafe_trunc(Int8, x), INT8) : x < typemax(Int16) ? cb(code, tlen, unsafe_trunc(Int16, x), INT16) : x < typemax(Int32) ? cb(code, tlen, unsafe_trunc(Int32, x), INT32) : cb(code, tlen, x, INT64)
_widen(T) = widen(T)
_widen(::Type{Int128}) = Float64
_widen(::Type{Float64}) = nothing

@noinline function _parseany(T, buf, pos, len, opts)::Parsers.Result{Any}
    return Parsers.xparse(T, buf, pos, len, opts, Any)
end

@inline function detect(cb, buf, pos, len, opts, ensure_full_buf_consumed=true, downcast=false, row=0, col=0)
    int = Parsers.xparse(Int, buf, pos, len, opts)
    code, tlen = int.code, int.tlen
    if Parsers.invalidquotedfield(code)
        fatalerror(buf, pos, tlen, code, row, col)
    end
    if Parsers.sentinel(code) && code > 0
        return cb(code, tlen, missing, NEEDSTYPEDETECTION)
    end
    if Parsers.ok(code) && (!ensure_full_buf_consumed || (ensure_full_buf_consumed == ((pos + tlen - 1) == len)))
        return downcast ? smallestint(cb, code, tlen, int.val) : cb(code, tlen, int.val, Int === Int64 ? INT64 : INT32)
    end
    float = Parsers.xparse(Float64, buf, pos, len, opts)
    code, tlen = float.code, float.tlen
    if Parsers.ok(code) && (!ensure_full_buf_consumed || (ensure_full_buf_consumed == ((pos + tlen - 1) == len)))
        return cb(code, tlen, float.val, FLOAT64)
    end
    if opts.dateformat === nothing
        date = Parsers.xparse(Date, buf, pos, len, opts, Date)
        code, tlen = date.code, date.tlen
        if Parsers.ok(code) && (!ensure_full_buf_consumed || (ensure_full_buf_consumed == ((pos + tlen - 1) == len)))
            return cb(code, tlen, date.val, DATE)
        end
        datetime = Parsers.xparse(DateTime, buf, pos, len, opts)
        code, tlen = datetime.code, datetime.tlen
        if Parsers.ok(code) && (!ensure_full_buf_consumed || (ensure_full_buf_consumed == ((pos + tlen - 1) == len)))
            return cb(code, tlen, datetime.val, DATETIME)
        end
        time = Parsers.xparse(Time, buf, pos, len, opts)
        code, tlen = time.code, time.tlen
        if Parsers.ok(code) && (!ensure_full_buf_consumed || (ensure_full_buf_consumed == ((pos + tlen - 1) == len)))
            return cb(code, tlen, time.val, TIME)
        end
    else
        # use user-provided dateformat
        DT = timetype(opts.dateformat)
        if DT === Date
            dt = Parsers.xparse(DT, buf, pos, len, opts)
            code, tlen = dt.code, dt.tlen
            if Parsers.ok(code) && (!ensure_full_buf_consumed || (ensure_full_buf_consumed == ((pos + tlen - 1) == len)))
                return cb(code, tlen, dt.val, DATE)
            end
        elseif DT === Time
            dt2 = Parsers.xparse(DT, buf, pos, len, opts)
            code, tlen = dt2.code, dt2.tlen
            if Parsers.ok(code) && (!ensure_full_buf_consumed || (ensure_full_buf_consumed == ((pos + tlen - 1) == len)))
                return cb(code, tlen, dt2.val, TIME)
            end
        elseif DT === DateTime
            dt3 = Parsers.xparse(DT, buf, pos, len, opts)
            code, tlen = dt3.code, dt3.tlen
            if Parsers.ok(code) && (!ensure_full_buf_consumed || (ensure_full_buf_consumed == ((pos + tlen - 1) == len)))
                return cb(code, tlen, dt3.val, DATETIME)
            end
        end
    end
    bool = Parsers.xparse(Bool, buf, pos, len, opts)
    code, tlen = bool.code, bool.tlen
    if Parsers.ok(code) && (ensure_full_buf_consumed == ((pos + tlen - 1) == len))
        return cb(code, tlen, bool.val, BOOL)
    end
    return cb(code, tlen, nothing, STRING)
end

# a ReversedBuf takes a byte vector and indexes backwards;
# used for the footerskip keyword argument, which starts at the bottom of the file
# and skips lines backwards
struct ReversedBuf <: AbstractVector{UInt8}
    buf::Vector{UInt8}
end

Base.size(a::ReversedBuf) = size(a.buf)
Base.IndexStyle(::Type{ReversedBuf}) = Base.IndexLinear()
Base.getindex(a::ReversedBuf, i::Int) = a.buf[end + 1 - i]
Base.pointer(a::ReversedBuf, pos::Integer=1) = pointer(a.buf, length(a.buf) + 1 - pos)

memset!(ptr, value, num) = ccall(:memset, Ptr{Cvoid}, (Ptr{Cvoid}, Cint, Csize_t), ptr, value, num)

struct PointerString
    ptr::Ptr{UInt8}
    len::Int
end

function Base.hash(s::PointerString, h::UInt)
    h += Base.memhash_seed
    ccall(Base.memhash, UInt, (Ptr{UInt8}, Csize_t, UInt32), s.ptr, s.len, h % UInt32) + h
end

import Base: ==
function ==(x::AbstractString, y::PointerString)
    sizeof(x) == y.len && ccall(:memcmp, Cint, (Ptr{UInt8}, Ptr{UInt8}, Csize_t), pointer(x), y.ptr, y.len) == 0
end
function ==(x::PointerString, y::PointerString)
    x.len == y.len && ccall(:memcmp, Cint, (Ptr{UInt8}, Ptr{UInt8}, Csize_t), x.ptr, y.ptr, y.len) == 0
end
==(y::PointerString, x::AbstractString) = x == y

Base.ncodeunits(s::PointerString) = s.len
@inline function Base.codeunit(s::PointerString, i::Integer)
    @boundscheck checkbounds(s, i)
    GC.@preserve s unsafe_load(s.ptr + i - 1)
end

_unsafe_string(p, len) = ccall(:jl_pchar_to_string, Ref{String}, (Ptr{UInt8}, Int), p, len)
Base.String(x::PointerString) = _unsafe_string(x.ptr, x.len)
WeakRefStrings.PosLenString(x::PointerString) = PosLenString(unsafe_wrap(Array, x.ptr, x.len), PosLen(1, x.len), 0x00)

struct Arg
    x::Any
end
Base.getindex(x::Arg) = x.x

macro refargs(ex)
    ex isa Expr || throw(ArgumentError("must pass an expressiong to @refargs"))
    (ex.head == :call || ex.head == :function) || throw(ArgumentError("@refargs ex must be function call or definition"))
    if ex.head == :call
        for i = 2:length(ex.args)
            ex.args[i] = Expr(:call, :(CSV.Arg), ex.args[i])
        end
        return esc(ex)
    else # ex.head == :function
        refs = Expr(:block)
        fargs = ex.args[1].args
        for i = 2:length(fargs)
            arg = fargs[i]
            (arg isa Symbol || arg.head == :(::)) || throw(ArgumentError("unsupported argument expression: `$arg`"))
            nm = arg isa Symbol ? arg : arg.args[1]
            T = arg isa Symbol ? :Any : arg.args[2]
            push!(refs.args, Expr(:(=), Expr(:(::), nm, T), Expr(:ref, nm)))
            fargs[i] = Expr(:(::), nm, :(CSV.Arg))
        end
        pushfirst!(ex.args[2].args, refs)
        return esc(ex)
    end
end

ext: jl


text:
using MLPreprocessing
using DataFrames
using Base.Test

tests = [
    "tst_expand.jl"
    "tst_center.jl"
    "tst_standardize.jl"
    "tst_fixedrangescaler.jl"
]

for t in tests
    @testset "$t" begin
        include(t)
    end
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

"""
    GitConfig(path::AbstractString, level::Consts.GIT_CONFIG=Consts.CONFIG_LEVEL_APP, force::Bool=false)

Create a new `GitConfig` by loading configuration information from the file at
`path`. See [`addfile`](@ref) for more information about the `level` and `force`
options.
"""
function GitConfig(path::AbstractString,
                   level::Consts.GIT_CONFIG = Consts.CONFIG_LEVEL_APP,
                   force::Bool=false)
    # create new config object
    cfg_ptr_ptr = Ref{Ptr{Cvoid}}(C_NULL)
    @check ccall((:git_config_new, :libgit2), Cint, (Ptr{Ptr{Cvoid}},), cfg_ptr_ptr)
    cfg = GitConfig(cfg_ptr_ptr[])
    try
        addfile(cfg, path, level, force)
    catch ex
        close(cfg)
        rethrow(ex)
    end
    return cfg
end

"""
    GitConfig(repo::GitRepo)

Get the stored configuration for the git repository `repo`. If `repo` does not
have a specific configuration file set, the default git configuration will be
used.
"""
function GitConfig(repo::GitRepo)
    cfg_ptr_ptr = Ref{Ptr{Cvoid}}(C_NULL)
    @check ccall((:git_repository_config, :libgit2), Cint,
                  (Ptr{Ptr{Cvoid}}, Ptr{Cvoid}), cfg_ptr_ptr, repo.ptr)
    return GitConfig(repo, cfg_ptr_ptr[])
end

"""
    GitConfig(level::Consts.GIT_CONFIG=Consts.CONFIG_LEVEL_DEFAULT)

Get the default git configuration by loading the global and system configuration
files into a prioritized configuration. This can be used to access default configuration
options outside a specific git repository.
"""
function GitConfig(level::Consts.GIT_CONFIG = Consts.CONFIG_LEVEL_DEFAULT)
    cfg_ptr_ptr = Ref{Ptr{Cvoid}}(C_NULL)
    @check ccall((:git_config_open_default, :libgit2), Cint,
                  (Ptr{Ptr{Cvoid}},), cfg_ptr_ptr)
    cfg = GitConfig(cfg_ptr_ptr[])
    if level != Consts.CONFIG_LEVEL_DEFAULT
        glb_cfg_ptr_ptr = Ref{Ptr{Cvoid}}(C_NULL)
        tmpcfg = cfg
        try
            @check ccall((:git_config_open_level, :libgit2), Cint,
                         (Ptr{Ptr{Cvoid}}, Ptr{Cvoid}, Cint),
                          glb_cfg_ptr_ptr, cfg.ptr, Cint(level))
            cfg = GitConfig(glb_cfg_ptr_ptr[])
        finally
            close(tmpcfg)
        end
    end
    return cfg
end

"""
    addfile(cfg::GitConfig, path::AbstractString, level::Consts.GIT_CONFIG=Consts.CONFIG_LEVEL_APP, force::Bool=false)

Add an existing git configuration file located at `path` to the current
`GitConfig` `cfg`. If the file does not exist, it will be created.
`level` sets the git configuration priority level and is determined by
[`Consts.GIT_CONFIG`](@ref). If `force` is `false` and a configuration for
the given priority level already exists, `addfile` will error. If `force` is
`true`, the existing configuration will be replaced by the one in the file at
`path`.
"""
function addfile(cfg::GitConfig, path::AbstractString,
                 level::Consts.GIT_CONFIG = Consts.CONFIG_LEVEL_APP,
                 force::Bool=false)
    @check ccall((:git_config_add_file_ondisk, :libgit2), Cint,
                  (Ptr{Ptr{Cvoid}}, Cstring, Cint, Cint),
                   cfg.ptr, path, Cint(level), Cint(force))
end

function get(::Type{<:AbstractString}, c::GitConfig, name::AbstractString)
    buf_ref = Ref(Buffer())
    @check ccall((:git_config_get_string_buf, :libgit2), Cint,
                 (Ptr{Buffer}, Ptr{Cvoid}, Cstring), buf_ref, c.ptr, name)
    buf = buf_ref[]
    str = unsafe_string(buf.ptr, buf.size)
    free(buf_ref)
    str
end

function get(::Type{Bool}, c::GitConfig, name::AbstractString)
    val_ptr = Ref(Cint(0))
    @check ccall((:git_config_get_bool, :libgit2), Cint,
          (Ptr{Cint}, Ptr{Cvoid}, Cstring), val_ptr, c.ptr, name)
    return Bool(val_ptr[])
end

function get(::Type{Int32}, c::GitConfig, name::AbstractString)
    val_ptr = Ref(Cint(0))
    @check ccall((:git_config_get_int32, :libgit2), Cint,
          (Ptr{Cint}, Ptr{Cvoid}, Cstring), val_ptr, c.ptr, name)
    return val_ptr[]
end

function get(::Type{Int64}, c::GitConfig, name::AbstractString)
    val_ptr = Ref(Cintmax_t(0))
    @check ccall((:git_config_get_int64, :libgit2), Cint,
          (Ptr{Cintmax_t}, Ptr{Cvoid}, Cstring), val_ptr, c.ptr, name)
    return val_ptr[]
end

function get(c::GitConfig, name::AbstractString, default::T) where T
    res = default
    try res = get(T,c,name) end
    return res
end

function getconfig(r::GitRepo, name::AbstractString, default)
    with(GitConfig, r) do cfg
        get(cfg, name, default)
    end
end

function getconfig(rname::AbstractString, name::AbstractString, default)
    with(GitRepo, rname) do r
        with(GitConfig, r) do cfg
            get(cfg, name, default)
        end
    end
end

function getconfig(name::AbstractString, default)
    with(GitConfig) do cfg
        get(cfg, name, default)
    end
end

function set!(c::GitConfig, name::AbstractString, value::AbstractString)
    @check ccall((:git_config_set_string, :libgit2), Cint,
                  (Ptr{Cvoid}, Cstring, Cstring), c.ptr, name, value)
end

function set!(c::GitConfig, name::AbstractString, value::Bool)
    bval = Int32(value)
    @check ccall((:git_config_set_bool, :libgit2), Cint,
                  (Ptr{Cvoid}, Cstring, Cint), c.ptr, name, bval)
end

function set!(c::GitConfig, name::AbstractString, value::Int32)
    @check ccall((:git_config_set_int32, :libgit2), Cint,
                  (Ptr{Cvoid}, Cstring, Cint), c.ptr, name, value)
end

function set!(c::GitConfig, name::AbstractString, value::Int64)
    @check ccall((:git_config_set_int64, :libgit2), Cint,
                  (Ptr{Cvoid}, Cstring, Cintmax_t), c.ptr, name, value)
end

function GitConfigIter(cfg::GitConfig)
    ci_ptr = Ref{Ptr{Cvoid}}(C_NULL)
    @check ccall((:git_config_iterator_new, :libgit2), Cint,
                  (Ptr{Ptr{Cvoid}}, Ptr{Cvoid}), ci_ptr, cfg.ptr)
    return GitConfigIter(ci_ptr[])
end

function GitConfigIter(cfg::GitConfig, name::AbstractString)
    ci_ptr = Ref{Ptr{Cvoid}}(C_NULL)
    @check ccall((:git_config_multivar_iterator_new, :libgit2), Cint,
                  (Ptr{Ptr{Cvoid}}, Ptr{Cvoid}, Cstring, Cstring),
                  ci_ptr, cfg.ptr, name, C_NULL)
    return GitConfigIter(ci_ptr[])
end

function GitConfigIter(cfg::GitConfig, name::AbstractString, value::Regex)
    ci_ptr = Ref{Ptr{Cvoid}}(C_NULL)
    @check ccall((:git_config_multivar_iterator_new, :libgit2), Cint,
                  (Ptr{Ptr{Cvoid}}, Ptr{Cvoid}, Cstring, Cstring),
                  ci_ptr, cfg.ptr, name, value.pattern)
    return GitConfigIter(ci_ptr[])
end

function GitConfigIter(cfg::GitConfig, name::Regex)
    ci_ptr = Ref{Ptr{Cvoid}}(C_NULL)
    @check ccall((:git_config_iterator_glob_new, :libgit2), Cint,
                  (Ptr{Ptr{Cvoid}}, Ptr{Cvoid}, Cstring),
                  ci_ptr, cfg.ptr, name.pattern)
    return GitConfigIter(ci_ptr[])
end

function Base.start(ci::GitConfigIter)
    entry_ptr_ptr = Ref{Ptr{ConfigEntry}}(C_NULL)
    err = ccall((:git_config_next, :libgit2), Cint,
                 (Ptr{Ptr{ConfigEntry}}, Ptr{Cvoid}), entry_ptr_ptr, ci.ptr)
    if err == Cint(Error.GIT_OK)
        state = unsafe_load(entry_ptr_ptr[])
    elseif err == Cint(Error.ITEROVER)
        state = nothing
    else
        throw(GitError(err))
    end
    return state
end

Base.done(ci::GitConfigIter, state) = state === nothing

function Base.next(ci::GitConfigIter, state)
    entry = notnothing(state)
    entry_ptr_ptr = Ref{Ptr{ConfigEntry}}(C_NULL)
    err = ccall((:git_config_next, :libgit2), Cint,
                 (Ptr{Ptr{ConfigEntry}}, Ptr{Cvoid}), entry_ptr_ptr, ci.ptr)
    if err == Cint(Error.GIT_OK)
        state = unsafe_load(entry_ptr_ptr[])
    elseif err == Cint(Error.ITEROVER)
        state = nothing
    else
        throw(GitError(err))
    end
    return (entry, state)
end

Base.IteratorSize(::Type{GitConfigIter}) = Base.SizeUnknown()

max_stars_repo_name: djsegal/julia-fork

max_issues_repo_name: djsegal/julia-fork

max_forks_repo_name: djsegal/julia-fork


text:
# Supporting parameters and functions for skew-normal simulation example

# Stuff for mixture weights
m = 10 # maximum number of components
a = 1/m # Gamma parameters
G = Gamma(a,1/1) # prior on weights (before conditioning on s>0)
lambda = 1.0  # Poisson parameter of the limiting prior distribution on k as m -> Inf
c = quantile(G,1-lambda/m)  # cutoff (Note: This choice makes p(k) \propto Binomial(k|m,lambda/m)I(k>0).)
sigma = 0.25 # scale for weight proposals
log_v_prior(v_i) = logpdf(G,v_i) # prior on latent weights v_i

# Stuff for mixture components
normpdf(x,m,l) = sqrt(l/(2*pi))*exp(-0.5*l*(x-m).*(x-m))
normlogpdf(x,m,s) = -0.5*log(2*pi) - log(s) - 0.5*(x-m)*(x-m)/(s*s)
likelihood(x,t) = normpdf(x,t[1],t[2])  # t=[mean,precision]
# prior (base measure) on component params, [mean,log(precision)]
m0m,s0m = 0,5
m0l,s0l = 0,2
log_theta_prior(t) = normlogpdf(t[1],m0m,s0m) + normlogpdf(log(t[2]),m0l,s0l)
# proposal distribution for component param moves
sm,sl = 0.2*s0m, 0.2
theta_prop(t) = Float64[t[1]+sm*randn(), t[2]*exp(randn()*sl)]
log_theta_prop(t,tp) = normlogpdf(tp[1],t[1],sm) + normlogpdf(log(tp[2]),log(t[2]),sl)
new_thetas(m) = [Float64[m0m,exp(m0l)] for i=1:m]  # initialize params for m new components

# Densities, etc. for skew-normal
normpdf(x) = exp(-0.5*x.*x)/sqrt(2.0*pi)
normcdf(x) = 0.5*(1.0+erf(x/sqrt(2.0)))
skewpdf(x,loc,scale,shape) = (y=(x-loc)/scale; 2.0*normpdf(y).*normcdf(shape*y)/scale)
skewrnd(n,loc,scale,shape) = [(shape*z<randn() ? loc-scale*z : loc+scale*z) for z=randn(n)]

# Params of skew-normal mixture to use in simulation
p1 = 0.5  # probability of component 1
l1,s1,a1 = -4,1,5  # location, scale, and shape of component 1
l2,s2,a2 = -1,2,5  # location, scale, and shape of component 2

# Sampler code
include("core.jl")

max_stars_repo_name: jwmi/CoarsenedPosterior

max_issues_repo_name: jwmi/CoarsenedPosterior

max_forks_repo_name: jwmi/CoarsenedPosterior


text:
# The target is to develop an efficient FEM for quadrilateral mesh, for research purposes. Here, the quadrilaternal mesh is axis parallel, so it is easy to describe its elements, boundary and neighboring geometry
# for more complicated meshing, one may resort to third party software to get the geometry. Once geometry obtained, FEM is easy.

using SparseArrays
using LinearAlgebra
using Logging
import Base.Threads: nthreads, @threads

# some lesson: exact calculation of integral of tent functions leads to best accuracy!
# can either pull out some functions (e.g. a, f, or u) out of the integral and then use exact formula for tent functions, or can write a, f, u as linear combination of tent functions and obtain exact integral for the interactions between tent functions
# the latter approach is safer. The former approach needs to pay attention to boundaries, where basis functions change
# thus, combine accuracy and convenience: use stiffness and mass matrices, but for the boundary condition, can use the latter approach (do not need to introduce boundary element; this is why we always implement neuman and robin boundary condition in an open interval, in which case there is no boundary of boundary issues)

struct FEM_2dUnifQuadMesh{Ti,Tf}
    Ne::Ti # number of elements in each dimension
    Grid_x::Vector{Tf} # uniform grid in x axis, boundary included
    Grid_y::Vector{Tf} # uniform grid in y axis, boundary included
    Bdy_coordinates::Matrix{Tf} # bdy coordinates, used for getting boundary data
    Bdy_indices::Vector{Ti} # bdy index among all the nodes; ordering: one runs over x first then y
    ElemNode_loc2glo::Function # local index of an element to global index of its node; one runs over i first then j
end

# store the stiffness and mass matrix
struct FEM_store{Ti,Tf}
    A::SparseMatrixCSC{Tf,Ti}
    M::SparseMatrixCSC{Tf,Ti}
end

# constructor
function FEM_2dUnifQuadMesh(Ne)
    x = collect(LinRange(0, 1, Ne+1))
    y = copy(x)

    Bdy_coordinates = vcat(
        [x[1:end] y[1]*ones(Ne+1)],
        [x[1]*ones(Ne) y[2:end]],
        [x[end]*ones(Ne) y[2:end]],
        [x[2:end-1] y[end]*ones(Ne-1)]
    )

    Bdy_indices = reduce(
        vcat,collect.(
            [1:Ne+1,
            Ne+2:Ne+1:(Ne+1)*(Ne+1),
            2*(Ne+1):Ne+1:(Ne+1)*(Ne+1),
            Ne*Ne+Ne+2:(Ne+1)*(Ne+1)-1]
        )
    )

    function ElementNode_loc2glo(N, i, j, ind_node)
        # N: number of elements in one direction
        # i,j is the location of the element
        # edge marks the index of the edge (from 1 to 4, from bottom to top)
        if ind_node <= 2 # bottom
            global_idx = (N+1)*(j-1) + i + ind_node - 1;
        else # top
            global_idx = (N+1)*j + i + 4 - ind_node;
        end
        return global_idx
    end

    @info "[Mesh generation] mesh generated, $(Ne+1) nodes in each dimension"
    return FEM_2dUnifQuadMesh(Ne,x,y,Bdy_coordinates, Bdy_indices, ElementNode_loc2glo)
end

function FEM_StiffnMassAssembly(FEMparam,PDEparam)
    d = 2
    Ne = FEMparam.Ne

    # for sparse stiffness mtx A, and mass matrix M construction
    Irow= zeros(16*Ne^2) # 16*Ne^2: Ne^2 elements, each one leads to 4^2 edge interactions
    Jcol= copy(Irow)
    Aval = copy(Irow)
    Mval = copy(Irow)

    println("[multithreading] using ", Threads.nthreads(), " threads")
    # run over all elements
    for i = 1:Ne
        @threads for j = 1:Ne
            # local stiffness and local mass matrix
            local_A, local_M = FEM_LocalAssembly(FEMparam,PDEparam, i, j)

            for p = 1:4
                global_p = FEMparam.ElemNode_loc2glo(Ne, i, j, p);
                for q = 1:4
                    index = 16*Ne*(i-1)+16*(j-1)+4*(p-1)+q
                    global_q = FEMparam.ElemNode_loc2glo(Ne, i, j, q);
                    Irow[index] = global_p
                    Jcol[index] = global_q
                    Aval[index] = local_A[p, q]
                    Mval[index] = local_M[p, q]
                end
            end
        end
    end
    A = sparse(Irow,Jcol,Aval,(Ne+1)^d,(Ne+1)^d)
    M = sparse(Irow,Jcol,Mval,(Ne+1)^d,(Ne+1)^d)
    @info "[Assembly] finish assembly of stiffness and mass matrice"
    return FEM_store(A,M)
end

function FEM_BdyRhsAssembly(FEMparam,PDEparam,FEMstore)

    A = copy(FEMstore.A) # stiffness matrix for solving the systen
    M = FEMstore.M

    # right hand side
    Ne = FEMparam.Ne
    x = FEMparam.Grid_x
    y = FEMparam.Grid_y
    F = [PDEparam.rhs(x[i],y[j]) for j in 1:Ne+1 for i in 1:Ne+1]
    F = M * F

    # boundary location
    bdy_loc = FEMparam.Bdy_indices
    bdy_points = FEMparam.Bdy_coordinates
    bdy_type = PDEparam.bdy_type.(bdy_points[:,1],bdy_points[:,2])

    # Dirichlet boundary
    Diri_loc = findall(x->x==1,bdy_type)
    if length(Diri_loc) > 0
        Diri_bdy = PDEparam.bdy_Diri.(bdy_points[Diri_loc,1],bdy_points[Diri_loc,2])
        Diri_loc = bdy_loc[Diri_loc] # global index
        A[Diri_loc,:] .= 0
        F[Diri_loc] .= Diri_bdy
        A[Diri_loc,Diri_loc] .= sparse(I, length(Diri_loc),length(Diri_loc));
    end

    Neum_loc = findall(x->x==2,bdy_type)
    if length(Neum_loc) > 0
        Neum_bdy = PDEparam.bdy_Neum.(bdy_points[Neum_loc,1],bdy_points[Neum_loc,2])
        Neum_loc = bdy_loc[Neum_loc] # global index
        F[Neum_loc] += 1/(Ne)*Neum_bdy
    end

    Robin_loc = findall(x->x==3,bdy_type)
    if length(Robin_loc) > 0
        Robin_bdy = [PDEparam.bdy_Robin.(bdy_points[Robin_loc[i],1],bdy_points[Robin_loc[i],2]) for i in 1:length(Robin_loc)]
        # Robin_bdy = reduce(vcat, Robin_bdy)

        Robin_bdy1 = vcat([Robin_bdy[i][1] for i in 1:length(Robin_loc)])
        Robin_bdy2 = vcat([Robin_bdy[i][2] for i in 1:length(Robin_loc)])
        Robin_loc = bdy_loc[Robin_loc] # global index
        A[Robin_loc,Robin_loc] = A[Robin_loc,Robin_loc] + diagm(Robin_bdy1)*1/(Ne)
        F[Robin_loc] += 1/(Ne)*Robin_bdy2
    end

    @info "[Assembly] finish incorporating boundary data"
    return A, F
end

function FEM_LocalAssembly(FEMparam, PDEparam, i, j)
    # the stiffness matrix required locally to compute local nodal basis
    # for the bilinear boundary value with 1 at node
    Ne = FEMparam.Ne
    xlow = FEMparam.Grid_x[i];
    xhigh = FEMparam.Grid_x[i+1];
    ylow = FEMparam.Grid_y[j];
    yhigh = FEMparam.Grid_y[j+1];

    xmid = (xlow+xhigh)/2;
    ymid = (ylow+yhigh)/2;

    local_A = zeros(4,4)
    local_M = copy(local_A)

    # mid-point of a, then exact quadrature for tent functions
    # f is written as sum of nodal values * tent functions
    for i = 1:4
        for j = 1:4
            if i==j
                local_A[i,j] = 2/3*PDEparam.a(xmid,ymid);
                local_M[i,j] = 1/9/Ne^2 # exact
            elseif i==j+2 || i==j-2
                local_A[i,j] = -1/3*PDEparam.a(xmid,ymid);
                local_M[i,j] = 1/36/Ne^2
            else
                local_A[i,j] = -1/6*PDEparam.a(xmid,ymid);
                local_M[i,j] = 1/18/Ne^2
            end
        end
    end

    return local_A, local_M
end

# domain [0,1]*[0,1], Ne is the number of interior DOF in each dimension
function FEM_Solver(FEMparam,PDEparam)
    # FEM for VarElliptic, uniform quadrilaternal mesh, so we only need one parameter Ne that is the number of interior DOF in each dimension, to fully characterize the mesh
    FEMstore = FEM_StiffnMassAssembly(FEMparam,PDEparam)
    A, F = FEM_BdyRhsAssembly(FEMparam,PDEparam,FEMstore)
    sol = A\F
    @info "[Linear solver] linear system solved"

    # result=reshape(sol,Ne+1,Ne+1)'
    # the ordering is that x moves first y second
    return sol, FEMstore
end

## subsequent solve
function FEM_SubsequentSolve(FEMparam,PDEparam,FEMstore)
    A, F = FEM_BdyRhsAssembly(FEMparam,PDEparam,FEMstore)
    sol = A\F
    @info "[Linear solver] linear system solved"
    return sol
end


text:
using TimeseriesPrediction
using Test
import Random

Random.seed!(42)
include("system_defs.jl")

@testset "Coupled Henon 1D" begin
    M = 100
    N_train = 2000
    p = 5
    U, V = coupled_henon1D(M,N_train+p+100, rand(M), rand(M))
    #Reconstruct this #
    utrain = U[1:N_train]
    vtrain = V[1:N_train]
    utest  = U[N_train:N_train+p]
    vtest  = V[N_train:N_train+p]

    @testset "ALM γ=$γ, B=$B" for γ=2:3, B=1:2
        em = cubic_shell_embedding(utrain,γ,1,B,1,ConstantBoundary(10.))
        upred = temporalprediction(utrain,em, p)

        @test upred[1] == utrain[end]
        @test sum(abs.(utest[2]-upred[2]))/M/p < 0.05
        ε = [sum(abs.(utest[i]-upred[i])) for i=1:p+1]
        @test sum(ε)/M / p < 0.15
    end
    @testset "1D Henon with offset start of pred" begin
        γ=3; B=2
        em = cubic_shell_embedding(utrain,γ,1,B,1,ConstantBoundary(10.))
        ustart = U[N_train+100-4:N_train+100]
        utest_off  = U[N_train+100:N_train+p+100]
        upred = temporalprediction(utrain,em, p; initial_ts=ustart)

        @test upred[1] == ustart[end]
        @test sum(abs.(utest_off[2]-upred[2]))/M/p < 0.05
        ε = [sum(abs.(utest_off[i]-upred[i])) for i=1:p+1]
        @test sum(ε)/M / p < 0.15
    end

    @testset "LLM γ=$γ, B=$B" for γ=2:3, B=1:2
        method = LinearLocalModel(TimeseriesPrediction.ω_safe, 0.001, 1.)
        em = cubic_shell_embedding(utrain, γ, 1, B,1, ConstantBoundary(10.))
        upred = temporalprediction(utrain,em,p; method=method)

        @test upred[1] == utrain[end]
        @test sum(abs.(utest[2]-upred[2]))/M/p < 0.05
        ε = [sum(abs.(utest[i]-upred[i])) for i=1:p+1]
        @test sum(ε)/M / p < 0.15
    end
end

@testset "Coupled Henon 2D" begin
    #Size
    X=5
    Y=5
    N_train = 1000
    p = 5
    U,V = coupled_henon2D(X,Y,N_train+p,rand(X,Y), rand(X,Y))
    #Reconstruct this #
    utrain = U[1:N_train]
    vtrain = V[1:N_train]
    utest  = U[N_train:N_train+p]
    vtest  = V[N_train:N_train+p]

    @testset "γ=$γ, B=$B" for γ=2:3, B=1:2
        em = cubic_shell_embedding(utrain,γ,1,B,1,ConstantBoundary(10.))
        upred = temporalprediction(utrain,em,p)

        @test upred[1] == utrain[end]
        ε = [sum(abs.(utest[i]-upred[i])) for i=1:p+1]
        @test sum(ε)/X/Y / p < 0.15
    end
end


text:
module Simulations

    using ProgressMeter

    using ..CommonFunctions
    using ..Equations
    using ..Integrators
    using ..Solutions

    import ..Solutions: DEFAULT_NSAVE, DEFAULT_NWRITE

    export Simulation, ParallelSimulation
    export run!

    include("simulations/simulation.jl")
    include("simulations/parallel_simulation.jl")

end


text:
"""
    FockSpace <: HilbertSpace

[`HilbertSpace`](@ref) defining a Fock space for bosonic operators.
See also: [`Destroy`](@ref), [`Create`](@ref)
"""
struct FockSpace{S} <: ConcreteHilbertSpace
    name::S
end
Base.:(==)(h1::T,h2::T) where T<:FockSpace = (h1.name==h2.name)

"""
    Destroy <: QSym

Bosonic operator on a [`FockSpace`](@ref) representing the quantum harmonic
oscillator annihilation operator.
"""
struct Destroy{H<:HilbertSpace,S,A,M} <: QSym
    hilbert::H
    name::S
    aon::A
    metadata::M
    function Destroy{H,S,A,M}(hilbert::H, name::S, aon::A, metadata::M) where {H,S,A,M}
        @assert has_hilbert(FockSpace,hilbert,aon)
        new(hilbert,name,aon,metadata)
    end
end

"""
    Create <: QSym

Bosonic operator on a [`FockSpace`](@ref) representing the quantum harmonic
oscillator creation operator.
"""
struct Create{H<:HilbertSpace,S,A,M} <: QSym
    hilbert::H
    name::S
    aon::A
    metadata::M
    function Create{H,S,A,M}(hilbert::H, name::S, aon::A, metadata::M) where {H,S,A,M}
        @assert has_hilbert(FockSpace,hilbert,aon)
        new(hilbert,name,aon,metadata)
    end
end

for T in (:Create, :Destroy)
    @eval Base.isequal(a::$T, b::$T) = isequal(a.hilbert, b.hilbert) && isequal(a.name, b.name) && isequal(a.aon, b.aon)
end

for f in [:Destroy,:Create]
    @eval $(f)(hilbert::H, name::S, aon::A; metadata::M=NO_METADATA) where {H,S,A,M} = $(f){H,S,A,M}(hilbert,name,aon,metadata)
    @eval $(f)(hilbert::FockSpace, name; metadata=NO_METADATA) = $(f)(hilbert,name,1; metadata)
    @eval function $(f)(hilbert::H, name::S, aon::A; metadata::M=NO_METADATA) where {H<:ProductSpace,S,A<:Int,M}
        if hilbert.spaces[aon] isa ClusterSpace
            hilbert.spaces[get_i(aon)].op_name[] = name
            op = $(f)(hilbert.spaces[aon].original_space, name; metadata)
            return _cluster(hilbert, op, aon) #todo OK?
        else
            return $(f){H,S,A,M}(hilbert,name,aon,metadata)
        end
    end
    @eval function $(f)(hilbert::ProductSpace, name; metadata=NO_METADATA)
        i = findall(x->isa(x,FockSpace) || isa(x,ClusterSpace{<:FockSpace}),hilbert.spaces)
        if length(i)==1
            return $(f)(hilbert, name, i[1]; metadata)
        else
            isempty(i) && error("Can only create $($(f)) on FockSpace! Not included in $(hilbert)")
            length(i)>1 && error("More than one FockSpace in $(hilbert)! Specify on which Hilbert space $($(f)) should be created with $($(f))(hilbert,name,i)!")
        end
    end
    @eval function Base.hash(op::T, h::UInt) where T<:($(f))
        hash($(f), hash(op.hilbert, hash(op.name, hash(op.aon, h))))
    end
end

Base.adjoint(op::Destroy) = Create(op.hilbert,op.name,acts_on(op); op.metadata)
Base.adjoint(op::Create) = Destroy(op.hilbert,op.name,acts_on(op); op.metadata)

# Commutation relation in simplification
function *(a::Destroy,b::Create)
    check_hilbert(a,b)
    aon_a = acts_on(a)
    aon_b = acts_on(b)
    if aon_a == aon_b
        return b*a + 1
    elseif aon_a < aon_b
        return QMul(1, [a,b])
    else
        return QMul(1, [b,a])
    end
end
ismergeable(::Destroy,::Create) = true


text:
# This file contains the functions for the GP model

# load packages
using Yeppp # for fast exp function

# type(s):

"""
GPModel is the type for the Gaussian processes model with qudaric mean function
and some covariance function.
"""
type GPModel
  est_method::String
  Β::Array
  logΦ::Array # log([σ, σ_kernel lengthscales])
  cov_inv::Array
  diff::Array # y-XΒ
  x::Array # training data
  indecies_to_include::Vector
end

# functions:

doc"""
    meanfunction(x::Array, intercept::Bool, indecies_to_include)

meanfunction computes the design matrix for the quadratic mean function.
"""
function meanfunction(x::Array, intercept::Bool, indecies_to_include)

  nbr_obs = size(x,1)
  dim = size(x,2)

  n_col = dim + sum(dim:-1:1)
  X = zeros(nbr_obs,n_col)

  X[:,1:dim] = x
  X[:,dim+1:2*dim] = x.^2

  idx = 2*dim
  for i = 1:dim-1
    for j = i+1:dim
      idx = idx + 1
      X[:,idx] = x[:,i].*x[:,j]
    end
  end

  if intercept
    X = [ones(size(x,1),1) X]
    return X[:,indecies_to_include]
  else
    #indecies_to_include  = indecies_to_include[2:end]
    return X[:,indecies_to_include]
  end

end

doc"""
    predict(x_pred::Array{Float64}, gp::GPModel, noisy_pred::Bool=true)

predict computes the predictions at the values in x_pred.
"""
function predict(x_pred::Array{Float64}, gp::GPModel, noisy_pred::Bool=true)
  mu_pred = zeros(size(x_pred,2))
  sigma_pred = zeros(size(x_pred,2))
  prediction_sample = zeros(size(x_pred,2))

  for i = 1:size(x_pred,2)
    (mu_pred[i], sigma_pred[i],prediction_sample[i]) = prediction(x_pred[:,i], gp,noisy_pred)
  end
  return mu_pred, sigma_pred, prediction_sample
end

doc"""
    predict(x_pred::Array{Float64}, gp::GPModel, pred_method::String, est_method::String,noisy_pred::Bool=true)

predict computes the predictions at the values in x_pred.
"""
function predict(x_pred::Array{Float64}, gp::GPModel, pred_method::String, est_method::String,noisy_pred::Bool=true, return_var::Bool=false)

  mu_pred = zeros(size(x_pred,2))
  sigma_pred = zeros(size(x_pred,2))
  prediction_sample = zeros(size(x_pred,2))

  if est_method == "ml"

    for i = 1:size(x_pred,2)
      (mu_pred[i], sigma_pred[i],prediction_sample[i]) = prediction(x_pred[:,i], gp,noisy_pred)
    end

  else

    error("The two-stage method is not implemented.")

  end

  if !return_var
    if pred_method == "mean"
      return mu_pred
    else
      return prediction_sample
    end
  else
    if pred_method == "mean"
      return mu_pred, sigma_pred
    else
      return prediction_sample, sigma_pred
    end
  end

end

doc"""
    prediction(x_pred::Array{Float64}, gp::GPModel, noisy_pred::Bool = true)

prediction computes the prediction at x_pred.
"""
function prediction(x_pred::Array{Float64}, gp::GPModel, noisy_pred::Bool = true)

    # pre-allocate arrays
    c = zeros(size(gp.x,2))
    covtemp = zeros(size(gp.x,2))
    mu_pred = zeros(1,1)
    sigma_pred = zeros(1,1)
    temp1 = zeros(1,1)
    temp2 = copy(gp.diff)

    # set parameters
    l_vec = -1/2*1./(exp(gp.logΦ[3:end]).^2)

    # calc variance vectro, i.e. calc K(x_pred, X)

    # compute for each dimension
    for i = 1:length(l_vec)
      covtemp = l_vec[i]*(x_pred[i] - gp.x[i,:]).^2
      c = c + vec(covtemp)
    end

    # take exponent
    Yeppp.exp!(c) # still slow.... migth be possible to fix this by calling the c function for exp

    # scale matrix with σ_kernel
    scale!(c, exp(gp.logΦ[2])^2)
    #c = exp(gp.logΦ[2])^2*c

    if noisy_pred
      # this is the slow part...
      LinAlg.A_mul_B!(mu_pred,meanfunction(x_pred',true,gp.indecies_to_include),gp.Β)
      LinAlg.A_mul_B!(temp2,gp.cov_inv,gp.diff)
      LinAlg.A_mul_B!(temp1,c',temp2)
      LinAlg.axpy!(1.,temp1,mu_pred)
      #mu_pred = meanfunction(x_pred',true,gp.indecies_to_include)*gp.Β + c'*gp.cov_inv*gp.diff
      LinAlg.A_mul_B!(temp2,gp.cov_inv,c)
      LinAlg.A_mul_B!(sigma_pred,c',temp2)
      sigma_pred = exp(gp.logΦ[1])^2 + exp(gp.logΦ[2])^2 + sigma_pred
      #LinAlg.axpy!(1.,exp(gp.logΦ[1])^2,sigma_pred)
      #LinAlg.axpy!(1.,exp(gp.logΦ[2])^2,sigma_pred)
      #sigma_pred = exp(gp.logΦ[2])^2 + exp(gp.logΦ[1])^2 - c'*gp.cov_inv*c #+ exp(gp.logΦ[1])^2
    else
      # this is the slow part...
      LinAlg.A_mul_B!(mu_pred,meanfunction(x_pred',true,gp.indecies_to_include),gp.Β)
      LinAlg.A_mul_B!(temp2,gp.cov_inv,gp.diff)
      LinAlg.A_mul_B!(temp1,c',temp2)
      LinAlg.axpy!(1.,temp1,mu_pred)
      #mu_pred = meanfunction(x_pred',true,gp.indecies_to_include)*gp.Β + c'*gp.cov_inv*gp.diff
      LinAlg.A_mul_B!(temp2,gp.cov_inv,c)
      LinAlg.A_mul_B!(temp1,c',temp2)
      LinAlg.axpy!(1.,temp1,sigma_pred)
      sigma_pred =  exp(gp.logΦ[2])^2 + sigma_pred

      #LinAlg.axpy!(1.,exp(gp.logΦ[2])^2,sigma_pred)

      #mu_pred = meanfunction(x_pred',true,gp.indecies_to_include)*gp.Β + c'*gp.cov_inv*gp.diff
      #sigma_pred = exp(gp.logΦ[2])^2 - c'*gp.cov_inv*c #+ exp(gp.logΦ[1])^2
    end
    if sigma_pred[1] < 0
      return mu_pred[1], sigma_pred[1], NaN #rand(Normal(mu_pred[1],sqrt(sigma_pred[1])))
    else
      return mu_pred[1], sigma_pred[1], rand(Normal(mu_pred[1],sqrt(sigma_pred[1])))
    end

  #=
  # pre-allocate arrays
  c = zeros(size(gp.x,2))
  covtemp = zeros(size(gp.x,2))

  # set parameters
  l_vec = -1/2*1./(exp(gp.logΦ[3:end]).^2)

  # calc variance vectro, i.e. calc K(x_pred, X)

  # compute for each dimension
  for i = 1:length(l_vec)
    covtemp = l_vec[i]*(x_pred[i] - gp.x[i,:]).^2
    # problem v0.4 !!!!!!!!!
    # was: c = c + covtempc
    c = c + vec(covtemp)
  end

  # take exponent
  Yeppp.exp!(c) # still slow.... migth be possible to fix this by calling the c function for exp

  # scale matrix with σ_kernel
  c = exp(gp.logΦ[2])^2*c

  # compute prediction
  res = zeros(1)
  res2 = zeros(1)
  d = At_mul_B(c,gp.cov_inv)
  A_mul_B!(res,d,gp.diff)
  A_mul_B!(res2,d,c)

  if noisy_pred
    mu_pred = meanfunction(x_pred',true,gp.indecies_to_include)*gp.Β + res #c'*gp.cov_inv*gp.diff
    sigma_pred = exp(gp.logΦ[2])^2 + exp(gp.logΦ[1])^2 - res2 # c'*gp.cov_inv*c #+ exp(gp.logΦ[1])^2
  else
    mu_pred = meanfunction(x_pred',true,gp.indecies_to_include)*gp.Β + res #c'*gp.cov_inv*gp.diff
    sigma_pred = exp(gp.logΦ[2])^2 - res2 #c'*gp.cov_inv*c #+ exp(gp.logΦ[1])^2
  end
  if sigma_pred[1] < 0
    return mu_pred[1], sigma_pred[1], NaN #rand(Normal(mu_pred[1],sqrt(sigma_pred[1])))
  else
    return mu_pred[1], sigma_pred[1], rand(Normal(mu_pred[1],sqrt(sigma_pred[1])))
  end

  =#

end

doc"""
    covariancefunction(x::Array, Φ::Array, kernel::String)

covariancefunction computes the covariance function using the SE kernel.
"""
function covariancefunction(x::Array, Φ::Array, kernel::String)

  # set parameters
  σ2 = Φ[1]^2
  σ_kernel2  =  Φ[2]^2
  l_vec = -1/2*(1./(Φ[3:end].^2))

  # pre-allocate matricies
  covtemp = zeros(length(x[1,:]),length(x[1,:]))
  cov_m = zeros(length(x[1,:]),length(x[1,:]))
  sigma_m = eye(size(cov_m,1))

  # calc sigma_m
  scale!(sigma_m,σ2)

  # set power2(x) = x^2 function
  power2(x) = x^2

  # compute for each dimension
  for i = 1:length(l_vec)
    broadcast!(-, covtemp, x[i,:], x[i,:]')
    broadcast!(power2, covtemp, covtemp)
    scale!(covtemp,l_vec[i])
    broadcast!(+, cov_m,cov_m,covtemp)
  end

  # take exponent
  Yeppp.exp!(cov_m) # still slow.... migth be possible to fix this by calling the c function for exp

  # scale matrix with σ_kernel
  scale!(cov_m,σ_kernel2)

  # add sigma_m on the main diagonal
  broadcast!(+, cov_m,cov_m,sigma_m)

  return cov_m

end

doc"""
    leastsquaresestimation(X::Matrix{Float64}, y::Vector{Float64})

Comutes regression estiamtins of beta for X*beta = y, also computes the confidence
bands for the beta's.
"""
function leastsquaresestimation(X::Matrix{Float64}, y::Vector{Float64})

  # compute betas
  beta = \(X'*X,X')*y

  # compute hat matrix
  H = X*\(X'*X,X')

  # compute mse
  sse = y'*(eye(size(H,1),size(H,2))-H)*y
  mse = sse/(size(X,1)-size(X,2))

  # compute variance matrix
  s = mse.*\(X'*X,eye(size(X'*X,1),size(X'*X,1)))

  conf_limits = 1.96*diag(s) # approximate t(1-0.05/2,n-p) with 1.96

  return beta, conf_limits

end

# not storing calc

doc"""
    calculate_common!(theta::Array, buffer::Array{Float64}, x::Array{Float64})

Computes the common calculations for the objectiv function, i.e. the negative ml function
and the gradient of the objective function.
"""
function calculate_common!(theta::Array, buffer::Array{Float64}, x::Array{Float64})

  # set parameters
  σ2 = exp(theta[end-size(x,1)-1])^2
  σ_kernel2  =  exp(theta[end-size(x,1)])^2
  l_vec = -1/2./(exp(theta[end-size(x,1)+1:end]).^2)

  # pre-allocate matricies
  sigma_m = eye(length(x[1,:]),length(x[1,:]))
  cov_m = zeros(length(x[1,:]),length(x[1,:]))
  difftemp = zeros(size(cov_m))
  diff_m = zeros(size(cov_m))
  cov_m = zeros(size(cov_m))

  # calc sigma_m
  scale!(sigma_m,σ2)

  # set power2(x) = x^2 function
  power2(x) = x^2

  # compute for each dimension
  for i = 1:length(l_vec)
    broadcast!(-, difftemp, x[i,:], x[i,:]')
    broadcast!(power2, difftemp, difftemp)
    scale!(difftemp,l_vec[i])
    broadcast!(+, diff_m,diff_m,difftemp)
  end

  # take exponent
  Yeppp.exp!(diff_m) # still slow.... migth be possible to fix this by calling the c function for exp

  # set cov matrix
  cov_m = diff_m

  # scale matrix with σ_kernel
  scale!(cov_m,σ_kernel2)

  broadcast!(+, cov_m,cov_m,sigma_m)
  # v0.4 problem
  # issymmetric(cov_m) is not implemented in v4.0
  if cov_m == cov_m'  && isposdef(cov_m) # check that cov matrix is symmetric and positive semi def.
    cov_m_inv = inv(cov_m)
    buffer[:,:] = [diff_m; cov_m; cov_m_inv]
  else
    # set buffer to something here
    cov_m_inv = zeros(size(cov_m))
    buffer[:,:] = [diff_m; cov_m; cov_m_inv]
  end

end

doc"""
    objectivfunction(x::Array{Float64},X::Array{Float64},y::Vector{Float64},
  kernel::String, theta::Array, buffer::Array{Float64})

The negative maximum likelihood function for the GP model.
"""
function  objectivfunction(x::Array{Float64},X::Array{Float64},y::Vector{Float64},
  kernel::String, theta::Array, buffer::Array{Float64})

  # calc common parts of the objective function and gradient
  calculate_common!(theta,buffer,x) # buffer = [diff_m; cov_m; cov_m_inv]

  # set parameters
  beta = theta[1:size(X,2)] # theta = [beta Phi]

  # set matrices that are stored on buffre
  diff_m = buffer[1:length(x[1,:]),:]
  cov_m = buffer[length(x[1,:])+1:2*length(x[1,:]),:]
  cov_m_inv = buffer[2*length(x[1,:])+1:3*length(x[1,:]),:]

  # calc chol factorizations
  if sum(cov_m_inv) != 0 # if the cov matrix is symmetrix and pos semi def compute negloglik

    R = chol(cov_m)
    # calc negloglik
    negloglik =  2*sum(log(diag(R))) + (y-X*beta)'*cov_m_inv*(y-X*beta)
    negloglik = negloglik[1]

  else  # set negloglik to some value

    println("non-symmetric covariance matrix at logΦ:")
    println(theta[end-4:end])
    negloglik = Inf

  end

  return  negloglik

end

doc"""
    g!(y::Array{Float64},x::Array{Float64},X::Array{Float64}, stor::Array,
  buffer::Array{Float64}, theta::Array)

Computes the gradient of the negative log-likelihood function.
"""
function g!(y::Array{Float64},x::Array{Float64},X::Array{Float64}, stor::Array,
  buffer::Array{Float64}, theta::Array)

  # calc common parts of the objective function and gradient
  calculate_common!(theta,buffer,x)

  # set parameters
  σ2 = exp(theta[end-size(x,1)-1])^2
  σ_kernel2  =  exp(theta[end-size(x,1)])^2
  l_vec = 1./(exp(theta[end-size(x,1)+1:end]).^3)
  beta = theta[1:size(X,2)]

  # pre-allocate matricies
  dCov_b_d_Phi_i_val = zeros(1)
  dCov_b_d_Phi_i = zeros(length(x[1,:]),length(x[1,:]))
  dCov_b_d_Phi_i_inner_deriv = zeros(length(x[1,:]),length(x[1,:]))

  # set matrices that are stored on buffre
  diff_m = buffer[1:length(x[1,:]),:]
  cov_m = buffer[length(x[1,:])+1:2*length(x[1,:]),:]
  cov_m_inv = buffer[2*length(x[1,:])+1:3*length(x[1,:]),:]

  power2(x) = x^2

  # calc chol factorizations
  if sum(cov_m_inv) != 0

    # calc gradient for beta
    stor[1:length(beta)] = -2*X'*cov_m_inv*(y-X*beta)

    # gradient for Phi
    for i = length(beta)+1:length(stor)

      if i == length(beta)+1
        # dCov_b_d_\sigma
        dCov_b_d_Phi_i = 2*sqrt(σ2)*eye(length(x[1,:]),length(x[1,:]))
      elseif i == length(beta)+2
        # dCov_b_d_\sigma_kernel
        dCov_b_d_Phi_i = 2*sqrt(σ_kernel2)*diff_m

      else
        # dCov_b_d_l_i
        idx = i - (length(beta)+2)
        dCov_b_d_Phi_i = scale(σ_kernel2, diff_m)
        broadcast!(-, dCov_b_d_Phi_i_inner_deriv, x[idx,:], x[idx,:]')
        broadcast!(power2, dCov_b_d_Phi_i_inner_deriv, dCov_b_d_Phi_i_inner_deriv)
        scale!(dCov_b_d_Phi_i_inner_deriv,l_vec[idx])
        dCov_b_d_Phi_i = broadcast!(*,dCov_b_d_Phi_i,dCov_b_d_Phi_i,dCov_b_d_Phi_i_inner_deriv)

      end

      trace_Phi_i = trace(cov_m_inv*dCov_b_d_Phi_i)
      dCov_b_d_Phi_i_val  = -(y-X*beta)'*cov_m_inv*dCov_b_d_Phi_i*cov_m_inv*(y-X*beta) + trace_Phi_i
      stor[i] = dCov_b_d_Phi_i_val[1]
    end

  else

    # set gradient to something
    stor[:] = 1000
    println("Non-symmetric covariance matrix at logΦ:")
    show(theta[end-size(x,1)-1:end])
    negloglik = Inf

  end

end

doc"""
    ml_est(gp::GPModel, data_training::Array{Float64}, kernel::String="SE",lasso::Bool=true)

Computes the maximum likelihood estimation of the parameters of the GP model,
by minimizing the negative log-likelihood function.
"""
function ml_est(gp::GPModel, data_training::Array{Float64}, kernel::String="SE",lasso::Bool=true, percentage::Float64=0.05, sided::String="both")

  @printf "Starting ml estimation"

  # v0.4 problem
  y = vec(data_training[end,:]) # need this to fix issue with types..
  x = data_training[1:end-1,:]
  (y,x) = removeoutlier(y,x,Int64(length(y)*percentage),sided)

  @printf "Length of training data: %d" size(x,2)
  dim = size(x,1)
  n_col = dim + sum(dim:-1:1)

  X = meanfunction(x',true,1:n_col+1)
  buffer = zeros(3*size(X,1), size(X,1))
  indeices_keep = collect(1:n_col+1)
  beta_0 = zeros(n_col+1)
  logΦ_0 = log([std(y); std(y)/sqrt(2); std(x,2)])
  theta_0 = zeros(n_col+1+length(logΦ_0))

  if !lasso
    beta_0 = inv(X'*X)*X'*y # start value for the the beta:s should be the regression est of the beta:s
  else
    X = meanfunction(x',false,collect(1:n_col))
    lassopath = Lasso.fit(LassoPath,X,y, λminratio = 0.02,nλ = 500,intercept=true)
    ind_best = indmax(lassopath.pct_dev)
    beta_0[1] = lassopath.b0[ind_best]

    for i = 2:length(beta_0)
      beta_0[i] = lassopath.coefs[i-1,ind_best]
    end

    indeices_keep = find(x -> x != 0.00,beta_0)
    beta_0 = beta_0[indeices_keep]
    X = [ones(size(X,1)) X] # add intercept
    X = X[:,indeices_keep]
  end

  # set start values
  theta_0 = [beta_0; logΦ_0]

  println("Start values for optimization:")
  show(theta_0)

  # set closures for objectivfunction and g!
  function objectivfunctionclosure(theta)
    return objectivfunction(x,X,y,kernel, theta, buffer)
  end

  function gclosure!(theta,stor)
    return g!(y,x,X, stor,buffer, theta)
  end

  # settings for the optimization algorithm
  options = Optim.Options(show_trace=true,iterations=1000,time_limit=60*10)

  # optimize parameters using conjugate gradient descent
  opt_res = Optim.optimize(objectivfunctionclosure,
                          gclosure!,
                          theta_0,
                          ConjugateGradient(),
                          options)

  # parameters found
  theta_hat = Optim.minimizer(opt_res)

  # set GPModel
  gp.est_method = "ml"
  gp.Β = theta_hat[1:length(beta_0)]
  gp.logΦ = theta_hat[end-length(logΦ_0)+1:end]
  gp.cov_inv = inv(covariancefunction(x,exp(theta_hat[end-length(logΦ_0)+1:end]),kernel))
  gp.diff = y-X*theta_hat[1:length(beta_0)]
  gp.x = x
  gp.indecies_to_include = indeices_keep

  println("Estimated parameters of the GP model:")
  println("Estimated Β:s:")
  println(theta_hat[1:length(beta_0)])
  println("Estimated log Φ:s:")
  println(theta_hat[end-length(logΦ_0)+1:end])

end

doc"""
    removeoutlier(y::Vector{Float64}, x::Array{Float64}, nbr_outliers::Int64=10, sided::String="both")

removeoutlier removes nbr_outliers outliers either from both tails or from left/rigth tail.
"""
function removeoutlier(y::Vector{Float64}, x::Array{Float64}, nbr_outliers::Int64=10, sided::String="both")

  idx_rm = Vector{Int64}
  idx = sortperm(y, rev=true)

  if sided == "both"
    # remove nbr_outliers on both tails
    idx_rm = [idx[1:nbr_outliers]; idx[end-(nbr_outliers-1):end]]
  elseif sided == "left"
    # remove nbr_outliers on left tail
    idx_rm = idx[end-(nbr_outliers-1):end]
  else
    # remove nbr_outliers on rigth tail
    idx_rm = idx[1:nbr_outliers]
  end

  return y[setdiff(1:end,idx_rm)], x[:,setdiff(1:end,idx_rm)]

end

doc"""
    RMSE(targets::Vector{Float64}, predictions::Vector{Float64})

Computes the root mean square error.

"""
function RMSE(targets::Vector{Float64}, predictions::Vector{Float64})
  return sqrt(mean((targets-predictions).^2))
end

max_stars_repo_name: SamuelWiqvist/adamcmcpaper

max_issues_repo_name: SamuelWiqvist/adamcmcpaper

max_forks_repo_name: SamuelWiqvist/adamcmcpaper


text:
## preallocation

function _preallocate_spmatrix(::Type{T}, sizes::Tuple, sizehint) where {T}
    res = spzeros(T, sizes...)
    @static if VERSION >= v"1.7"
        sizehint!(res, first(sizes) * sizehint)
    end

    return res
end

__an_elt(χ::Character) = first(first(conjugacy_classes(χ)))
__an_elt(α::AlgebraElement) = first(basis(parent(α)))

_projection_size(p::PermutationGroups.AbstractPerm) = (d = degree(p); (d, d))
_projection_size(m::AbstractMatrix) = size(m)

_projection_size(::Nothing, χ::Character) = _projection_size(__an_elt(χ))
_projection_size(::Nothing, α::AlgebraElement) = _projection_size(__an_elt(α))
_projection_size(hom::InducedActionHomomorphism, χ::Character) =
    _projection_size(induce(hom, __an_elt(χ)))
_projection_size(hom::InducedActionHomomorphism, α::AlgebraElement) =
    _projection_size(induce(hom, __an_elt(α)))

_hint(χ::Character) = length(conjugacy_classes(χ))
_hint(α::AlgebraElement) = count(!iszero, StarAlgebras.coeffs(α))

preallocate(::Type{T}, χ::Union{Character,AlgebraElement}) where {T} =
    preallocate(T, nothing, χ)

function preallocate(
    hom::InducedActionHomomorphism,
    χ::Union{Character,AlgebraElement},
)
    T = Base._return_type(*, (coeff_type(hom), eltype(χ)))
    return preallocate(T, hom, χ)
end

function preallocate(
    ::Type{T},
    a::Union{Nothing,InducedActionHomomorphism},
    χ::Union{Character,AlgebraElement},
) where {T}
    sizes = _projection_size(a, χ)
    return _preallocate_spmatrix(T, sizes, _hint(χ))
end

## matrix projection [irreducible]
"""
    matrix_projection([hom::InducedActionHomomorphism, ]χ::Character{T})
Compute matrix projection associated to character `χ`.

Returned `M<:AbstractMatrix{T}` of size `(d, d)` where the degree `d` of the
projecion is determined by elements in `conjugacy_classes(χ)`. E.g. `d` could
be equal to the `degree` when conjugacy classes consist of `AbstractPerms`.
If the homomorphism is passed, the dimension will be derived in similar
manner from the elements of the image of the homomorphism.

The precise type of `M` can be altered by overloading

```
preallocate(::Type{T}, [hom::InducedActionHomomorphism, ]χ::Character)
```
"""
matrix_projection(χ::Character{T}) where {T} =
    _mproj_outsT!(preallocate(T, χ), χ)

function matrix_projection(
    hom::InducedActionHomomorphism,
    χ::Character{T},
) where {T}
    return _mproj_outsT!(preallocate(hom, χ), hom, χ)
end

function matrix_projection(χ::Character{T}) where {T<:Rational}
    if all(isone ∘ Cyclotomics.conductor, table(χ))
        return _mproj_fitsT!(preallocate(T, χ), χ)
    else
        return _mproj_outsT!(preallocate(T, χ), χ)
    end
end

matrix_projection(χ::Character{T}) where {T<:Union{Cyclotomic,Complex}} =
    _mproj_fitsT!(preallocate(T, χ), χ)
function matrix_projection(
    hom::InducedActionHomomorphism,
    χ::Character{T},
) where {T<:Union{Cyclotomic,Complex}}
    return _mproj_fitsT!(preallocate(hom, χ), hom, χ)
end

function matrix_projection(χ::Character{T}) where {T<:AbstractFloat}
    if all(isreal, table(χ))
        return _mproj_fitsT!(preallocate(T, χ), χ)
    else
        return _mproj_outsT!(preallocate(T, χ), χ)
    end
end

function matrix_projection(
    hom::InducedActionHomomorphism,
    χ::Character{T},
) where {T<:AbstractFloat}
    if all(isreal, table(χ))
        return _mproj_fitsT!(preallocate(hom, χ), hom, χ)
    else
        return _mproj_outsT!(preallocate(hom, χ), hom, χ)
    end
end

_mproj_fitsT!(args...) = _mproj_outsT!(args...)

function _mproj_outsT!(mproj::AbstractMatrix{T}, χ::Character) where {T}
    mproj .= sum(
        c .* matrix_projection_irr(Character(table(χ), i)) for
        (i, c) in pairs(constituents(χ)) if !iszero(c)
    )
    return mproj
end

function _mproj_outsT!(
    mproj::AbstractMatrix{T},
    hom::InducedActionHomomorphism,
    χ::Character,
) where {T}
    mproj .= sum(
        c .* matrix_projection_irr(hom, Character(table(χ), i)) for
        (i, c) in pairs(constituents(χ)) if !iszero(c)
    )
    return mproj
end

"""
    matrix_projection_irr([::Type, ]χ::Character)
    matrix_projection_irr(hom::InducedActionHomomorphism, χ::Character)
Compute matrix projection associated to an *irreducible* character `χ`.

The returned matrix defines so called *isotypical* projection.

See also [matrix_projection](@ref).
"""
matrix_projection_irr(χ::Character{T}) where {T} = matrix_projection_irr(T, χ)
matrix_projection_irr(::Type{T}, χ::Character) where {T} =
    matrix_projection_irr_acc!(preallocate(T, χ), χ, 1)

matrix_projection_irr(hom::InducedActionHomomorphism, χ::Character) =
    matrix_projection_irr_acc!(preallocate(hom, χ), hom, χ, 1)

function matrix_projection_irr_acc!(
    result::AbstractMatrix,
    χ::Character,
    weight,
)
    @assert isirreducible(χ)
    LinearAlgebra.checksquare(result)
    ccls = conjugacy_classes(χ)
    vals = collect(values(χ))
    weight *= degree(χ) // sum(length, ccls)
    matrix_projection_irr_acc!(result, vals, ccls, weight)
    return result
end

function matrix_projection_irr_acc!(
    result::AbstractMatrix,
    vals,
    ccls::AbstractVector{<:AbstractOrbit{<:PermutationGroups.AbstractPerm}},
    weight,
) where {T}
    iszero(weight) && return result
    for (val, cc) in zip(vals, ccls)
        iszero(val) && continue
        w = weight * val
        for g in cc
            for i in 1:size(result, 1)
                result[i, i^g] += w
            end
        end
    end

    return result
end

function matrix_projection_irr_acc!(
    result::AbstractMatrix,
    vals,
    ccls::AbstractVector{<:AbstractOrbit{<:AbstractMatrix}},
    weight,
)
    # TODO: call to inv(Matrix(g)) is a dirty hack, since if `g`
    # is given by a sparse matrix `inv(g)` will fail.
    for (val, cc) in zip(vals, ccls)
        iszero(val) && continue
        for g in cc
            res .+= (weight * val) .* inv(convert(Matrix, g))
        end
    end
    return res
end

function matrix_projection_irr_acc!(
    result::AbstractMatrix,
    hom::InducedActionHomomorphism,
    χ::Character,
    weight,
)
    @assert isirreducible(χ)
    LinearAlgebra.checksquare(result)
    vals = collect(values(χ))
    ccls = conjugacy_classes(χ)
    weight *= degree(χ) // sum(length, ccls)
    matrix_projection_irr_acc!(result, hom, vals, ccls, weight)
    return result
end

function matrix_projection_irr_acc!(
    result::AbstractMatrix,
    hom::InducedActionHomomorphism{<:ByPermutations},
    class_values,
    conjugacy_cls,
    weight,
)
    iszero(weight) && return result
    for (val, ccl) in zip(class_values, conjugacy_cls)
        iszero(val) && continue
        w = weight * val
        for g in ccl
            h = induce(hom, g)
            @assert h isa PermutationGroups.Perm
            for i in 1:size(result, 1)
                result[i, i^h] += w
            end
        end
    end
    return result
end

function matrix_projection_irr_acc!(
    result::AbstractMatrix,
    hom::InducedActionHomomorphism{<:ByLinearTransformation},
    class_values,
    conjugacy_cls,
    weight,
)
    iszero(weight) && return result
    for (val, cc) in zip(class_values, conjugacy_cls)
        iszero(val) && continue
        w = weight * val
        for g in cc
            result .+= w .* induce(hom, inv(g))
        end
    end
    return result
end

"""
    matrix_representation([::Type{T}, ]α::AlgebraElement)
    matrix_representation([::Type{T}, ]hom::InducedActionHomomorphism, α::AlgebraElement)
Compute matrix representative of the given algebra element `α`.

If `α` is an element of the group (or monoid) algebra `RG` of `G`, then every
linear representation `π : G → GL(V)` gives rise to an algebra homomorphism
`π : RG → B(V)` (bounded operators on `V`). This function evaluates this
representation given either by
 * the natural action of `G` on some `V` (e.g. every permutation group of degree
 `d` acts on `d`-dimensional `V` by basis vector permutation), or
 * by action homomorphism `hom`.

See also [matrix_projection](@ref).
"""
matrix_representation(α::AlgebraElement) = matrix_representation(eltype(α), α)
matrix_representation(::Type{T}, α::AlgebraElement) where {T} =
    matrix_representation_acc!(preallocate(T, α), α)

matrix_representation(hom::InducedActionHomomorphism, α::AlgebraElement) =
    matrix_representation_acc!(preallocate(hom, α), hom, α)

function matrix_representation_acc!(
    result::AbstractMatrix,
    α::AlgebraElement{
        <:StarAlgebra{<:PermutationGroups.AbstractPermutationGroup},
    },
)
    b = basis(parent(α))
    for (idx, val) in StarAlgebras._nzpairs(StarAlgebras.coeffs(α))
        g = b[idx]
        iszero(val) && continue
        for i in 1:size(result, 1)
            result[i, i^g] += val
        end
    end

    return result
end

function matrix_representation_acc!(
    result::AbstractMatrix,
    hom::InducedActionHomomorphism{<:ByPermutations},
    α::AlgebraElement,
)
    b = basis(parent(α))
    for (idx, val) in StarAlgebras._nzpairs(StarAlgebras.coeffs(α))
        iszero(val) && continue
        g = induce(hom, b[idx])
        @assert g isa PermutationGroups.Perm
        for i in 1:size(result, 1)
            result[i, i^g] += val
        end
    end

    return result
end

function matrix_representation_acc!(
    result::AbstractMatrix,
    hom::InducedActionHomomorphism{<:ByLinearTransformation},
    α::AlgebraElement,
)
    b = basis(parent(α))
    for (idx, val) in StarAlgebras._nzpairs(StarAlgebras.coeffs(α))
        iszero(val) && continue
        result .+= val .* induce(hom, inv(b[idx]))
    end

    return result
end

## Finding basis of the row-space (right image) of an AbstractMatrix

"""
    image_basis(A::AbstractMatrix)
    image_basis([hom::InducedActionHomomorphism, ]χ::Character)
    image_basis([hom::InducedActionHomomorphism, ]α::AlgebraElement)
Return basis of the row-space of a matrix.

For characters or algebra elements return a basis of the row-space of
`matrix_projection([hom, ]χ)` or `matrix_representation([hom, ]α)`.

Two methods are employed to achieve the goal.
* By default (symbolic, exact) row echelon form is produced, and therefore there is no guarantee on the orthogonality of the returned basis vectors (rows).
* If `eltype(A) <: AbstractFloat` this is followed by a call to `svd` and the appropriate rows of its `Vt` are returned, thus the basis is (numerically) orthonormal.

# Examples:
```julia
julia> a = rand(-3:3, 3,3).//rand(1:4, 3,3);

julia> a[3, :] .= 2a[1, :] .- 1a[2, :]; a
3×3 Matrix{Rational{Int64}}:
  3//2   0//1  1//1
 -1//2  -2//1  1//2
  7//2   2//1  3//2

julia> ib = SymbolicWedderburn.image_basis(a)
2×3 Matrix{Rational{Int64}}:
 1//1  0//1   2//3
 0//1  1//1  -5//12

julia> ibf = SymbolicWedderburn.image_basis(float.(a))
2×3 Matrix{Float64}:
 -0.666651  0.416657  -0.618041
  0.529999  0.847998  -8.85356e-17

```
"""
image_basis(A::AbstractMatrix) =
    ((m, p) = image_basis!(deepcopy(A)); m[1:length(p), :])

function image_basis(χ::Character)
    mpr = isirreducible(χ) ? matrix_projection_irr(χ) : matrix_projection(χ)
    image, pivots = image_basis!(mpr)
    return image[1:length(pivots), :]
end

function image_basis(hom::InducedActionHomomorphism, χ::Character)
    mpr =
        isirreducible(χ) ? matrix_projection_irr(hom, χ) :
        matrix_projection(hom, χ)
    image, pivots = image_basis!(mpr)
    return image[1:length(pivots), :]
end

function image_basis(α::AlgebraElement)
    image, pivots = image_basis!(matrix_representation(α))
    return image[1:length(pivots), :]
end

function image_basis(hom::InducedActionHomomorphism, α::AlgebraElement)
    image, pivots = image_basis!(matrix_representation(hom, α))
    return image[1:length(pivots), :]
end

image_basis!(A::AbstractMatrix) = row_echelon_form!(A)
function image_basis!(A::AbstractSparseMatrix)
    A, p = row_echelon_form!(A)
    dropzeros!(A)
    return A, p
end

function image_basis!(A::AbstractMatrix{T}) where {T<:Union{AbstractFloat,Complex}}
    A, p = row_echelon_form!(A)
    fact = svd!(convert(Matrix{T}, @view(A[1:length(p), :])))
    m = T <: Complex ? 2eps(real(T)) : eps(T)
    A_rank = sum(fact.S .> maximum(size(A)) * m)
    @assert A_rank == length(p)
    return fact.Vt, 1:length(p)
end