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

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

include("tools.jl")

@testset "Widgets" begin
    ## label
    l = label("Hello")
    @test observable(l) == l.observable
    @test observable(observable(l)) == l.observable
    @test get_gtk_property(l, "label", String) == "Hello"
    l[] = "world"
    @test get_gtk_property(l, "label", String) == "world"
    @test string(l) == string("Gtk.GtkLabelLeaf with ", string(observable(l)))
    # Test other elements of the Observables API
    counter = Ref(0)
    ofunc = on(l) do _
        counter[] += 1
    end
    @test counter[] == 0
    l[] = "changeme"
    @test counter[] == 1
    off(ofunc)
    l[] = "and again"
    @test counter[] == 1
    ldouble = map(l) do str
        str*str
    end
    @test ldouble[] == "and againand again"

    ## checkbox
    w = Window("Checkbox")
    check = checkbox(label="click me")
    push!(w, check)
    Gtk.showall(w)
    @test check[] == false
    @test Gtk.G_.active(check.widget) == false
    check[] = true
    @test 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 tgl[] == false
    @test Gtk.G_.active(tgl.widget) == false
    tgl[] = true
    @test 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"
    txt[] = "ok"
    @test get_gtk_property(txt, "text", String) == "ok"
    set_gtk_property!(txt, "text", "other direction")
    signal_emit(widget(txt), "activate", Nothing)
    @test txt[] == "other direction"
    @test get_gtk_property(num, "text", String) == "5"
    @test_throws ArgumentError num[] = 11
    num[] = 8
    @test get_gtk_property(num, "text", String) == "8"
    meld = map(txt, num) do t, n
        join((t, n), 'X')
    end
    @test meld[] == "other directionX8"
    num[] = 4
    @test meld[] == "other directionX4"
    txt[] = "4"
    @test meld[] == "4X4"
    destroy(win)

    ## textarea (aka TextView)
    v = textarea("Type something longer")
    win = Window(v)
    Gtk.showall(win)
    @test v[] == "Type something longer"
    v[] = "ok"
    @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 s[] == 8
    s[] = 3
    @test s[] == 3
    s3 = slider(IdentityRange(-3:3))
    sleep(0.01)
    @test s3[] == 0
    s3[] = -3
    @test s3[] == -3

    # Use a single observable for two widgets
    s2 = slider(1:15, observable=observable(s), orientation='v')
    @test s2[] == 3
    s2[] = 11
    @test s[] == 11
    destroy(s2)
    destroy(s)

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

    ## dropdown
    dd = dropdown(("Strawberry", "Vanilla", "Chocolate"))
    @test dd[] == "Strawberry"
    dd[] = "Chocolate"
    @test get_gtk_property(dd, "active", Int) == 2
    destroy(dd.widget)

    r = Ref(0)
    dd = dropdown(["Five"=>x->x[]=5,
                   "Seven"=>x->x[]=7])
    on(dd.mappedsignal) do f
        f(r)
    end
    dd[] = dd[]  # fire the callback for the initial value
    @test dd[] == "Five"
    @test r[] == 5
    dd[] = "Seven"
    @test dd[] == "Seven"
    @test r[] == 7
    dd[] = "Five"
    @test r[] == 5
    destroy(dd.widget)

    ## spinbutton
    s = spinbutton(1:15)
    sleep(0.01)    # For the Gtk eventloop
    @test s[] == 1
    s[] = 3
    @test 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 s[] == 1
    s[] = 1:7, 5
    @test s[] == 5

    ## cyclicspinbutton
    a = spinbutton(1:10, value = 5)
    carry_up = Observable(false)
    on(carry_up) do up
        a[] = a[] - (-1)^up
    end
    b = cyclicspinbutton(1:3, carry_up)
    @test a[] == 5
    @test b[] == 1
    b[] = 2
    @test a[] == 5
    @test b[] == 2
    b[] = 0
    @test a[] == 4
    @test b[] == 3
    b[] = 4
    @test a[] == 5
    @test 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 = Observable(t)
    tw = timewidget(t, observable=s)
    @test tw[] == s[] == t
    t = Dates.Time(2,2,2)
    tw[] = t
    @test tw[] == s[] == t
    t = Dates.Time(3,3,3)
    s[] = t
    @test tw[] == s[] == t

    # datetimewidget
    t = DateTime(1,1,1,1,1,1)
    s = Observable(t)
    tw = datetimewidget(t, observable=s)
    @test tw[] == s[] == t
    t = DateTime(2,2,2,2,2,2)
    tw[] = t
    @test tw[] == s[] == t
    t = DateTime(3,3,3,3,3,3)
    s[] = t
    @test tw[] == s[] == t

    # progressbar
    pb = progressbar(1..10)
    @test pb[] == 1
    pb[] = 5
    @test pb[] == 5
    pb = progressbar(2:8)
    @test pb[] == 2

end

const counter = Ref(0)

@testset "Button" begin
    ## button
    w = Window("Widgets")
    b = button("Click me")
    push!(w, b)
    action = map(b) do val
        counter[] += 1
    end
    Gtk.showall(w)
    cc = counter[]  # map seems to fire it once, so record the "new" initial value
    click(b::GtkObservables.Button) = ccall((:gtk_button_clicked,Gtk.libgtk),Cvoid,(Ptr{Gtk.GObject},),b.widget)
    GC.gc(true)
    click(b)
    if VERSION >= v"1.2.0"
        @test counter[] == cc+1
    else
        @test_broken counter[] == cc+1
    end
    destroy(w)

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

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 = Observable(1)
        p = player(s, 1:8)
        win = Window("Compound", 400, 100) |> (g = Grid())
        g[1,1] = p
        Gtk.showall(win)
        btn_fwd = p.widget.step_forward
        @test s[] == 1
        btn_fwd[] = nothing
        @test s[] == 2
        p.widget.play_forward[] = nothing
        for i = 1:7
                    sleep(0.1)
        end
        @test s[] == 8
        destroy(win)

        p = player(1:1000)
        win = Window("Compound 2", 400, 100)
        push!(win, frame(p))
        Gtk.showall(win)
        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)
    reveal(c, true)
    sleep(0.3)
    can_test_width = !(VERSION.minor < 3 && Sys.iswindows())
    can_test_width && @test Graphics.width(c) == 208
    @test Graphics.height(c) == 207
    @test isa(c, GtkObservables.Canvas{DeviceUnit})
    destroy(win)
    c = canvas(UserUnit, 208, 207)
    win = Window(c)
    Gtk.showall(win)
    reveal(c, true)
    sleep(1.0)
    @test isa(c, GtkObservables.Canvas{UserUnit})
    corner_dev = (DeviceUnit(208), DeviceUnit(207))
    can_test_coords = (VERSION < v"1.3" || get(ENV, "CI", nothing) != "true" || !Sys.islinux()) &&
                      can_test_width
    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)
        if can_test_coords
            # FIXME: the new JLL-based version fails on Travis.
            # Unfortunately this is difficult to debug because it doesn't replicate
            # locally or on a local headless server. See #91.
            @test GtkObservables.convertunits(DeviceUnit, c, corner_dev...) == corner_dev
            @test GtkObservables.convertunits(DeviceUnit, c, corner_usr...) == corner_dev
            @test GtkObservables.convertunits(UserUnit, c, corner_dev...) == corner_usr
            @test GtkObservables.convertunits(UserUnit, c, corner_usr...) == corner_usr
        end
    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)
    lastevent[] = "nothing"
    @test lastevent[] == "nothing"
    signal_emit(widget(c), "button-press-event", Bool, eventbutton(c, BUTTON_PRESS, 1))
    sleep(0.1)
    # FIXME: would prefer that this works on all Julia versions
    VERSION >= v"1.2.0" && @test lastevent[] == "press"
    signal_emit(widget(c), "button-release-event", Bool, eventbutton(c, GtkObservables.BUTTON_RELEASE, 1))
    sleep(0.1)
    sleep(0.1)
    VERSION >= v"1.2.0" && @test lastevent[] == "release"
    signal_emit(widget(c), "scroll-event", Bool, eventscroll(c, UP))
    sleep(0.1)
    sleep(0.1)
    VERSION >= v"1.2.0" && @test lastevent[] == "scroll"
    signal_emit(widget(c), "motion-notify-event", Bool, eventmotion(c, 0, UserUnit(20), UserUnit(15)))
    sleep(0.1)
    sleep(0.1)
    VERSION >= v"1.2.0" && @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)
    @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 = Observable(20), Observable(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)
    xsig[] = 100
    sleep(1)
    # Check that the displayed image is as expected
    if get(ENV, "CI", nothing) != "true" || !Sys.islinux() || VERSION < v"1.3" # broken on Travis
        fn = joinpath(tempdir(), "circled.png")
        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]
    end
    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 = GtkObservables.zoom(zr, 0.5)
    @test zrz.currentview.x == 26..75
    @test zrz.currentview.y == 21..60
    zrp = GtkObservables.pan_x(zrz, 0.2)
    @test zrp.currentview.x == 36..85
    @test zrp.currentview.y == 21..60
    zrp = GtkObservables.pan_x(zrz, -0.2)
    @test zrp.currentview.x == 16..65
    @test zrp.currentview.y == 21..60
    zrp = GtkObservables.pan_y(zrz, -0.2)
    @test zrp.currentview.x == 26..75
    @test zrp.currentview.y == 13..52
    zrp = GtkObservables.pan_y(zrz, 0.2)
    @test zrp.currentview.x == 26..75
    @test zrp.currentview.y == 29..68
    zrp = GtkObservables.pan_x(zrz, 1.0)
    @test zrp.currentview.x == 51..100
    @test zrp.currentview.y == 21..60
    zrp = GtkObservables.pan_y(zrz, -1.0)
    @test zrp.currentview.x == 26..75
    @test zrp.currentview.y == 1..40
    zrz2 = GtkObservables.zoom(zrz, 2.0001)
    @test zrz2 == zr
    zrz2 = GtkObservables.zoom(zrz, 3)
    @test zrz2 == zr
    zrz2 = GtkObservables.zoom(zrz, 1.9)
    @test zrz2.currentview.x == 4..97
    @test zrz2.currentview.y == 3..78
    zrz = GtkObservables.zoom(zr, 0.5, GtkObservables.XY{DeviceUnit}(50.5, 40.5))
    @test zrz.currentview.x == 26..75
    @test zrz.currentview.y == 21..60
    zrz = GtkObservables.zoom(zr, 0.5, GtkObservables.XY{DeviceUnit}(60.5, 30.5))
    @test zrz.currentview.x == 31..80
    @test zrz.currentview.y == 16..55
    zrr = GtkObservables.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 = Observable(zr)
    zrsig[] = (3:5, 4:7)
    zr = 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
    zrsig[] = XY(1..2, 3..4)
    zr = 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

@testset "More zoom/pan" begin
    ### Simulate the mouse clicks, etc. to trigger zoom/pan
    win = Window() |> (c = canvas(UserUnit))
    zr = Observable(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, 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))
    signal_emit(widget(c), "motion-notify-event", Bool,
                eventmotion(c, mask(1), UserUnit(10), UserUnit(4)))
    signal_emit(widget(c), "button-release-event", Bool,
                eventbutton(c, GtkObservables.BUTTON_RELEASE, 1, UserUnit(10), UserUnit(4)))
    @test zr[].currentview.x == 5..10
    @test zr[].currentview.y == 3..4
    # Ensure that the rubber band damage has been repaired
    if get(ENV, "CI", nothing) != "true" || !Sys.islinux() || VERSION < v"1.3" # broken on Travis
        fn = tempname()
        Cairo.write_to_png(getgc(c).surface, fn)
        imgout = load(fn)
        rm(fn)
        @test all(x->x==colorant"blue", imgout)
    end

    # Pan-drag
    signal_emit(widget(c), "button-press-event", Bool,
                eventbutton(c, BUTTON_PRESS, 1, UserUnit(6), UserUnit(3), 0))
    signal_emit(widget(c), "motion-notify-event", Bool,
                eventmotion(c, mask(1), UserUnit(7), UserUnit(2)))
    @test zr[].currentview.x == 4..9
    @test 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))
    @test zr[].currentview.x == 1..20
    @test zr[].currentview.y == 1..11

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

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

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

    destroy(win)
end

@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 = GtkObservables.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

examplepath = joinpath(dirname(dirname(@__FILE__)), "examples")
include(joinpath(examplepath, "imageviewer.jl"))
include(joinpath(examplepath, "widgets.jl"))
include(joinpath(examplepath, "drawing.jl"))


text:
struct OverlapChunker
    ρ::Float64
    w_max::Int
end

function pack_stripe(A::SparseMatrixCSC{Tv, Ti}, method::OverlapChunker, args...; n_nets = nothing, kwargs...) where {Tv, Ti}
    @inbounds begin
        # matrix notation...
        # i = 1:m rows, j = 1:n columns
        m, n = size(A)

        ρ = method.ρ
        w_max = method.w_max

        A_pos = A.colptr
        A_idx = A.rowval

        hst = zeros(Int, m)

        spl = Vector{Int}(undef, n + 1) # Column split locations
        if n_nets isa Nothing
            n_nets = Ref(Vector{Int}(undef, n)) # n_nets[j] is the corresponding number of distinct nonzero entries in the part
        else
            @assert n_nets isa Ref{Vector{Int}}
            n_nets[] = Vector{Int}(undef, n) # n_nets[j] is the corresponding number of distinct nonzero entries in the part
        end

        d = A_pos[2] - A_pos[1] #The number of distinct values in the part
        c = A_pos[2] - A_pos[1] #The cardinality of the first column in the part
        j = 1
        K = 0
        spl[1] = 1
        for q in A_pos[1] : A_pos[2] - 1
            i = A_idx[q]
            hst[i] = 1
        end
        for j′ = 2:n
            c′ = A_pos[j′ + 1] - A_pos[j′] #The cardinality of the candidate column
            d′ = d #Becomes the number of distinct values in the candidate part
            cc′ = 0 #The cardinality of the intersection between column j and j′
            for q = A_pos[j′] : A_pos[j′ + 1] - 1
                i = A_idx[q]
                h = hst[i]
                if abs(h) == j
                    cc′ += 1
                    hst[i] = -j′
                elseif j < h
                    hst[i] = j′
                elseif h < -j
                    cc′ += 1
                    hst[i] = -j′
                else
                    d′ += 1
                    hst[i] = j′
                end
            end
            w = j′ - j #Current block size
            if w == w_max || cc′ < ρ * min(c, c′)
                K += 1
                spl[K + 1] = j′
                n_nets[][K] = d
                j = j′
                d = c′
            else
                d = d′
            end
        end
        K += 1
        n_nets[][K] = d
        spl[K + 1] = n + 1
        resize!(spl, K + 1)
        resize!(n_nets[], K)
        return SplitPartition{Ti}(K, spl)
    end
end

ext: jl


text:
using ScanByte
using Test

include("byteset.jl")
include("scanbyte.jl")


text:
# Authors: Clemens Brunner
# License: BSD (3-clause)

module XDF

export read_xdf

using Logging: @info, @debug

CHUNK_TYPE = Dict(1=>"FileHeader",
                  2=>"StreamHeader",
                  3=>"Samples",
                  4=>"ClockOffset",
                  5=>"Boundary",
                  6=>"StreamFooter")

DATA_TYPE = Dict("int8"=>Int8,
                 "int16"=>Int16,
                 "int32"=>Int32,
                 "int64"=>Int64,
                 "float32"=>Float32,
                 "double64"=>Float64,
                 "string"=>String)

"""
    read_xdf(filename::AbstractString, sync::Bool=true)

Read XDF file and optionally sync streams (default true).
"""
function read_xdf(filename::AbstractString, sync::Bool=true)
    streams = Dict{Int, Any}()
    counter = Dict(zip(keys(CHUNK_TYPE), zeros(Int, length(CHUNK_TYPE))))  # count chunks

    open(filename) do io
        String(read(io, 4)) == "XDF:" || error("invalid magic bytes sequence")

        # first pass, determine array size for each stream
        while !eof(io)
            len = read_varlen_int(io)
            tag = read(io, UInt16)
            counter[tag] += 1
            @debug "Chunk $(sum(values(counter))): $(CHUNK_TYPE[tag]) ($tag), $len bytes"
            len -= sizeof(UInt16)

            if tag in (2, 3, 4, 6)  # read stream ID
                id = Int(read(io, UInt32))  # convert because Julia displays UInt32 as hex
                len -= sizeof(UInt32)
                @debug "    StreamID: $id"
            end

            if tag == 1  # FileHeader
                xml = String(read(io, len))
                @debug "    $xml"
            elseif tag == 2  # StreamHeader
                xml = String(read(io, len))
                @debug "    $xml"
                streams[id] = Dict("nchannels"=>findtag(xml, "channel_count", Int),
                                   "srate"=>findtag(xml, "nominal_srate", Float32),
                                   "name"=>findtag(xml, "name"),
                                   "type"=>findtag(xml, "type"))
                streams[id]["dtype"] = DATA_TYPE[findtag(xml, "channel_format")]
                streams[id]["data"] = 0
                streams[id]["time"] = Array{Float64}(undef, 0)
                streams[id]["clock"] = Float64[]
                streams[id]["offset"] = Float64[]
            elseif tag == 3  # Samples
                mark(io)
                nchannels = streams[id]["nchannels"]
                dtype = streams[id]["dtype"]
                nsamples = read_varlen_int(io)
                @debug "    nchans: $nchannels, nsamples: $nsamples, dtype: $dtype"
                streams[id]["data"] += nsamples
                reset(io)
                skip(io, len)
            elseif tag == 4  # ClockOffset
                push!(streams[id]["clock"], read(io, Float64))
                push!(streams[id]["offset"], read(io, Float64))
            elseif tag == 6  # StreamFooter
                xml = String(read(io, len))
                @debug "    $xml"
            else  # unknown chunk type
                skip(io, len)
            end
        end

        # second pass, read actual data for each stream into pre-allocated arrays
        index = Dict()
        for (id, stream) in streams
            dtype = stream["dtype"]
            nsamples = stream["data"]
            nchannels = stream["nchannels"]
            stream["data"] = Array{dtype}(undef, nsamples, nchannels)
            stream["time"] = Array{Float64}(undef, nsamples)
            index[id] = 1  # current sample index
        end
        seek(io, 4)  # go back to start of file (but skip the magic bytes)
        while !eof(io)
            len = read_varlen_int(io)
            tag = read(io, UInt16)
            len -= sizeof(UInt16)

            if tag != 3
                skip(io, len)
            else
                id = Int(read(io, UInt32))
                len -= sizeof(UInt32)
                mark(io)
                nchannels = streams[id]["nchannels"]
                dtype = streams[id]["dtype"]
                nsamples = read_varlen_int(io)
                for i in 1:nsamples
                    if read(io, UInt8) == 8  # optional timestamp available
                        streams[id]["time"][index[id]] = read(io, Float64)
                    else
                        delta = 1 / streams[id]["srate"]
                        previous = index[id] == 1 ? 0 : streams[id]["time"][index[id] - 1]
                        streams[id]["time"][index[id]] = previous + delta
                    end
                    if streams[id]["dtype"] === String
                        streams[id]["data"][index[id], :] .= String(read(io, read_varlen_int(io)))
                    else
                        streams[id]["data"][index[id], :] = reinterpret(dtype, read(io, sizeof(dtype) * nchannels))
                    end
                    index[id] += 1
                end
            end
        end
    end

    total = sum(values(counter))  # total number of chunks
    width = length(string(total))
    width_chunk = maximum([length(v) for v in values(CHUNK_TYPE)])
    msg = "File $filename contains $total chunks:\n"
    for (key, value) in sort(collect(counter))
        msg *= "- $(rpad(CHUNK_TYPE[key], width_chunk)) $(lpad(value, width))\n"
    end
    @info msg

    if sync
        for (id, stream) in streams
            offsets = hcat(stream["clock"], stream["offset"])
            stream["time"] = sync_clock(stream["time"], offsets)
        end
    end

    return streams
end

"Read variable-length integer."
function read_varlen_int(io::IO)
    nbytes = read(io, Int8)
    if nbytes == 1
        read(io, UInt8)
    elseif nbytes == 4
        read(io, UInt32)
    elseif nbytes == 8
        read(io, UInt64)
    end
end

"Find XML tag and return its content (optionally converted to specified type)."
function findtag(xml::AbstractString, tag::AbstractString, type=String::DataType)
    m = match(Regex("<$tag>(.*)</$tag>"), xml)
    content = isnothing(m) ? nothing : m[1]
    return isnothing(content) || type == String ? content : parse(type, content)
end

"Synchronize clock values by their given offsets."
function sync_clock(time::Array{Float64,1}, offsets::Array{Float64,2})
    x = hcat(ones(size(offsets, 1), 1), offsets[:, 1])
    y = offsets[:, 2]
    coefs = x \ y
    return time .+ (coefs[1] .+ coefs[2] .* time)
end

end


text:
module QuantumTensors

include("sparse_tensors.jl")

end # module


text:
# Monoid / Semigroup
# ==================

"""
    neutral
    neutral(::Type)
    neutral(_default_return_value) = neutral

Neutral element for `⊕`, also called "identity element". `neutral` is a function which can give you
the neutral element for a concrete type, or alternatively you can use it as a singleton value which combines with everything.

By default `neutral(type)` will return the generic `neutral` singleton. You can override it for your specific type to have a more specific neutral value.

We decided for name `neutral` according to https://en.wikipedia.org/wiki/Identity_element. Alternatives seem inappropriate
  - "identity" is already taken
  - "identity_element" seems to long
  - "I" is too ambiguous
  - "unit" seems ambiguous with physical units

Following Laws should hold
--------------------------

Left Identity

      ⊕(neutral(T), t::T) == t

Right Identity

      ⊕(t::T, neutral(T)) == t
"""
function neutral end
neutral(T::Type) = neutral  # we use the singleton type itself as the generic neutral value
neutral(a) = neutral(typeof(a))

"""
    combine(::T, ::T)::T  # overload this
    ⊕(::T, ::T)::T  # alias \\oplus
    combine(a, b, c, d, ...)  # using combine(a, b) internally

Associcative combinator operator.

The symbol `⊕` (\\oplus) following http://hackage.haskell.org/package/base-unicode-symbols-0.2.3/docs/Data-Monoid-Unicode.html

# Following Laws should hold

Associativity

    ⊕(::T, ⊕(::T, ::T)) == ⊕(⊕(::T, ::T), ::T)
"""
function combine end
const ⊕ = combine

combine(a, b, c, more...) = foldl(⊕, more, init=(a⊕b)⊕c)

# supporting `neutral` as generic neutral value
combine(::typeof(neutral), b) = b
combine(a, ::typeof(neutral)) = a
combine(::typeof(neutral), ::typeof(neutral)) = neutral

for reduce ∈ [:reduce, :foldl, :foldr]
  reduce_monoid = Symbol(reduce, "_monoid")
  _reduce_monoid = Symbol("_", reduce_monoid)

  @eval begin
    """
        $($reduce_monoid)(itr; init=TypeClasses.neutral)

    Combines all elements of `itr` using the initial element `init` if given and `TypeClasses.combine`.
    """
    function $reduce_monoid(itr; init=neutral)
      # taken from Base._foldl_impl

      # Unroll the while loop once to hopefully infer the element type at compile time
      y = iterate_named(itr)
      isnothing(y) && return init
      v = combine(init, y.value)
      while true
          y = iterate_named(itr, y.state)
          isnothing(y) && break
          v = combine(v, y.value)
      end
      return v
    end
  end
end

# Alternative
# ===========

"""
    orelse(a, b)  # overload this
    ⊘(a, b)  # alias \\oslash
    orelse(a, b, c, d, ...)  # using orelse(a, b) internally

Implements an alternative logic, like having two options a and b, taking the first valid one.
We decided for "orelse" instead of "alternatives" to highlight the intrinsic asymmetry in choosing.

The operator ⊘ (\\oslash) is choosen to have an infix operator which is similar to \\oplus, however clearly distinguishable, asymmetric, and somehow capturing a choice semantics.
The slash actually is used to indicate choice (at least in some languages, like German), and luckily \\oslash exists (and is not called \\odiv).
"""
function orelse end
const ⊘ = orelse

orelse(a, b, c, more...) = foldl(⊘, more, init=(a⊘b)⊘c)


text:
#=
Represents cells in the neighborhood of a face.
Has:
- cell indices
- cell center coords
- specified variable values (optional)
=#
struct Neighborhood
    cells::Array{Array{Int,1},1}        # indices for cells on left and right: [[L1,L2,...], [R1,R2,...]]
    centers::Array{Array{Float64,2},1}  # cell center coords [[L1x,L2x; L1y,L2y; L1z,L2z], [similar]]
    values::Union{Array{Array{Float64,2},1}, Array{Array{Float64,1},1}, Nothing}   # variable values in those cells or nothing

    Neighborhood(ind, center) = new(ind, center, nothing);
    Neighborhood(ind, center, var::Variable) = new(ind, center, get_neighborhood_variable_values(ind, var));
    Neighborhood(ind, center, var::Array{Array{Float64,1},1}) = new(ind, center, var);
    Neighborhood(ind, center, var::Array{Array{Float64,2},1}) = new(ind, center, var);
end

function get_neighborhood_variable_values(cells::Array{Array{Int,1},1}, var::Array{Variable,1})
    nvars = 0;
    for vi=1:length(var)
        nvars += size(var[vi].values,1);
    end
    vals = [zeros(nvars, length(cells[1])), zeros(nvars, length(cells[2]))];

    st = 1;
    for vi=1:length(var)
        en = st + size(var[vi].values,1) - 1;
        vals[1][st:en,:] = var[vi].values[:,cells[1]];
        vals[2][st:en,:] = var[vi].values[:,cells[2]];
        st = en+1;
    end

    return vals;
end
function get_neighborhood_variable_values(cells::Array{Array{Int,1},1}, var::Variable)
    return [var.values[:,cells[1]], var.values[:,cells[2]]];
end

Base.length(n::Neighborhood) = 1;

max_stars_repo_name: paralab/Finch

max_issues_repo_name: paralab/Finch

max_forks_repo_name: paralab/Finch


text:
import Catalyst
import Catalyst:@reaction_network
import ModelingToolkit
using Transducers

function gene_expression_system(; mean_s=50, corr_time_s=1.0, corr_time_x=0.1, u0=SA[mean_s, mean_s], dtimes=0:0.1:2.0)
    sn = @reaction_network begin
        κ, ∅ --> S
        λ, S --> ∅
    end κ λ

    xn = @reaction_network begin
        ρ, S --> X + S
        μ, X --> ∅
    end ρ μ

    λ = 1 / corr_time_s
    κ = mean_s * λ
    μ = 1 / corr_time_x
    ρ = μ

    ps = [κ, λ]
    px = [ρ, μ]

    SimpleSystem(sn, xn, u0, ps, px, dtimes)
end

function chemotaxis_system(;
        mean_L=20,
        num_receptors=10000,
        Y_tot=5000,
        L_timescale=1.0,
        LR_timescale=0.01,
        LR_ratio=0.5,
        Y_timescale=0.1,
        Y_ratio=1/6,
        q=0,
        dtimes=0:0.1:2.0
    )
    mean_LR = num_receptors * LR_ratio
    mean_R = num_receptors - mean_LR

    mean_Yp = Y_tot * Y_ratio
    mean_Y = Y_tot - mean_Yp

    eq_L = mean_L * exp(-q)

    sn = @reaction_network begin
        κ, ∅ --> L
        λ, L --> ∅
    end κ λ

    rn = @reaction_network begin
        ρ, L + R --> L + LR
        μ, LR --> R
    end ρ μ

    xn = @reaction_network begin
        δ, LR + Y --> Yp + LR
        χ, Yp --> Y
    end δ χ

    u0 = round.([mean_L, mean_R, mean_LR, mean_Y, mean_Yp])
    ps = [mean_L, 1 / L_timescale]

    ρ = inv(eq_L * LR_timescale * (1 + mean_R / mean_LR))
    μ = eq_L * ρ * mean_R / mean_LR

    pr = [ρ, μ]

    χ = inv(Y_timescale * (1 + Y_ratio))
    δ = χ * Y_ratio / mean_LR

    px = [δ, χ]

    ComplexSystem(sn, rn, xn, u0, ps, pr, px, dtimes)
end

function chemotaxis_parameters(;
    E₀ = 3.0,
    Kₐ = 500,
    Kᵢ = 25,
    δg = log(Kₐ/Kᵢ),
    δf = -1.5,
    k⁺ = 0.05,
    k⁺ₐ = k⁺,
    k⁺ᵢ = k⁺,
    k⁻ₐ = Kₐ * k⁺ₐ,
    k⁻ᵢ = Kᵢ * k⁺ᵢ,
    a_star = 0.5,
    γ = 1 / 10,
    k_B = (1 - a_star) * γ / abs(δf),
    k_R =      a_star  * γ / abs(δf)
)
    [
        E₀,
        δg,
        δf,
        k⁺ₐ,
        k⁺ᵢ,
        k⁻ₐ,
        k⁻ᵢ,
        k_B,
        k_R
    ]
end

function cooperative_chemotaxis_system(;
    lmax = 3,
    mmax = 9,
    n_clusters = 100,
    n_chey = 10_000,

    mean_l = 50,
    tau_l = 1.0,

    phosphorylate = 2000.0 / (n_chey * n_clusters),
    dephosphorylate = 2000.0 / (n_chey),
    dtimes = 0:0.1:20.0,
    varargs...
)
    sn = @reaction_network begin
        κ, ∅ --> L
        λ, L --> ∅
    end κ λ

    rn = Catalyst.make_empty_network()

    @Catalyst.parameters t E0 δg δf lba lbi lda ldi mda mbi ρ
    @Catalyst.variables L(t) Y(t) Yp(t)

    Catalyst.addspecies!(rn, L)

    Catalyst.addparam!(rn, E0)
    Catalyst.addparam!(rn, δg)
    Catalyst.addparam!(rn, δf)
    Catalyst.addparam!(rn, lba)
    Catalyst.addparam!(rn, lbi)
    Catalyst.addparam!(rn, lda)
    Catalyst.addparam!(rn, ldi)
    Catalyst.addparam!(rn, mda)
    Catalyst.addparam!(rn, mbi)

    xn = @reaction_network begin
        μ, Yp --> Y
    end μ

    Catalyst.addparam!(xn, ρ)

    spmap = Dict()
    for l=0:lmax, m=0:mmax
        receptor_species = ModelingToolkit.Num(ModelingToolkit.Variable{ModelingToolkit.FnType{Tuple{Any},Real}}(Symbol("R_", l, "_", m)))(t)

        spmap[(l, m)] = receptor_species

        Catalyst.addspecies!(xn, receptor_species)
        Catalyst.addspecies!(rn, receptor_species)
    end

    p_active(l, m) = 1 / (1 + exp(E0 + l*δg + m*δf))

    for l=0:lmax, m=0:mmax

        if l > 0
            ligand_bind = Catalyst.Reaction((lba * p_active(l-1, m) + lbi * (1 - p_active(l-1, m))) * (lmax + 1 - l), [spmap[(l-1, m)], L], [spmap[(l, m)], L])
            ligand_unbind = Catalyst.Reaction((lda * p_active(l, m) + ldi * (1 - p_active(l, m))) * l, [spmap[(l, m)]], [spmap[(l-1, m)]])
            Catalyst.addreaction!(rn, ligand_bind)
            Catalyst.addreaction!(rn, ligand_unbind)
        end

        if m > 0
            demethylate_active = Catalyst.Reaction(mda * p_active(l, m), [spmap[(l, m)]], [spmap[(l, m-1)]])
            methylate_inactive = Catalyst.Reaction(mbi * (1 - p_active(l, m-1)), [spmap[(l, m-1)]], [spmap[(l, m)]])

            Catalyst.addreaction!(rn, demethylate_active)
            Catalyst.addreaction!(rn, methylate_inactive)
        end

        # every receptor phosphorylates Y with rate ρ if the receptor is active
        # if the receptor is inactive, no phosphorylation can happen
        receptor = spmap[(l, m)]
        phosphorylation = Catalyst.Reaction(ρ * p_active(l, m), [Y, receptor], [Yp, receptor])
        Catalyst.addreaction!(xn, phosphorylation)
    end

    joint = merge(merge(sn, rn), xn)

    u0 = zeros(Int, length(spmap) + 3)
    u0[1] = round(Int, mean_l)
    u0[Catalyst.speciesmap(joint)[spmap[(0, 0)]]] = n_clusters
    u0[Catalyst.speciesmap(joint)[Y]] = n_chey

    ps = [mean_l/tau_l, 1.0/tau_l]
    pr = chemotaxis_parameters(; varargs...)
    px = [dephosphorylate, phosphorylate]

    ComplexSystem(sn, rn, xn, u0, ps, pr, px, dtimes; aggregator=DiffEqJump.RSSACR())
end

getname(sym) = String(ModelingToolkit.operation(sym).name)

function parse_receptor(species)
    sname = getname(species)
    mtch = match(r"([A-Z])_(\d+)_(\d+)", sname)
    if mtch !== nothing
        l = parse(Int, mtch.captures[2])
        m = parse(Int, mtch.captures[3])
        (l, m)
    else
        nothing
    end
end

function receptor_states(rs::ReactionSystem)
    smap = Catalyst.speciesmap(rs)
    xf = KeepSomething() do (species, index)
        result = parse_receptor(species)
        if result === nothing
            nothing
        else
            result => index
        end
    end
    smap |> xf
end

receptor_states(system::ComplexSystem) = receptor_states(reaction_network(system))

function active_receptors(conf::SRXconfiguration, system::ComplexSystem)
    rstates = receptor_states(system)

    E0 = system.pr[1]
    δg = system.pr[2]
    δf = system.pr[3]

    p_active(l, m) = 1 / (1 + exp(E0 + l*δg + m*δf))

    p_a = rstates |> Map(((l, m), i)::Pair -> i => p_active(l, m)) |> collect

    f = function(u)
        sum(p_a) do (i, p)
            p * u[i]
        end
    end

    merge_trajectories(conf.s_traj, conf.r_traj, conf.x_traj) |> Map((u,t,i)::Tuple -> (ensurevec(f(u)),t,i)) |> collect_trajectory
end

precompile(gene_expression_system, ())
precompile(chemotaxis_system, ())
precompile(cooperative_chemotaxis_system, ())


text:
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

using UnitCommitment, Test, LinearAlgebra
import UnitCommitment: _Violation, _offer, _query

@testset "find_violations" begin
    instance = UnitCommitment.read_benchmark("test/case14")
    for line in instance.lines, t in 1:instance.time
        line.normal_flow_limit[t] = 1.0
        line.emergency_flow_limit[t] = 1.0
    end
    isf = UnitCommitment._injection_shift_factors(
        lines = instance.lines,
        buses = instance.buses,
    )
    lodf = UnitCommitment._line_outage_factors(
        lines = instance.lines,
        buses = instance.buses,
        isf = isf,
    )
    inj = [1000.0 for b in 1:13, t in 1:instance.time]
    overflow = [0.0 for l in instance.lines, t in 1:instance.time]
    violations = UnitCommitment._find_violations(
        instance = instance,
        net_injections = inj,
        overflow = overflow,
        isf = isf,
        lodf = lodf,
        max_per_line = 1,
        max_per_period = 5,
    )
    @test length(violations) == 20
end


text:
using HTSEngineAPI
using Base.Test

include("htsengine.jl")


text:
###############################################################################
#
#   Zumbroich basis implementation

function J(k::Integer, p::Integer)
    if p == 2
        return ifelse(iszero(k), 0:0, 0:1)
    else
        pm1 = convert(Int, p) - 1
        return ifelse(iszero(k), 1:pm1, -pm1>>1:pm1>>1)
    end
end

function zumbroich_plain(n::Integer, m::Integer = 1)

    # following the description of
    # https://www.gap-system.org/Manuals/doc/ref/chap60_mj.html#X7F52BEA0862E06F2

    @assert iszero(last(divrem(n, m)))
    isone(n) && return [0]

    factors_n = Primes.factor(n)
    factors_m = Primes.factor(m)

    exps = Vector{Int}[]
    sizehint!(exps, 2 * length(factors_n))

    for (p, ν) in factors_n
        for k = factors_m[p]:ν-1
            push!(exps, [div(n, p^(k + 1)) * j for j in J(k, p)])
        end
    end

    w = (sum.(Iterators.product(exps...))) .% n
    return sort!(vec(w))
end

function zumbroich_direct(n::Integer)
    isone(n) && return [0]
    basis = [0]

    factor_n = Primes.factor(n)

    if iseven(n)
        for k = 1:factor_n[2]-1
            basis = union!(basis, basis .+ (n >> (k + 1)))
        end
    end

    for (p, ν) in factor_n
        p == 2 && continue
        for k = 0:ν-1
            ran = div(n, p^(k + 1)) .* J(k, p)
            new_basis = typeof(basis)()
            sizehint!(new_basis, length(basis) * length(ran))
            for b in basis
                append!(new_basis, ran .+ b)
            end
            basis = new_basis
        end
    end

    for (i, x) in enumerate(basis)
        basis[i] = mod(basis[i], n)
    end
    unique!(basis)
    sort!(basis)
    return basis
end

###############################################################################
#
# Forbidden residues, i.e. the complement of the Zumbroich basis

function _forbidden_residues(n, p, ν, q = p^ν)
    k = div(n, q)
    if isodd(p)
        a = (p^(ν - 1) - 1) >> 1
        return BitSet((k*(-a+q):k:k*(a+q)) .% q)
    else # iseven(p)
        return BitSet((k*(q>>1):k:k*(q-1)) .% q)
    end
end

struct ForbiddenResidues{I}
    primes_powers::Vector{Tuple{I,I,BitSet}}
end

function ForbiddenResidues(
    n::Integer,
    pf::Primes.Factorization{I} = Primes.factor(n),
) where {I}
    l = length(pf)
    primes_powers = Vector{Tuple{I,I,BitSet}}(undef, l)
    for (i, (p, ν)) in enumerate(pf)
        pν = p^ν
        primes_powers[i] = (p, pν, _forbidden_residues(n, p, ν, pν))
    end
    return ForbiddenResidues(primes_powers)
end

Base.length(fr::ForbiddenResidues) = length(fr.primes_powers)

Base.iterate(fr::ForbiddenResidues, s = 1) =
    (s > length(fr) ? nothing : (fr.primes_powers[s], s + 1))

@inline function Base.in(n::Integer, fr::ForbiddenResidues)
    @assert n ≥ 0
    for (_, q, residues) in fr
        n % q ∈ residues && return true
    end
    return false
end

function zumbroich_viacomplement(n::Integer, factor_n::Primes.Factorization)

    @debug "memoizing Zumbroich basis for ℚ(ζ_$n)"

    # following the wonderful documenting comments at the top of
    # https://github.com/gap-system/gap/blob/master/src/cyclotom.c

    forbidden = ForbiddenResidues(n, factor_n)

    isone(n) && return BitSet([0]), forbidden

    # exps = Vector{typeof(n)}(undef, Primes.totient(factor_n))
    # exps = BitSet(i for i in 1:n if !(i ∈ forbidden))

    exps = BitSet()
    sizehint!(exps, Primes.totient(factor_n))
    count = 0
    for i = 0:n-1
        i ∈ forbidden && continue
        count += 1
        # exps[count] = i
        push!(exps, i)
    end
    @assert count == length(exps)
    return exps, forbidden
end

import Memoize
import LRUCache

@static if VERSION >= v"1.2.0"
Memoize.@memoize LRUCache.LRU{Tuple{Int}, Tuple{BitSet, ForbiddenResidues{Int}}}(maxsize=1000) zumbroich_viacomplement(n::Int) =
    zumbroich_viacomplement(n, Primes.factor(n))
end

zumbroich_viacomplement(n::Integer) = zumbroich_viacomplement(n, Primes.factor(n))

zumbroich_basis(n::Integer) = first(zumbroich_viacomplement(n))


text:
#= Tests for Basket values =#
@testset "Basket" begin

# Simple arithmetic
basket_a = Basket(100USD)
basket_b = Basket(20EUR)
basket_c = basket_a + basket_b
basket_d = 4 * basket_c
basket_f = basket_a - basket_b
basket_g = basket_f / 4

# Mixed arithmetic
basket_h = basket_f + 20JPY
basket_i = basket_h - 100USD
basket_j = -20EUR - basket_i
basket_k = Basket() + basket_j
basket_l = (Basket(20USD) * 2 - 20CAD - 40USD) / 2 + 10CAD

@testset "Arithmetic" begin
    @test Basket([100USD, 200USD]) == Basket(300USD)
    @test basket_c == Basket([100USD, 20EUR])
    @test basket_d == Basket([400USD, 80EUR])
    @test basket_f == Basket([100USD, -20EUR])
    @test basket_g == Basket([25USD, -5EUR])

    # false positive tests
    @test basket_c ≠ Basket()
    @test basket_c ≠ Basket([100USD, 20EUR, 20GBP])
    @test Basket([100USD, 20EUR, 20GBP]) ≠ basket_c

    @test basket_h == Basket([100USD, -20EUR, 20JPY])
    @test basket_i == Basket([-20EUR, 20JPY])
    @test basket_j == -20JPY
    @test basket_j == basket_k
    @test basket_l == Basket()
    @test basket_l == 0EUR
    @test isempty(basket_l)

    # false positive tests
    @test basket_i ≠ -20EUR
    @test basket_j ≠ -20USD
    @test basket_i ≠ basket_j
    @test basket_l ≠ 1JPY
    @test basket_l ≠ Basket(1JPY)
end

# Basket constructor & zero
basket_m = Basket([basket_i, basket_j, 100JPY])
basket_n = zero(Basket)
basket_o = zero(basket_m)
basket_p = zero(Basket)
basket_q = Basket([basket_o, Basket([10USD, 20USD])])

@testset "Constructors" begin
    @test basket_m == Basket([-20EUR, 100JPY])
    @test basket_n == basket_o == Basket() == basket_p == zero(JPY)
    @test basket_q == 30USD
    @test one(Basket) ≡ 1

    @testset "Basket — Illegal Construction" begin
        @test_throws Exception Basket([1, 2, 3])
        @test_throws Exception Basket([1USD, (1USD, 2USD, 3)])
    end
end

# Iteration & access
@testset "As Collection" begin
    @test isempty(collect(Basket([100USD, -100USD])))
    @test basket_g[:EUR] == -5EUR
    @test haskey(Basket([10USD, -10CAD]), :CAD)
    @test !haskey(Basket([10USD, -10USD]), :USD)
    @test length(collect(Basket([1USD, 2USD]))) == 1
    @test length(collect(Basket([1USD, 1CAD]))) == 2
end

# Dynamic
basket_dyn = Basket() + basket_g

@testset "Dynamic" begin
    basket_dyn[:CAD] = 10CAD
    @test basket_dyn == Basket([25USD, -5EUR, 10CAD])
    @test haskey(basket_dyn, :USD)
    @test !haskey(basket_dyn, :JPY)
    push!(basket_dyn, 15JPY)
    @test basket_dyn == Basket([25USD, -5EUR, 10CAD, 15JPY])
    push!(basket_dyn, -10EUR)
    @test basket_dyn == Basket([25USD, -15EUR, 10CAD, 15JPY])
    push!(basket_dyn, -10CAD)
    @test !haskey(basket_dyn, :CAD)  # zero keys should act invisible
    @test length(collect(basket_dyn)) == 3
    push!(basket_dyn, Basket([25USD, 25EUR]))
    @test basket_dyn == Basket([50USD, 10EUR, 15JPY])
end

# Errors
@testset "Type Safety" begin
    @test_throws MethodError basket_a * basket_b
    @test_throws MethodError basket_a % basket_b
    @test_throws MethodError basket_a ÷ basket_b
    @test_throws MethodError basket_a ÷ 20USD
    @test_throws MethodError basket_a / basket_b
    @test_throws MethodError basket_a > basket_b
    @test_throws ArgumentError basket_dyn[:USD] = 100CAD
end

@testset "Conversion" begin
    b = Basket(100CAD)
    @test convert(Monetary, b) == 100CAD
    @test convert(Monetary{:CAD}, b) == 100CAD
    @test_throws Exception convert(Monetary{:USD}, b)
    @test_throws Exception convert(Monetary{:USD}, b + 100USD)
    @test_throws Exception convert(Monetary, b + 100USD)

    z = Basket()
    @test_throws Exception convert(Monetary, z)
    @test convert(Monetary{:CAD}, z) == 0CAD
end

end  # testset Basket


text:
"""
	branch_likelihood(t1, t2, L1, L2)

Given two branches of length `t1` and `t2`, with segments of length `L1` and `L2`, return the log-ratio of:
- the probability that the two branch lengths `t1` and `t2` are the same and equal to `(t1*L1 + t2*L2) / (L1 + L2)`.
- the probability that the two branch lenghts are different and equal to respectively `t1` and `t2`.
"""
function branch_likelihood(t1::Real, t2::Real, L1, L2)
	n1s = (t1*L1 + t2*L2) / (L1 + L2) * L1
	n2s = (t1*L1 + t2*L2) / (L1 + L2) * L2
	n1 = t1*L1
	n2 = t2*L2
	L = 0.
	if n1 != 0
		L += n1 - n1s + n1 * log(n1s/n1)
	else
		L += -n1s
	end
	if n2 != 0
		L += n2 - n2s + n2 * log(n2s/n2)
	else
		L += -n2s
	end
	return L
end

branch_likelihood(t1::Missing, t2::Real, L1, L2) = missing
branch_likelihood(t1::Real, t2::Missing, L1, L2) = missing
branch_likelihood(t1::Missing, t2::Missing, L1, L2) = 0.
branch_likelihood(t1, t2, L1, L2) = 0.

function conf_likelihood(conf::Array{Bool,1}, g::Graph, μ, trees; v=false, mode=:mutations)
	conf_likelihood_(TreeTools.divtime, conf, g, μ, trees, v=v)
end

"""
	conf_likelihood_(divfunction, conf::Array{Bool,1}, g::Graph, μ, trees; v=false)

For each leaf `n` remaining in `conf`, find the first non-trivial ancestors in graph `g`. For each pair of ancestors `a1` and `a2`, compare the likelihood  that the real branch length from `n` to `a1` and `a2` is the same, with the likelihood that it is different.
"""
function conf_likelihood_(divfunction, conf::Array{Bool,1}, g::Graph, seq_lengths, trees; v=false)
	L = 0.
	Z = 0.
	if length(g.leaves) + length(g.internals) == 1
		return L
	end
	for (i,s) in enumerate(conf)
		if s
			# Similar to compute_energy
			for k1 in 1:g.K
				# Tree node of k1 corresponding to graph leaf i
				tn1 = trees[k1].lleaves[g.labels[i]]
				# Going up in SplitGraph and tree at the same time
				a1 = g.leaves[i].anc[k1]
				ta1 = tn1.anc
				while is_trivial_split(a1.conf, conf) && !a1.isroot
					a1 = a1.anc::SplitNode
					ta1 = ta1.anc
				end
				for k2 in (k1+1):g.K
					# Tree node of k2 corresponding to graph leaf i
					tn2 = trees[k2].lleaves[g.labels[i]]
					# Going up in SplitGraph and tree at the same time
					a2 = g.leaves[i].anc[k2]
					ta2 = tn2.anc
					while is_trivial_split(a2.conf, conf) && !a2.isroot
						a2 = a2.anc
						ta2 = ta2.anc
					end
					# If a1.conf and a2.conf are equal, we just inferred that this branch is common to trees k1 and k2
					if (!get_resolve() && are_equal(a1.conf, a2.conf, conf)) ||
						(get_resolve() && are_equal_with_resolution(g, a1.conf, a2.conf, conf, k1, k2))
						tau1 = divfunction(tn1, ta1)
						tau2 = divfunction(tn2, ta2)
						# v && println("No inconsistency for leaf $i ($(g.labels[i]))")
						# v && println("Inferring that times $(tau1) on tree $k1 and $(tau2) on tree $k2 are the same.")
						dL = branch_likelihood(tau1, tau2, seq_lengths[k1], seq_lengths[k2])
						L += dL
						Z += 1.
						# v && println("Corresponding likelihood is $dL.")
					end
				end
			end
		else # branch above `i` is not shared, since we removes `i` from the trees
			for k1 in 1:g.K
				# Tree node of k1 corresponding to graph leaf i
				tn1 = trees[k1].lleaves[g.labels[i]]
				ta1 = tn1.anc
				for k2 in 1:g.K
					# Tree node of k2 corresponding to graph leaf i
					tn2 = trees[k2].lleaves[g.labels[i]]
					ta2 = tn2.anc
					# Removing log-likelihood ratio, since it's L(shared) / L(non shared)
					tau1 = divfunction(tn1, ta1)
					tau2 = divfunction(tn2, ta2)
					dL = -branch_likelihood(tau1, tau2, seq_lengths[k1], seq_lengths[k2])
					L += dL
					Z += 1.
				end
			end
		end
	end
	# v && println("--> Final likelihood $(L/max(Z,1))")
	return L/max(Z,1)
end
"""
"""
function conf_likelihood_times(conf::Array{Bool,1}, g::Graph, μ, trees; v=false)
	L = 0.
	for (i,s) in enumerate(conf)
		if s
			# Similar to compute energy
			for k1 in 1:g.K
				# Tree node of k1 corresponding to graph leaf i
				tn1 = trees[k1].lleaves[g.labels[i]]
				# Going up in SplitGraph and tree at the same time
				a1 = g.leaves[i].anc[k1]
				ta1 = tn1.anc
				while is_trivial_split(a1.conf, conf) && !a1.isroot
					a1 = a1.anc::SplitNode
					ta1 = ta1.anc
				end
				for k2 in (k1+1):g.K
					# Tree node of k2 corresponding to graph leaf i
					tn2 = trees[k2].lleaves[g.labels[i]]
					# Going up in SplitGraph and tree at the same time
					a2 = g.leaves[i].anc[k2]
					ta2 = tn2.anc
					while is_trivial_split(a2.conf, conf) && !a2.isroot
						a2 = a2.anc
						ta2 = ta2.anc
					end
					# If a1.conf and a2.conf are equal, we just inferred that this branch is common to trees k1 and k2
					if are_equal(a1.conf, a2.conf, conf)
						τ1 = TreeTools.divtime(tn1, ta1)
						τ2 = TreeTools.divtime(tn2, ta2)
						dL = branch_likelihood(μ[k1]*τ1, μ[k2]*τ2, μ[k1], μ[k2])
						# in dL, μ[k1]*τ1 should ultimately be replaced by n1 which is the observed number of mutations (same goes for 2)
						# What is currently there only makes sense for artificial data, where t is known exactly.
						L += dL
						v && println("No inconsistency for leaf $i ($(g.labels[i]))")
						v && println("Inferring that times $(τ1) on tree $k1 and $(τ2) on tree $k2 are the same.")
						v && println("Corresponding likelihood is $dL.")
					end
				end
			end
		end
	end
	return L
end

max_stars_repo_name: PierreBarrat/OldTreeKnit

max_issues_repo_name: PierreBarrat/OldTreeKnit

max_forks_repo_name: PierreBarrat/OldTreeKnit


text:
module RandomLineChart
using Plots

export gen

function gen(c=10,l=2)
    x = 1:c; y = rand(c,l);
    return plot(x,y, title="Random Line Chart")

end

end


text:
#using LinearAlgebra
#using DelimitedFiles
#using Statistics
#using Plots
#using DelimitedFiles
#using Random
#using SparseArrays
#using Arpack

#include("Geometry_twisted.jl")
#include("tight-binding_twisted.jl")
#include("valley_operator_twisted.jl")
#using PyCall
#plt = pyimport("matplotlib.pyplot")

#using Dates

function get_k_points(b1,b2,n_k)
    k_points=[b2*i/n_k+(b1-b2/2)/2 for i=0:n_k]
    return k_points
end

function get_band_twisted(H::get_hamiltonian,P::get_hamiltonian,k_points::Array{Array{Float64,1},1},n_bands,dim=2) # dim=1 for ribbon geometry
    n_k=size(k_points,1)
    momenta=Array{Float64}(undef,Int(n_k*n_bands))
    energies=Array{Float64}(undef,Int(n_k*n_bands))
    valley_polarization=Array{Float64}(undef,Int(n_k*n_bands))
    index=1
    for i=1:n_k
        k=k_points[i]
        if dim==2
            Hk=get_Hk(H,k)
            Pk=get_Pk(P,k)#get_layer_operator_twisted_bilayer(R_unitcell)#
        else
            Hk=get_Hk_ribbon(H,k)
            Pk=get_Pk_ribbon(P,k)#get_layer_operator_twisted_bilayer(R_unitcell)#
        end

        eigvals,eigvecs=eigs(Hk,nev=n_bands,sigma=0.0001,which=:LM,maxiter=10000)
        #F=eigen(Matrix(H))
        #eigvals,eigvecs=F.values,F.vectors
        #eigvals=real(eigvals)
        #n_bands=size(eigvals,1)
        for i_band=1:n_bands
            momenta[index]=i
            energies[index]=real(eigvals[i_band])
            psi=eigvecs[:,i_band]
            valley_polarization[index]=real(adjoint(psi)*Pk*psi)
            index+=1
        end
    end
    return momenta, energies, valley_polarization
end

function get_band_twisted_2(H::get_hamiltonian,P::get_hamiltonian,k_points::Array{Array{Float64,1},1},n_bands,dim=2) # dim=1 for ribbon geometry
    n_k=size(k_points,1)
    momenta=Array{Float64}(undef,Int(n_k*n_bands))
    energies=Array{Float64}(undef,Int(n_k*n_bands))
    valley_polarization=Array{Float64}(undef,Int(n_k*n_bands))
    index=1
    for i=1:n_k
        k=k_points[i]
        if dim==2
            Hk=get_Hk(H,k)
            Pk=get_Pk(P,k)#get_layer_operator_twisted_bilayer(R_unitcell)#
        else
            Hk=get_Hk_ribbon(H,k)
            Pk=get_Pk_ribbon(P,k)#get_layer_operator_twisted_bilayer(R_unitcell)#
        end
        Hk=Hk+0.05*Pk
        eigvals,eigvecs=eigs(Hk,nev=n_bands,sigma=0.0001,which=:LM,maxiter=10000)
        #F=eigen(Matrix(H))
        #eigvals,eigvecs=F.values,F.vectors
        #eigvals=real(eigvals)
        #n_bands=size(eigvals,1)
        for i_band=1:n_bands
            momenta[index]=i
            energies[index]=real(eigvals[i_band])
            psi=eigvecs[:,i_band]
            valley_polarization[index]=real(adjoint(psi)*Pk*psi)
            index+=1
        end
    end
    return momenta, energies, valley_polarization
end

function get_eigen_twisted(H::get_hamiltonian,k_points::Array{Array{Float64,1},1},n_bands)
    n_k=size(k_points,1)
    energies=Array{Complex}(undef,n_k*n_bands)
    wavefuncs=Array{Any}(undef,n_k*n_bands)
    t0=now()
    ii=1
    for i=1:n_k
        k=k_points[i]

        Hk=get_Hk(H,k)

        eigvals,eigvecs=eigs(Hk,nev=n_bands,sigma=0.0001,which=:LM,maxiter=10000)

        for i_band=1:n_bands
            energies[ii]=eigvals[i_band]
            wavefuncs[ii]=eigvecs[:,i_band]
            ii+=1
        end
    end

    t1=now()
    open("log.txt","a") do io
        println(io,"summed over", n_k," k points, time spent: ",t1-t0)
    end

    return energies, wavefuncs
end

function get_DOS(energies,omega,eta)
    DOS=0
    for E in energies
        DOS=DOS+1/(E-omega+im*eta)
    end
    DOS=-1/pi*imag(DOS)
    return real(DOS)
end

function get_LDOS(energies,wavefuncs,omega,eta)
    n_sites=size(wavefuncs[1],1)
    n_energies=size(energies,1)

    LDOS=zeros(Complex{Float64},n_sites)
    for i=1:n_energies
        psi=wavefuncs[i]
        E=energies[i]
        for j=1:n_sites
            LDOS[j]=LDOS[j]+abs2(psi[j])/(omega-E+im*eta)
        end
    end
    LDOS=real(-1/pi*imag(LDOS))
    return LDOS
end

function plot_LDOS(R_unitcell,LDOS)
    plt.figure(figsize=(6,6),dpi=80)
    N=size(R_unitcell,1)
    x=Array{Float64}(undef,N)
    y=Array{Float64}(undef,N)
    for i=1:N
        x[i]=R_unitcell[i][1]
        y[i]=R_unitcell[i][2]
    end
    #println(size(x,1),size(y,1))
    colors=[(1,0,0),(0,1,0),(0,0,1)]
    cm=colormap.LinearSegmentedColormap.from_list("my_list", colors)
    sc=plt.scatter(x, y,c=LDOS,cmap=cm)
    plt.colorbar(sc)
    plt.axis([-40,40,-40,40])
    plt.xlabel("")
    plt.ylabel("")
    #plt.savefig()
    plt.show()
end


text:
import SynapseClient: Activity

function _savecoverageplot(layer,format,outPath,args...)
    if format==:png
        args = (Theme(background_color=colorant"white"), args...)
    end

    pl = plot(layer,args...)

    filename = string(outPath,'.',format)
    if format==:png
        draw(PNG(filename, 29.7cm, 21cm), pl)
    elseif format==:pdf
        draw(PDF(filename, 29.7cm, 21cm), pl)
    elseif format==:svg
        draw(SVG(filename, 29.7cm, 21cm), pl)
    else
        error("Unknown image format $format")
    end
    filename
end

# outputs list of files created
function _coverageplots(samplePaths, sampleNames, outputFolder, outName, outFormats, log)
    print(log, "Computing read coverage...")
    coverage = readcoverage(samplePaths)
    println(log, "Done.")

    # check if these samples have one or multiple segments
    nbrSegments = length(coverage)

    segmentNames = [""]
    if nbrSegments>1
        seqs = sequences(BamFile(samplePaths[1])) # array of (segmentName,length)
        segmentNames = String[x[1] for x in seqs]
    end

    outFiles = Vector{String}()
    outSegment = Vector{String}()

    print(log, "Plotting...")
    for (i,segmentName) in enumerate(segmentNames)
        # Create DataFrame with all the data for the plot
        cov = coverage[i]
        maximum(cov)==0 && (cov[1]=1) # Gadfly fix for logscale when no sample has data.
        nbrSamples = size(cov,1)
        nbrPos     = size(cov,2)

        pos = repeat( (1:nbrPos)', nbrSamples, 1)
        name = repeat(sampleNames, 1, nbrPos)
        df = DataFrame(position=pos[:], coverage=cov[:], name=name[:])

        xTickStep = Int( 10^ceil(log10(nbrPos/2))/10 )
        xTicks = xTickStep:xTickStep:nbrPos

        covLayer = layer(df,x=:position,y=:coverage,color=:name,Geom.line)
        covPlotParams = (Scale.y_log10, Coord.cartesian(xmin=1,xmax=nbrPos,ymin=0), Guide.xticks(ticks=collect(xTicks)))

        filename = outName
        if !isempty(segmentName)
            filename = "$(filename)_$segmentName"
        end

        for f in outFormats
            push!(outFiles, _savecoverageplot(covLayer,f,joinpath(outputFolder,filename),covPlotParams...))
            push!(outSegment, segmentName)
        end
    end
    println(log, "Done.")

    outFiles, outSegment
end

# Create coverage plots for a single run
function coverageplots(syn, alignmentFolder, runName, doUpload, scriptFilename;
                       bamCache="bam", outputFolder="plots", outFormats=[:png, :pdf],
                       log=stdout, activityName = "Read Coverage Plots")
# --- Setup --------------------------------------------------------------------
    # Where to find sample .bam files. Synapse Folder ID or local folder. Synapse ID should point to "MyProject/Analysis/Alignment/Bam".
    bamPath = childpath(syn, alignmentFolder, "Bam")

    # Where to find sample .bam files. Synapse Folder ID or local folder. Synapse ID should point to "MyProject/Analysis/Alignment/Bam".
    scriptPath = childpath(syn, childpath(syn, alignmentFolder, "Scripts"), runName)

    alignLogFile = childpath(syn,scriptPath,"AlignUtils.log")
# --- Cleanup ------------------------------------------------------------------
    isdir(bamCache) || mkdir(bamCache)
    isdir(outputFolder) || mkdir(outputFolder)

# --- Make plots ---------------------------------------------------------------
    # find samples
    samplePaths, sampleNames = find_aligned(syn, bamPath, runName)
    references = referencefromlog(syn, alignLogFile, sampleNames)

    @assert !any(isempty,references)

    println(log, "Downloading samples")
    localPaths = map(x->localpath(syn,x,downloadLocation=bamCache,ifcollision="overwrite.local"), samplePaths)
    #sampleNamesShort = map(x->replace(x,Regex("^$(runName)_"),""), sampleNames) # remove runName_ from sample names

    filesToUpload = Vector{String}()
    fileSubFolder = Vector{String}()
    activities = Vector{Activity}() # Synapse Provenance

    # split samples by reference sequence
    for reference in unique(references)
        mask = references.==reference
        nbrSamples = count(mask)
        println(log, "Processing $nbrSamples samples for reference \"$reference\".")
        filesCreated, segments = _coverageplots(localPaths[mask], sampleNames[mask], outputFolder, "readcoverage_$(runName)_$reference", outFormats, log)
        for (f,s) in zip(filesCreated,segments)
            push!(filesToUpload, f)
            push!(fileSubFolder, s)
            act = Activity(activityName)
            used(act, samplePaths[mask])
            used(act, alignLogFile)
            push!(activities, act)
        end
    end

# --- Upload -------------------------------------------------------------------
    if doUpload
        basicAnalysesFolder = createchildfolder(syn, alignmentFolder, "Basic Analyses")
        runUploadFolder = createchildfolder(syn, basicAnalysesFolder, runName)
        plotUploadFolder = createchildfolder(syn, runUploadFolder, "Read Coverage")

        # upload script
        println(log, "Uploading script file.")
        scriptID = uploadiflocal(syn, scriptPath, scriptFilename, activityName=activityName)

        # upload plots
        println(log, "Uploading plots.")
        for (f, s, act) in zip(filesToUpload, fileSubFolder, activities)
            uploadFolder = isempty(s) ? plotUploadFolder : createchildfolder(syn, plotUploadFolder, s)
            executed(act, scriptID)
            uploadiflocal(syn, uploadFolder, f, act)
        end
    end

end


text:
###############################################################################
#
# Collection of methods generating random projectors.
#
# Improvements:
#  * random_matrix_4
#  * include sparse projectors
#  * rank-1 projectors
#
###############################################################################

"""
Random projector of size k x m
with a normal distribution
"""
function random_matrix_1(k, m)
  P = zeros(k, m) #init

  #normal distribution
  P = rand(k, m)

  return 1 / sqrt(k) * P
end

"""
Random projector of size k x m
with -1 or 1 both with probability 1/2
"""
function random_matrix_2(k, m)
  P = zeros(k, m) #init

  #-1 or 1 both with probability 1/2
  P = rand([-1, 1], (k, m))

  return 1 / sqrt(k) * P
end

"""
Random projector of size k x m
with -1,0,1 respectively with probability 1/6,4/6,1/6
"""
function random_matrix_3(k, m)
  P = zeros(k, m) #init

  #-1,0,1 respectively with probability 1/6,4/6,1/6
  P = rand([-1, 0, 0, 1], (k, m))

  return 1 / sqrt(k) * P
end

"""
Random projector of size k x m
with orthogonal projection on a random k-dimensional linear subspace of R^m
"""
function random_matrix_4(k, m)
  P = zeros(k, m) #init

  #orthogonal projection on a random k-dimensional
  #linear subspace of R^m
  #P = rand(k, m)
  throw("NotImplemented")
  return 1 / sqrt(k) * P
end

"""
Check the average sparsity of random projector over N random matrices of size k*m
"""
function random_projector_sparse(projector::Function, N::Int64, k::Int64, m::Int64)
  tot = 0 #nb of non-zero elements

  for i = 1:N
    tot += SparseArrays.nnz(SparseArrays.sparse(projector(k, m))) / (k * m)
  end

  return tot / N
end

"""
Check the average rank of random projector over N random matrices of size k*m
"""
function random_projector_rank(projector::Function, N::Int64, k::Int64, m::Int64)
  tot = 0 #average rank of the projector

  for i = 1:N
    tot += rank(projector(k, m))
  end

  return tot / N
end


text:
using Documenter
using DocumenterMarkdown
using GNSSTools

"""
    fix_eq_format(file_name)

Fix the format of the equations in the Markdown
files to be compatible with the GitLab equation
format requirement.
"""
function fix_eq_format(file_name)
    lines = readlines(file_name)
    for i in 1:length(lines)
        if occursin("\$\$", lines[i])
            lines[i] = "```math\n"
            i += 1
            while ~occursin("\$\$", lines[i])
                i += 1
            end
            lines[i] = "\n```"
        end
    end
    file = open(file_name, "w")
    for line in lines
        if line == ""
            write(file, "\n")
        else
            write(file, string(line, "\n"))
        end
    end
    close(file)
end

"""
    fix_format_in_files(directory, folders)

Fix the format of equations in each Markdown in
folders within `directory`.
"""
function fix_format_in_files(directory, folders)
    if isa(folders, Tuple)
        folers = [folders]
    end
    for folder in folders
        src = string(directory, "/", folder[1])
        dest = string(directory, "/", folder[2])
        files = readdir(src)
        for file in files
            if endswith(file, ".md")
                fix_eq_format(string(dest, "/", file))
            end
        end
    end
end

build_name = "GNSSTools.jl.wiki"

# Generate API and other page Markdown files
makedocs(sitename="GNSSTools Documentation",
         # format=Documenter.HTML(prettyurls=true)
         format=Markdown(),
         build=build_name,
         clean=false)

# Correct equation format in Markdown files
@info "Correcting equations."
directory = string(pwd(), "/docs")
folders = [("src", build_name),
           ("src/API", string(build_name, "/API"))]
fix_format_in_files(directory, folders)


text:
using Primes

tm(T::Type{<:Integer}) = typemax(T)
tm(::Type{BigInt}) = big"1000000000000000000000000000000000067"

@testset "construction and promotion" begin
    ZZ13 = Int / 13
    ZZ14 = Int / 14
    @test basetype(ZZmod{:p,Int}) == ZZ{Int}
    @test depth(ZZmod{13,BigInt}) == 2
    @test lcunit(ZZmod{13,Int}(7)) == 1
    @test lcunit(ZZmod{14,Int}(7)) == 1
    @test ZZmod{13,Int}(ZZ(8)) == ZZ13(8)
    @test ZZ13(1) + ZZ13(Int8(12)) == 0
    @test_throws MethodError ZZmod{13,Int}(1) - ZZmod{14,Int}(1)
    ZZp = ZZ(Int) / 13
    @test ZZmod{13,Int}(1) + ZZp(12) == 0
    @test typeof(ZZmod{13,Int}(1) + ZZp(1)) == ZZp
    @test typeof(promote(ZZ13(1), ZZp(1))) == Tuple{ZZp, ZZp}
    @test typeof(ZZp(1) + ZZ13(1)) == ZZp
    @test typeof(promote(ZZp(1), ZZ13(1))) == Tuple{ZZp, ZZp}
    @test ZZ(Int) / 13 == ZZmod{13,Int}
    @test BigInt/13 <: (ZZmod{X,BigInt} where X)
    @test convert(ZZp, ZZ13(Int8(7))) == ZZp(-6)
    @test_throws DomainError convert(ZZp, ZZmod{23,Int}(Int8(7)))
end

@testset "ZZmod{$m,$T}" for T in (UInt16, Int64, BigInt), m in (65, tm(T))

    while T != BigInt && isprime(m)
        m -= 2
    end
    m = T(m)

    if T != BigInt || m <= typemax(UInt128)
        phi = totient(m)
        p = factor(m).pe[end].first # the greatest prime factor of m
    else
        phi = (m-1)*4
        m = m * T(5)
    end

    n1 = T(19)
    while gcd(n1, m) != 1
        n1 += T(1)
    end
    n2 = tm(T) - T(16)
    while gcd(n2, m) != 1
        n2 += T(1)
    end

    ZZp = T / m
    @test typeof(modulus(ZZp)) == T
    isbitstype(T) && @test typeof(modulus(ZZmod{m,T})) == T
    @test typeof(modulus(ZZmod{3,T})) == T
    z = zero(ZZp)
    @test z == zero(z)
    @test iszero(z)
    @test !isone(z)
    o = one(z)
    @test isunit(o)
    @test isone(o)
    @test !iszero(o)
    z1 = ZZp(n1)
    z2 = ZZp(n2)
    zp = ZZp(p)
    @test isunit(z1)
    @test z1 + z1 == ZZp(T(2n1))

    @test z1 - z1 == z
    @test z1 - z2 == ZZp(m + mod(n1, m) - mod(n2, m))
    @test -z2 == ZZp(m - mod(n2, m))

    @test z1 * z1 == ZZp(T(n1 * n1))
    @test z1 * n1 == ZZp(T(n1 * n1))
    @test n1 * z1 == ZZp(T(n1 * n1))

    @test z1^2 == ZZp(n1^2)

    z3 = 2z1
    z4 = z3 + o
    @test z3 / z1 == ZZp(T(2))
    @test z4 / z1 == z4 * inv(z1)
    @test z1 \ z4 == z4 / z1

    @test z1^phi == o
    @test z2^(phi-1) == inv(z2)

    @test_throws DomainError inv(zp)
    @test_throws DomainError ZZmod{0,T}(0)

    if T != BigInt
        @test zero(ZZmod{m,Int}) == z
        @test one(ZZmod{m,Int}) == o
        @test ZZmod(n1, m) == ZZp(n1)
        @test hash(ZZmod(n1, m)) == hash(ZZp(n1))
    end

    @test ZZp(m-1) + ZZp(2) == ZZp(1)
    @test ZZp(1) - ZZp(2) == ZZp(m-1)
    @test ZZp(1) + ZZp(3) == ZZp(4)
    @test ZZp(3) - ZZp(2) == ZZp(1)

    @test ZZp(4)^-1 == inv(ZZp(4))
    @test ZZp(3)^1 == ZZp(3)
    @test ZZp(3)^0 == o

    @test copy(ZZp(5)) == ZZp(5)
    @test deg(z3) == 0
    @test div(z3, z3) == one(z3)
    @test rem(z3, z3) == zero(z3)

    @test sprint(show, z1) == "$(Int(n1))°"
    @test sprint(show, -z1) in ("$(big(m-n1))°", "-$(Int(n1))°")
end
@testset "constructors and type assertion" begin
    ZZp1 = ZZmod{17,Int8}
    p1 = ZZp1(-1)
    @test ZZp1(p1) === p1
    ZZp2 = ZZmod{17,Int}
    p2 = ZZp2(-1)
    @test ZZp1(p2) !== p2
    @test ZZp1(p2).val == p2.val
    @test ZZp1(128) == ZZp1(9)
    @test modulus(Int8/17) == 17
    @test modulus(BigInt/31) == 31
    @test (Int8/17)(1) == ZZp1(1)

    @test value((UInt8/17)(18)) === Int8(1)
    @test value.((ZZ/255).(126:128)) == [126, 127, -127]
    @test value.((ZZ/254).(126:128)) == [126, -127, -126]
end

@testset "auxiliary functions" begin
    @test CommutativeRings.invmod2(big"12", big"31") == big"13" # gcdx[2] > 0
    @test CommutativeRings.invmod2(big"15", big"31") == big"29" # gcdx[2] < 0
    @test CommutativeRings._unsigned(Int) == UInt
    @test CommutativeRings._unsigned(BigInt) == BigInt
    @test CommutativeRings._unsigned(big"-1") == big"-1"
end


text:
@testset "MontyHallEnv" begin

    rng = StableRNG(123)
    env = MontyHallEnv(; rng = rng)

    RLBase.test_interfaces!(env)
    RLBase.test_runnable!(env)

    n_win_car = 0
    N = 50_000

    for _ in 1:N
        a = rand(rng, action_space(env))
        env(a)
        env(a)
        if reward(env) == RLEnvs.REWARD_OF_CAR
            n_win_car += 1
        end
        reset!(env)
    end

    @test isapprox(n_win_car / N, 1 / 3; atol = 0.01)

end


text:
using Test
using ExtendableGrids

# This is the way how we can extend the grid type: define a new
# type, and either add the corresponding value or allow for
# lazy creation using the instantiate method

# Introduce an abstract type denoting the component
abstract type NodeCells <: AbstractGridAdjacency end

# Write an instantiate method which constructs the component if it is not available
ExtendableGrids.instantiate(grid, ::Type{NodeCells})=atranspose(grid[CellNodes])

function test_create_circle(;nref=5)
    nrad=10*2^nref
    grid=generate(
        points=reduce(hcat,[ [cos(2*pi*i/nrad), sin(2*pi*i/nrad)] for i=1:nrad]),
        bfaces=reduce(hcat,vcat([ [i,i+1] for i=1:nrad-1],[[nrad,1]])),
        bfaceregions=ones(nrad),
        regionpoints=[0.0 0.0;],
        regionnumbers=[1],
        regionvolumes=[0.1*2.0^(-2*nref)]
    )
end

function test_create_square(;nref=3)
    points=[[-2.0 -2];
            [2 -2];
            [2 2];
            [-2 2];
            [-1 -1];
            [1 -1];
            [1 1];
            [-1 1]]

    bfaces=[[0 1];
            [1 2];
            [2 3];
            [3 0];

            [4 5];
            [5 6];
            [6 7];
            [7 4]].+1

    bfaceregions=[1,2,3,4,5,5,5,5]

    regionpoints=[[-1.5 -1.5];
                  [0.0 0.0]];
    regionnumbers=[1,2]
    regionvolumes=[1.0*2.0^(-nref),1.0*2.0^(-nref)]

    grid=generate(points=points,
                  bfaces=bfaces,
                  bfaceregions=bfaceregions,
                  regionpoints=regionpoints,
                  regionnumbers=regionnumbers,
                  regionvolumes=regionvolumes,
                  flags="pqaAD")
end

function test_lazy(grid)
    nodecells=grid[NodeCells]
end

function sum_adj(a,n)
    sum=0
    for isource=1:num_sources(a)
        for itarget=1:num_targets(a,isource)
            sum+=a[itarget,isource]
        end
    end
    sum
end

function test_performance(grid)
    cellnodes=grid[CellNodes]
    n=num_sources(cellnodes)
    @time begin
        sum=0
        for isource=1:num_sources(cellnodes)
            for itarget=1:num_targets(cellnodes,isource)
                sum+=cellnodes[itarget,isource]
            end
        end
    end
    @time sum_adj(cellnodes,n)
    vadj=VariableTargetAdjacency(cellnodes)
    @time sum_adj(vadj,n)
    matrix=Matrix(cellnodes)
    @time sum_adj(matrix,n)
    nothing
end

function test_performance2(grid)
    sum=0
    @time begin
        sum=0
        for isource=1:num_sources(grid[CellNodes])
            for itarget=1:num_targets(grid[CellNodes],isource)
                sum+=grid[CellNodes][itarget,isource]
            end
        end
    end
    @show sum
    cellnodes=grid[CellNodes]
    @time begin
        sum=0
        for isource=1:num_sources(cellnodes)
            for itarget=1:num_targets(cellnodes,isource)
                sum+=cellnodes[itarget,isource]
            end
        end
    end
    @show sum
    nothing
end


text:
function git_make_commit(; commit_message::String)
    result = try
        cmd = `git commit -m $(commit_message)`
        p = pipeline(cmd; stdout=stdout, stderr=stderr)
        success(p)
    catch
        false
    end
    return result
end

function git_get_current_branch()::String
    return convert(String, strip(read(`git rev-parse --abbrev-ref HEAD`, String)))::String
end

function git_decide_master_branch(master_branch::DefaultBranch,
                                  default_branch::String)::String
    return convert(String, strip(default_branch))
end

function git_decide_master_branch(master_branch::AbstractString,
                                  default_branch::String)::String
    return convert(String, strip(master_branch))
end


text:
export rankx

"""
`rankx(A)` computes the rank of the exact matrix `A`
"""
function rankx(A::AbstractMatrix{T})::Int where {T}
    r, c = size(A)

    AA = deepcopy(A)
    try
        AA //= 1  # convert to rational if need be
        try
            rrefx!(AA)
        catch
            AA = big.(AA)
            rrefx!(AA)
        end
    catch
        rrefx!(AA)
    end

    count = 0
    for i = 1:r
        if any(AA[i, :] .!= 0)
            count += 1
        end
    end
    return count
end


text:
using POMDPGym, POMDPs, Test, PyCall

init_mujoco_render() # Required for visualization
pyimport("safety_gym")
mdp = GymPOMDP(Symbol("Safexp-PointGoal1"), version=:v0)

s = rand(initialstate(mdp))
@test size(s) == (60,)

info = Dict()
s, o, r = gen(mdp, s, 1, info=info)

info

max_forks_repo_name: ancorso/POMDPGym


text:
module AbstractModule
    function foo end
end