text:
using Makie: PriorityObservable, MouseButtonEvent, KeyEvent
using Makie: Not, And, Or

# rudimentary equality for tests
Base.:(==)(l::Exclusively, r::Exclusively) = l.x == r.x
Base.:(==)(l::Not, r::Not) = l.x == r.x
Base.:(==)(l::And, r::And) = l.left == r.left && l.right == r.right
Base.:(==)(l::Or, r::Or) = l.left == r.left && l.right == r.right

@testset "PriorityObservable" begin
    po = PriorityObservable(0)

    first = Observable(0.0)
    second = Observable(0.0)
    third = Observable(0.0)

    on(po, priority=1) do x
        sleep(0)
        first[] = time()
    end
    on(po, priority=0) do x
        sleep(0)
        second[] = time()
        return Consume(isodd(x))
    end
    on(po, priority=-1) do x
        sleep(0)
        third[] = time()
        return Consume(false)
    end

    x = setindex!(po, 1)
    @test x == true
    @test first[] < second[]
    @test third[] == 0.0

    x = setindex!(po, 2)
    @test x == false
    @test first[] < second[] < third[]
end

@testset "Events" begin
    @testset "Mouse and Keyboard state" begin
        events = Makie.Events()
        @test isempty(events.mousebuttonstate)
        @test isempty(events.keyboardstate)

        events.mousebutton[] = MouseButtonEvent(Mouse.left, Mouse.press)
        events.keyboardbutton[] = KeyEvent(Keyboard.a, Keyboard.press)
        @test events.mousebuttonstate == Set([Mouse.left])
        @test events.keyboardstate == Set([Keyboard.a])

        events.mousebutton[] = MouseButtonEvent(Mouse.right, Mouse.press)
        events.keyboardbutton[] = KeyEvent(Keyboard.b, Keyboard.press)
        @test events.mousebuttonstate == Set([Mouse.left, Mouse.right])
        @test events.keyboardstate == Set([Keyboard.a, Keyboard.b])

        events.mousebutton[] = MouseButtonEvent(Mouse.left, Mouse.release)
        events.keyboardbutton[] = KeyEvent(Keyboard.a, Keyboard.release)
        @test events.mousebuttonstate == Set([Mouse.right])
        @test events.keyboardstate == Set([Keyboard.b])

        events.mousebutton[] = MouseButtonEvent(Mouse.right, Mouse.release)
        events.keyboardbutton[] = KeyEvent(Keyboard.b, Keyboard.release)
        @test isempty(events.mousebuttonstate)
        @test isempty(events.keyboardstate)
    end

    @testset "ispressed" begin
        events = Makie.Events()
        @test isempty(events.mousebuttonstate)
        @test isempty(events.keyboardstate)

        events.mousebutton[] = MouseButtonEvent(Mouse.left, Mouse.press)
        events.keyboardbutton[] = KeyEvent(Keyboard.a, Keyboard.press)

        # Buttons
        @test ispressed(events, Keyboard.a)
        @test !ispressed(events, Keyboard.b)
        @test ispressed(events, Mouse.left)
        @test !ispressed(events, Mouse.right)

        # Collections
        @test ispressed(events, (Keyboard.a, ))
        @test ispressed(events, [Keyboard.a, ])
        @test ispressed(events, Set((Keyboard.a, )))

        @test !ispressed(events, (Keyboard.a, Keyboard.b))
        @test !ispressed(events, [Keyboard.a, Keyboard.b])
        @test !ispressed(events, Set((Keyboard.a, Keyboard.b)))

        @test ispressed(events, (Keyboard.a, Mouse.left))
        @test ispressed(events, [Keyboard.a, Mouse.left])
        @test ispressed(events, Set((Keyboard.a, Mouse.left)))

        # Boolean
        @test ispressed(events, Keyboard.a & Mouse.left)
        @test !ispressed(events, Keyboard.a & Mouse.right)
        @test ispressed(events, Keyboard.a & !Mouse.right)
        @test !ispressed(events, !Keyboard.a & Mouse.left)

        @test ispressed(events, Keyboard.a | Mouse.left)
        @test ispressed(events, Keyboard.a | Mouse.right)
        @test !ispressed(events, Keyboard.b | Mouse.right)
        @test ispressed(events, Keyboard.b | !Mouse.right)

        # Exclusively
        @test !ispressed(events, Exclusively(Keyboard.a))
        @test !ispressed(events, Exclusively(Mouse.left))

        @test ispressed(events, Exclusively((Keyboard.a, Mouse.left)))
        @test ispressed(events, Exclusively([Keyboard.a, Mouse.left]))
        @test ispressed(events, Exclusively(Set((Keyboard.a, Mouse.left))))

        @test Exclusively(Keyboard.a & Mouse.left) == Exclusively((Keyboard.a, Mouse.left))
        @test ispressed(events, Exclusively(Keyboard.a & Mouse.left))
        @test Exclusively(Keyboard.a | Mouse.left) == Makie.Or(Exclusively(Keyboard.a), Exclusively(Mouse.left))
        @test !ispressed(events, Exclusively(Keyboard.a | Mouse.left))

        expr = Mouse.left & (Keyboard.a | (Keyboard.a & Keyboard.b))
        lowered = Or(
            Exclusively((Mouse.left, Keyboard.a)),
            Exclusively((Mouse.left, Keyboard.a, Keyboard.b))
        )
        @test ispressed(events, expr)
        @test Exclusively(expr) == lowered
        @test ispressed(events, Exclusively(expr))

        events.keyboardbutton[] = KeyEvent(Keyboard.b, Keyboard.press)
        @test ispressed(events, expr)
        @test ispressed(events, Exclusively(expr))

        events.keyboardbutton[] = KeyEvent(Keyboard.c, Keyboard.press)
        @test ispressed(events, expr)
        @test !ispressed(events, Exclusively(expr))

        events.keyboardbutton[] = KeyEvent(Keyboard.a, Keyboard.release)
        @test !ispressed(events, expr)
        @test !ispressed(events, Exclusively(expr))

        # Bools
        @test ispressed(events, true)
        @test !ispressed(events, false)

        @test Exclusively(true) == true
        @test Exclusively(false) == false

        for x in (Keyboard.a, Mouse.left)
            @test true & x == x
            @test x & true == x
            @test false & x == false
            @test x & false == false

            @test true | x == true
            @test x | true == true
            @test false | x == x
            @test x | false == x
        end
    end

    # This testset is based on the results the current camera system has. If
    # cam3d! is updated this is likely to break.
    @testset "cam3d!" begin
        scene = Scene(resolution=(800, 600));
        e = events(scene)
        cam3d!(scene, fixed_axis=true, cad=false, zoom_shift_lookat=false)
        cc = cameracontrols(scene)

        # Verify initial camera state
        @test cc.lookat[]       == Vec3f(0)
        @test cc.eyeposition[]  == Vec3f(3)
        @test cc.upvector[]     == Vec3f(0, 0, 1)

        # Rotation
        # 1) In scene, in drag
        e.mouseposition[] = (400, 250)
        e.mousebutton[] = MouseButtonEvent(Mouse.left, Mouse.press)
        e.mouseposition[] = (600, 250)
        @test cc.lookat[]       ≈ Vec3f(0)
        @test cc.eyeposition[]  ≈ Vec3f(4.14532, -0.9035063, 3.0)
        @test cc.upvector[]     ≈ Vec3f(0, 0, 1)

        # 2) Outside scene, in drag
        e.mouseposition[] = (1000, 450)
        @test cc.lookat[]       ≈ Vec3f(0)
        @test cc.eyeposition[]  ≈ Vec3f(-2.8912058, -3.8524969, -1.9491522)
        @test cc.upvector[]     ≈ Vec3f(-0.5050875, -0.6730229, 0.5403024)

        # 3) not in drag
        e.mousebutton[] = MouseButtonEvent(Mouse.left, Mouse.release)
        e.mouseposition[] = (400, 250)
        @test cc.lookat[]       ≈ Vec3f(0)
        @test cc.eyeposition[]  ≈ Vec3f(-2.8912058, -3.8524969, -1.9491522)
        @test cc.upvector[]     ≈ Vec3f(-0.5050875, -0.6730229, 0.5403024)

        # Reset state so this is indepentent from the last checks
        scene = Scene(resolution=(800, 600));
        e = events(scene)
        cam3d!(scene, fixed_axis=true, cad=false, zoom_shift_lookat=false)
        cc = cameracontrols(scene)

        # Verify initial camera state
        @test cc.lookat[]       == Vec3f(0)
        @test cc.eyeposition[]  == Vec3f(3)
        @test cc.upvector[]     == Vec3f(0, 0, 1)

        # translation
        # 1) In scene, in drag
        e.mousebutton[] = MouseButtonEvent(Mouse.right, Mouse.press)
        e.mouseposition[] = (600, 250)
        @test cc.lookat[]       ≈ Vec3f(5.4697413, -3.3484206, -2.1213205)
        @test cc.eyeposition[]  ≈ Vec3f(8.469742, -0.34842062, 0.8786795)
        @test cc.upvector[]     ≈ Vec3f(0.0, 0.0, 1.0)

        # 2) Outside scene, in drag
        e.mouseposition[] = (1000, 450)
        @test cc.lookat[]       ≈ Vec3f(9.257657, -5.4392805, -3.818377)
        @test cc.eyeposition[]  ≈ Vec3f(12.257658, -2.4392805, -0.81837714)
        @test cc.upvector[]     ≈ Vec3f(0.0, 0.0, 1.0)

        # 3) not in drag
        e.mousebutton[] = MouseButtonEvent(Mouse.right, Mouse.release)
        e.mouseposition[] = (400, 250)
        @test cc.lookat[]       ≈ Vec3f(9.257657, -5.4392805, -3.818377)
        @test cc.eyeposition[]  ≈ Vec3f(12.257658, -2.4392805, -0.81837714)
        @test cc.upvector[]     ≈ Vec3f(0.0, 0.0, 1.0)

        # Reset state
        scene = Scene(resolution=(800, 600));
        e = events(scene)
        cam3d!(scene, fixed_axis=true, cad=false, zoom_shift_lookat=false)
        cc = cameracontrols(scene)

        # Verify initial camera state
        @test cc.lookat[]       == Vec3f(0)
        @test cc.eyeposition[]  == Vec3f(3)
        @test cc.upvector[]     == Vec3f(0, 0, 1)
        @test cc.zoom_mult[]    == 1f0

        # Zoom
        e.scroll[] = (0.0, 4.0)
        @test cc.lookat[]       ≈ Vec3f(0)
        @test cc.eyeposition[]  ≈ Vec3f(3)
        @test cc.upvector[]     ≈ Vec3f(0.0, 0.0, 1.0)
        @test cc.zoom_mult[]    ≈ 0.6830134f0

        # should not work outside the scene
        e.mouseposition[] = (1000, 450)
        e.scroll[] = (0.0, 4.0)
        @test cc.lookat[]       ≈ Vec3f(0)
        @test cc.eyeposition[]  ≈ Vec3f(3)
        @test cc.upvector[]     ≈ Vec3f(0.0, 0.0, 1.0)
        @test cc.zoom_mult[]    ≈ 0.6830134f0
    end

    @testset "mouse state machine" begin
        scene = Scene(resolution=(800, 600));
        e = events(scene)
        bbox = Observable(Rect2(200, 200, 400, 300))
        msm = addmouseevents!(scene, bbox, priority=typemax(Int8))
        eventlog = MouseEvent[]
        on(x -> begin push!(eventlog, x); false end, msm.obs)

        e.mouseposition[] = (0, 200)
        @test isempty(eventlog)

        # move inside
        e.mouseposition[] = (300, 200)
        @test length(eventlog) == 1
        @test eventlog[1].type == MouseEventTypes.enter
        @test eventlog[1].px == Point2f(300, 200)
        @test eventlog[1].prev_px == Point2f(0, 200)
        empty!(eventlog)

        # over
        e.mouseposition[] = (300, 300)
        @test length(eventlog) == 1
        @test eventlog[1].type == MouseEventTypes.over
        @test eventlog[1].px == Point2f(300, 300)
        @test eventlog[1].prev_px == Point2f(300, 200)
        empty!(eventlog)

        for button in (:left, :middle, :right)
            # click
            e.mousebutton[] = MouseButtonEvent(getfield(Mouse, button), Mouse.press)
            e.mousebutton[] = MouseButtonEvent(getfield(Mouse, button), Mouse.release)
            @test length(eventlog) == 3
            for (i, t) in enumerate((
                    getfield(MouseEventTypes, Symbol(button, :down)),
                    getfield(MouseEventTypes, Symbol(button, :click)),
                    getfield(MouseEventTypes, Symbol(button, :up))
                ))
                @test eventlog[i].type == t
                @test eventlog[i].px == Point2f(300, 300)
                @test eventlog[i].prev_px == Point2f(300, 300)
            end
            empty!(eventlog)

            # doubleclick
            e.mousebutton[] = MouseButtonEvent(getfield(Mouse, button), Mouse.press)
            e.mousebutton[] = MouseButtonEvent(getfield(Mouse, button), Mouse.release)
            @test length(eventlog) == 3
            for (i, t) in enumerate((
                    getfield(MouseEventTypes, Symbol(button, :down)),
                    getfield(MouseEventTypes, Symbol(button, :doubleclick)),
                    getfield(MouseEventTypes, Symbol(button, :up))
                ))
                @test eventlog[i].type == t
                @test eventlog[i].px == Point2f(300, 300)
                @test eventlog[i].prev_px == Point2f(300, 300)
            end
            empty!(eventlog)

            # triple click = click
            e.mousebutton[] = MouseButtonEvent(getfield(Mouse, button), Mouse.press)
            e.mousebutton[] = MouseButtonEvent(getfield(Mouse, button), Mouse.release)
            @test length(eventlog) == 3
            for (i, t) in enumerate((
                    getfield(MouseEventTypes, Symbol(button, :down)),
                    getfield(MouseEventTypes, Symbol(button, :click)),
                    getfield(MouseEventTypes, Symbol(button, :up))
                ))
                @test eventlog[i].type == t
                @test eventlog[i].px == Point2f(300, 300)
                @test eventlog[i].prev_px == Point2f(300, 300)
            end
            empty!(eventlog)

            # drag
            e.mousebutton[] = MouseButtonEvent(getfield(Mouse, button), Mouse.press)
            e.mouseposition[] = (500, 300)
            e.mouseposition[] = (700, 200)
            e.mousebutton[] = MouseButtonEvent(getfield(Mouse, button), Mouse.release)
            @test length(eventlog) == 6
            prev_px = Point2f[(300, 300), (300, 300), (300, 300), (500, 300), (700, 200), (700, 200)]
            px      = Point2f[(300, 300), (500, 300), (500, 300), (700, 200), (700, 200), (700, 200)]
            for (i, t) in enumerate((
                    getfield(MouseEventTypes, Symbol(button, :down)),
                    getfield(MouseEventTypes, Symbol(button, :dragstart)),
                    getfield(MouseEventTypes, Symbol(button, :drag)),
                    getfield(MouseEventTypes, Symbol(button, :drag)),
                    getfield(MouseEventTypes, Symbol(button, :dragstop)),
                    getfield(MouseEventTypes, :out),
                    # TODO this is kinda missing an "up outside"
                ))
                @test eventlog[i].type == t
                @test eventlog[i].px == px[i]
                @test eventlog[i].prev_px == prev_px[i]
            end
            e.mouseposition[] = (300, 300)
            empty!(eventlog)
        end

        # TODO: This should probably be:
        # left down > right down > right click > right up > left up
        e.mousebutton[] = MouseButtonEvent(Mouse.left, Mouse.press)
        e.mousebutton[] = MouseButtonEvent(Mouse.right, Mouse.press)
        e.mousebutton[] = MouseButtonEvent(Mouse.right, Mouse.release)
        e.mousebutton[] = MouseButtonEvent(Mouse.left, Mouse.release)
        @test length(eventlog) == 3
        @test eventlog[1].type == MouseEventTypes.leftdown
        @test eventlog[2].type == MouseEventTypes.leftclick
        @test eventlog[3].type == MouseEventTypes.leftup
        empty!(eventlog)

        # double left up? :(
        e.mousebutton[] = MouseButtonEvent(Mouse.left, Mouse.press)
        e.mousebutton[] = MouseButtonEvent(Mouse.right, Mouse.press)
        e.mousebutton[] = MouseButtonEvent(Mouse.left, Mouse.release)
        e.mousebutton[] = MouseButtonEvent(Mouse.right, Mouse.release)
        @test length(eventlog) == 4
        @test eventlog[1].type == MouseEventTypes.leftdown
        @test eventlog[2].type == MouseEventTypes.leftdoubleclick
        @test eventlog[3].type == MouseEventTypes.leftup
        @test eventlog[4].type == MouseEventTypes.leftup
        empty!(eventlog)

        # This should probably be a leftdragstop on right down
        e.mouseposition[] = (300, 300)
        empty!(eventlog)
        e.mousebutton[] = MouseButtonEvent(Mouse.left, Mouse.press)
        e.mouseposition[] = (350, 350)
        e.mousebutton[] = MouseButtonEvent(Mouse.right, Mouse.press)
        e.mouseposition[] = (350, 400)
        e.mousebutton[] = MouseButtonEvent(Mouse.right, Mouse.release)
        e.mouseposition[] = (400, 400)
        e.mousebutton[] = MouseButtonEvent(Mouse.left, Mouse.release)
        @test length(eventlog) == 7
        @test eventlog[1].type == MouseEventTypes.leftdown
        @test eventlog[2].type == MouseEventTypes.leftdragstart
        @test eventlog[3].type == MouseEventTypes.leftdrag
        @test eventlog[4].type == MouseEventTypes.leftdrag
        @test eventlog[5].type == MouseEventTypes.over
        @test eventlog[6].type == MouseEventTypes.leftdragstop
        @test eventlog[7].type == MouseEventTypes.leftup
        @test eventlog[1].px == Point2f(300, 300)
        @test eventlog[2].px == Point2f(350, 350)
        @test eventlog[3].px == Point2f(350, 350)
        @test eventlog[4].px == Point2f(350, 400)
        @test eventlog[5].px == Point2f(400, 400)
        @test eventlog[6].px == Point2f(400, 400)
        @test eventlog[7].px == Point2f(400, 400)
        empty!(eventlog)
    end
end


text:
# Pasquill-Gifford dispersion correlations
include("pasquill_gifford.jl")

# correlations for Monin-Obukhov length
include("monin_obukhov.jl")

# correlations for windspeed
include("wind_profile.jl")

ext: jl


text:
# Coefficients calculated with https://github.com/simonbyrne/Remez.jl

@inline function approx_sin8(x::Union{T,Vec{<:Any,T},VecUnroll{<:Any,<:Any,T}}) where {T <: Real}
  # poly(x) ≈ (xʳ = sqrt(x); sin((xʳ*π)/2)/xʳ)
    x² = x * x
    c0 = T(2.22144146907918312350794048535203995923494010677251491220479906920966593121882)
    c1 = T(-0.9135311874994298224944388934705417261765270518848695099428083902179199377101094)
    c2 = T(0.1127023928584587596980569269678174942915399051122642981118394498722218063783927)
    c3 = T(-0.006621000193853498898990183110992108352486751535892362800909323879419896057043918)
    c4 = T(0.0002268980994233557245363541171760472387529757765245978583128895641498725296271051)
    c5 = T(-5.089532691384021959110856232473979525292167742059549332987900223626864039349914e-06)
    c6 = T(8.049906344315649609313027324977744156866597923196983008950128144505665619892402e-08)
    c7 = T(-9.453796623737636858301034347145347814693537235132105505794304057287442064404052e-10)
    c8 = T(8.320735422342537824261297491878000532726851750329165722059039816086266315937799e-12)
    p = vfmadd(vfmadd(vfmadd(
        vfmadd(vfmadd(vfmadd(
            vfmadd(vfmadd(
                c8, x², c7),
            x², c6), x², c5), x², c4),
        x², c3), x², c2), x², c1), x², c0)
    p * x
end
# @inline function approx_sin12(x::Union{T,Vec{<:Any,T}}) where {T <: Real}
#     # poly(x) ≈ (xʳ = sqrt(x); sin((xʳ*π)/2)/xʳ)
#     x² = x * x
#     c0 = T(1.570796326794896619231321691639751442087433306473273974291471596002143089408967)
#     c1 = T(-0.6459640975062462536557565638714840878228221616991079162636960728776200926827234)
#     c2 = T(0.07969262624616704512050554673779356754386556916433562280307069303594093234088455)
#     c3 = T(-0.004681754135318688100685379129717344612020387136900677528796079362820120586116841)
#     c4 = T(0.0001604411847873598218714490487175497697128038098258443135990371062725913378293313)
#     c5 = T(-3.598843235212085330760986854968071769532101571696000087022804157171548649851493e-06)
#     c6 = T(5.692172921967922014173914534767998233640265149949183614622512228877789263783426e-08)
#     c7 = T(-6.688035109809916561166255796732305220204238240380007953680372148880958984800929e-10)
#     c8 = T(6.066935730769290440108765932783579488881358689623800108082605751157383050893581e-12)
#     c9 = T(-4.377065417731331420103035981625834800685920271728367262547966148291799276629881e-14)
#     c10 = T(2.571418016198708615875917881136145309875324227233510421118020948686662677979356e-16)
#     c11 = T(-1.253592449512705798908955136513569509617634496103293074276456663908804526347008e-18)
#     c12 = T(5.044383456268885650704416950405914330732446213362030631794638992715099861013542e-21)
#     p = vfmadd(vfmadd(vfmadd(
#         vfmadd(vfmadd(vfmadd(
#             vfmadd(vfmadd(vfmadd(
#                 vfmadd(vfmadd(vfmadd(
#                     c12, x², c11), x², c10),
#                 x², c9), x², c8), x², c7),
#             x², c6), x², c5), x², c4),
#         x², c3), x², c2), x², c1), x², c0)
#     p * x
# end

@inline suboneopenconst(::Type{Float32}) = 1.9999999f0
@inline suboneopenconst(::Type{Float64}) = 1.9999999999999998
@inline function randsincos(u, ::Type{T}) where {T}
  # return SLEEFPirates.sincos(mask(u, T))
  r = floatbitmask(u, T)
  ooc = oneopenconst(T)
  sininput = vsub(r, ooc)
  cosinput = vfnmadd(ooc, r, suboneopenconst(T))
  sc = data(approx_sin8(VecUnroll((sininput, cosinput))))
  s = copysign( getfield(sc,1), reinterpret(T, u))
  c = copysign( getfield(sc,2), reinterpret(T, u << 1))
  s, c
end

# @inline function log12_7(x) # each extra coef cuts max error by about 6.5
#     c0 = -3.245537891437475818527978529229908008038541532632077901681793316955253799627853
#     c1 = 7.133969761783261596308839380142764345357825207927592180328308957632143143761325
#     c2 = -7.494130176731051024066811749217581442825262641998555203442250312457063219780548
#     c3 = 5.781439120377507449162563050101602606694402138734039600998751617511391254917479
#     c4 = -2.985021102358628224549723815823770715482850258619037789336465429092098006063079
#     c5 = 0.9780948488102108131759607721708499684604503249092972679349980253346585174787013
#     c6 = -0.1839396423223307845519460189847030764334420832555490882187076088746804014578904
#     c7 = 0.01512535916800840093163249616452966347663616377265149854210868241945849663079766
#     vfmadd(vfmadd(vfmadd(vfmadd(vfmadd(vfmadd(vfmadd(c7, x, c6), x, c5), x, c4), x, c3), x, c2), x, c1), x, c0)
# end

# # @inline function log12_8(x)
# #     c0 = -3.425393083666334067790108512545303457895869147629969342312722110773152768787745
# #     c1 = 8.154804072298105118614242101352098813998371180480781803351886784650194633423925
# #     c2 = -10.00713525074289150060544073624417642119577303442852729267703418853804728270245
# #     c3 = 9.285983673564252215780879208501085478907382239737051965242871871829356007506347
# #     c4 = -6.013432610610964086674378270192946254591450512746857671796983783112863051993776
# #     c5 = 2.638773868907455181856795850097739510981561473919572777387590013143564751620262
# #     c6 = -0.7483574222662208128858646256056442436988142359531261888282163941214837199890769
# #     c7 = 0.1238457029355639890633309483614995366456841215797361229753033357177669422002932
# #     c8 = -0.009088908104546665130978765013488165959525861872199676157403513022150601871065725
# #     vfmadd(vfmadd(vfmadd(vfmadd(vfmadd(vfmadd(vfmadd(vfmadd(c8, x, c7), x, c6), x, c5), x, c4), x, c3), x, c2), x, c1), x, c0)
# # end

# @inline function log12_9(x)
#     c0 = -3.585298173957385989978703196069874884252520373044145691918937267709490632531495
#     c1 = 9.175512833125597682830796241005729674265545403887738180485766461445305540243307
#     c2 = -12.88095295153181316823685136857741572991739006443899617826355360043323954685944
#     c3 = 13.97023726212960122584115270603499243664357604160987633843608122694259681596965
#     c4 = -10.88490936111543295028953131289531668167628753782395329546272132276910097503374
#     c5 = 5.991068536274933212572646135712744769742009519680374067042945152617478421745127
#     c6 = -2.274948730892808615344167059260983651564370908737772545840737790860551065669143
#     c7 = 0.5675102681955644633507049003267428212647262631800341022730034095800899973034826
#     c8 = -0.0837672704777144218198482121635843275577262111753625917138697703746721702699812
#     c9 = 0.005547594786690080999583023292516375924224904367437611790985685156038375417640469
#     vfmadd(vfmadd(vfmadd(vfmadd(vfmadd(vfmadd(vfmadd(vfmadd(vfmadd(c9, x, c8), x, c7), x, c6), x, c5), x, c4), x, c3), x, c2), x, c1), x, c0)
# end

# @inline function log12_16(x)
#     c0 = -4.37978460046975823625062481303177042964655415404536866782931434738284915486727
#     c1 = 16.31871050954055333956396222100764137104164259191322821181471405491328356924757
#     c2 = -43.09849769930020991569153712762384436449878854640277987997630267083267961690192
#     c3 = 94.07576064264437747981356599549555160319548024066792684545068563206644131393042
#     c4 = -160.2802910490630628148404579429898607373312993988526153520977924333905240690591
#     c5 = 214.2893836574384742952624060466181446104100940691092327871064519814161128088174
#     c6 = -227.1155113171482030403928820851457771576391114870681368560167920657167422599979
#     c7 = 192.2033142217916816132668362174888403726852676823685986312496157227978471900725
#     c8 = -130.2749105774541912642407178793062852794871781993373055603363370412151378757816
#     c9 = 70.61291078124372952080517981043356347605704563698401990419283029704191724347773
#     c10 = -30.40547369086585240971017001009032050577317283007958566465348067875074495409737
#     c11 = 10.26725226239003763659962692748550577872238983137251963601200333368657729640871
#     c12 = -2.660799790256174577955863634381340009240338534864386807650530364213169177140881
#     c13 = 0.5109263929851040905835535339027051835170022363218384968860233384732631558306201
#     c14 = -0.0684866374459784970521584206556922390960261357021144271918860042253780684140548
#     c15 = 0.005721132779491026341043740576624751586837718238326317641335022781958996821344917
#     c16 = -0.0002242388100013828616223303469764320113892766876681769892427408765524456008665568
# end

# @inline function log12_5_5(x)
#     n0 = -6.109900562053389599719353325364255967820553763613149427966199577373249418751788
#     n1 = -47.21033890517859950528344185143337802180242634918153480624598815718145572048678
#     n2 = -23.72899742183309520410030528859144026073452004231769448510028309291765826686717
#     n3 = 52.7756604934258744605387726053489962661521090040028567523770922726855897228027
#     n4 = 23.016711670393619896797584282375610751405740717936585553335237271483200924596
#     n5 = 1.256864725247499354645915543743422456498233774923938832688724623689696616505147

#     d0 = 1.0
#     d1 = 17.89075270582490509520473828435679777182649663126891381927673222776567771555027
#     d2 = 50.90705741116772592703914300524183279313383357917248215823723397522708772511337
#     d3 = 35.99672550568958966594388544740561745244861045550134035742275301908396910863221
#     d4 = 6.325336279410182109220032402097269110224428178264585067800799249731808867726456
#     d5 = 0.1767766952966368811002110905496910114147496513674936569539785785207122898660294

#     n = vfmadd(n5, x, n4)
#     d = vfmadd(d5, x, d4)
#     # n = vfmadd(n, x, n4)
#     # d = vfmadd(d, x, d4)
#     n = vfmadd(n, x, n3)
#     d = vfmadd(d, x, d3)
#     n = vfmadd(n, x, n2)
#     d = vfmadd(d, x, d2)
#     n = vfmadd(n, x, n1)
#     d = vfmadd(d, x, d1)
#     n = vfmadd(n, x, n0)
#     d = vfmadd(d, x, d0)
#     n / d
# end

# unlikely to do anything, but avoids bias when more than 12 of the leading bits are 0
# @inline function shift_excess_zeros(u::AbstractSIMD{W,UInt64}, lz::AbstractSIMD{W,UInt64}) where {W}
@inline function shift_excess_zeros(u, lz)
    lzsub = reinterpret(Int64, lz) - 11
    bitmask = lzsub > 0
    lzshift = reinterpret(UInt64, lzsub)
    ifelse(bitmask, ( u << lzshift ), u)
end
# @inline function shift_excess_zerosv2(u::AbstractSIMD{W,UInt64}, lz::AbstractSIMD{W,UInt64}) where {W}
#     lzsub = reinterpret(Int64, lz) - 11
#     ( u << reinterpret(UInt64, max( vzero(Vec{W,Int64}), lzsub ) ))
# end

# @inline function nlog01v2(u::AbstractSIMD{W,UInt64}, ::Type{Float64}) where {W}
#     lz = leading_zeros( u )
#     # f = mask(u, Float64) # shift by lz
#     f = mask(shift_excess_zeros(u, lz), Float64) # shift by lz
#     # l2h = log12_9(f)
#     l2h = log12_5_5(f)
#     l2 = l2h - (lz + 1)
#     -0.6931471805599453 * l2
# end

@inline function log2_3q(v, e)
  T = eltype(v)
  m1 = v * v
  fma1 = muladd(m1, T(0.22119417504560815), T(0.22007686931522777))
  fma2 = muladd(fma1, m1, T(0.26237080574885147))
  fma3 = muladd(fma2, m1, T(0.32059774779444955))
  fma4 = muladd(fma3, m1, T(0.41219859454853247))
  fma5 = muladd(fma4, m1, T(0.5770780162997059))
  fma6 = muladd(fma5, m1, T(0.9617966939260809))
  m2 = v * T(2.8853900817779268)
  fma7 = VectorizationBase.vfmsub(v, T(2.8853900817779268), m2)
  a1 = e + m2
  s1 = e - a1
  a2 = m2 + s1
  a3 = fma7 + a2
  m3 = v * m1
  a4 = a1 + a3
  muladd(fma6, m3, a4)
end
@inline function nlog01(u, ::Type{T}) where {T}
  lz = reinterpret(Base.uinttype(T), leading_zeros( u ))
  f = floatbitmask(shift_excess_zeros(u, lz), T) # shift by lz
  f = ( f - T(1.3333333333333333) ) / ( f + T(1.3333333333333333) )
  # l2h = log12_9(f)
  l2 = log2_3q(f, T(-0.5849625007211561814537389439478165087598144076924810604557526545410982277943579) - lz)
  T(-0.6931471805599453) * l2
end
# TODO: Add support for Float32

# @inline function nlog01(u::AbstractSIMD{W,UInt64}, ::Type{Float32}) where {W}
#     -log(mask(u, Float32) - oneopenconst(Float32))
# end


text:
using CUDA, Glob, MPI, SeisNoise, Dates, Serialization

# initialize MPI with CUDA
# suggested to use same number of CPU as GPU
MPI.Init()
comm = MPI.COMM_WORLD
lcomm = MPI.Comm_split_type(comm, MPI.MPI_COMM_TYPE_SHARED, MPI.Comm_rank(comm))
CUDA.device!(MPI.Comm_rank(lcomm) % length(devices()))
rank = MPI.Comm_rank(comm)
sz = MPI.Comm_size(comm)

# directory structure
FFTDIR = "/n/holyscratch01/denolle_lab/tclements/LASSO/FFT/"
CORRDIR = "/n/holyscratch01/denolle_lab/tclements/LASSO/CORR/"

# parameters
maxlag = 200.
Nper = 73

# files and such
files = glob("*",FFTDIR)
stas = [split(basename(s),'.')[2] |> (y -> parse(Int,y)) for s in files]
ind = sortperm(stas)
files = files[ind]
Nfiles = length(files)
Ntot = sz * Nper
splits = ceil(Int,Nfiles / Ntot)

# make output directories
if rank == 0
    for sta in stas
         mkpath(joinpath(CORRDIR,"2A.$sta..DPZ"))
     end
 end

MPI.Barrier(comm)

include("MPI-CC-impl.jl")
function main(files,splits,Nper,CORRDIR)

    comm = MPI.COMM_WORLD
    rank = MPI.Comm_rank(comm)
    sz = MPI.Comm_size(comm)
    Nfiles = length(files)
    Ntot = sz * Nper

    # loop through each split
    for ii = 1:splits
        if ii == splits
            extra = Nfiles % Ntot
            Nlocal = ceil(Int,extra / sz)
        else
            Nlocal = Nper
        end

        # get files for this rank
        startind = (ii - 1) * Ntot + rank * Nlocal + 1
        endind = (ii - 1) * Ntot + (rank + 1) * Nlocal
        endind = min(endind,Nfiles)
        rankfiles = files[startind:endind]

	println("Loading split $ii on rank $rank $(now())")
        # load files on rank
        arr1 = loadFFT(rankfiles) .|> gpu

        # correlate all files on this rank
        data = Channel{AbstractArray}(1)
	println("Cross-correlating split $ii on rank $rank $(now())")
        @sync begin
            @async one2all(arr1,maxlag,CORRDIR)
            @async begin
                src = mod(rank-1,sz)
                startind = (ii - 1) * Ntot + src * Nlocal + 1
                endind = (ii - 1) * Ntot + (src + 1) * Nlocal
                endind = min(endind,Nfiles)
                srcfiles = files[startind:endind]
                put!(data,loadFFT(srcfiles))
            end
        end
        arr2 = take!(data) .|> gpu
        MPI.Barrier(comm)

        # loop through splits
        for jj = 2:sz
			data = Channel{AbstractArray}(1)
	    println("Cross-correlating $ii-$jj on rank $rank $(now())")
            @sync begin
                @async all2all(arr1,arr2,maxlag,CORRDIR)
                @async begin
                    src = mod(rank-jj,sz)
                    startind = (ii - 1) * Ntot + src * Nlocal + 1
                    endind = (ii - 1) * Ntot + (src + 1) * Nlocal
                    endind = min(endind,Nfiles)
                    srcfiles = files[startind:endind]
                    put!(data,loadFFT(srcfiles))
                end
            end
            # transfer CPU arr to gpu
			freemem!(arr2)
			arr2 = take!(data) .|> gpu
        end
		freemem!(arr2)

        # loop through each remaining split
        for jj = ii + 1:splits

            if jj == splits
                extra = Nfiles % Ntot
                Nlocal = ceil(Int,extra / sz)
            else
                Nlocal = Nper
            end

            # get files for this rank
            startind = (jj - 1) * Ntot + rank * Nlocal + 1
            endind = (jj - 1) * Ntot + (rank + 1) * Nlocal
            endind = min(endind,Nfiles)
            rankfiles = files[startind:endind]

			# load onto GPU
	    println("Loading split $jj on rank $rank $(now())")
            arr2 = loadFFT(rankfiles) .|> gpu

            for kk = 2:sz
				data = Channel{AbstractArray}(1)
				println("Cross-correlating $ii-$jj on rank $rank $(now())")
				@sync begin
	            	@async all2all(arr1,arr2,maxlag,CORRDIR)
					@async begin
	                    src = mod(rank-kk,sz)
	                    startind = (jj - 1) * Ntot + src * Nlocal + 1
	                    endind = (jj - 1) * Ntot + (src + 1) * Nlocal
	                    endind = min(endind,Nfiles)
	                    srcfiles = files[startind:endind]
	                    put!(data,loadFFT(srcfiles))
	                end
				end
				freemem!(arr2)
				arr2 = take!(data) .|> gpu

            end
			freemem!(arr2)
            MPI.Barrier(comm)
        end
        MPI.Barrier(comm)
    end
    return nothing
end

# run test first
main(files[1:16],4,2,CORRDIR)
if rank == 0
	testfiles = glob("*/*",CORRDIR)
	rm.(testfiles)
end
MPI.Barrier(comm)
t1 = now()
main(files,splits,Nper,CORRDIR)
t2 = now()
total = Dates.canonicalize(Dates.CompoundPeriod(t2-t1))
println("Total computation took $total")
open(expanduser("~/LASSO-MPI.txt"),"w") do io
    println(io,"Total computation took $total")
end
MPI.Finalize()


text:
@cache struct EulerCache{uType,rateType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  tmp::uType
  k::rateType
  fsalfirst::rateType
end

@cache struct SplitEulerCache{uType,rateType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  tmp::uType
  k::rateType
  fsalfirst::rateType
end

function alg_cache(alg::SplitEuler,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  SplitEulerCache(u,uprev,similar(u),zero(rate_prototype),zero(rate_prototype))
end

struct SplitEulerConstantCache <: OrdinaryDiffEqConstantCache end

alg_cache(alg::SplitEuler,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = SplitEulerConstantCache()

function alg_cache(alg::Euler,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  EulerCache(u,uprev,similar(u),zero(rate_prototype),zero(rate_prototype))
end

struct EulerConstantCache <: OrdinaryDiffEqConstantCache end

alg_cache(alg::Euler,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = EulerConstantCache()

@cache struct HeunCache{uType,rateType,uNoUnitsType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  tmp::uType
  atmp::uNoUnitsType
  k::rateType
  fsalfirst::rateType
end

@cache struct RalstonCache{uType,rateType,uNoUnitsType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  tmp::uType
  atmp::uNoUnitsType
  k::rateType
  fsalfirst::rateType
end

function alg_cache(alg::Heun,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  HeunCache(u,uprev,similar(u),similar(u, uEltypeNoUnits),zero(rate_prototype),
            zero(rate_prototype))
end

function alg_cache(alg::Ralston,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  RalstonCache(u,uprev,similar(u),similar(u, uEltypeNoUnits),zero(rate_prototype),
               zero(rate_prototype))
end

struct HeunConstantCache <: OrdinaryDiffEqConstantCache end

alg_cache(alg::Heun,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = HeunConstantCache()

struct RalstonConstantCache <: OrdinaryDiffEqConstantCache end

alg_cache(alg::Ralston,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = RalstonConstantCache()

@cache struct MidpointCache{uType,rateType,uNoUnitsType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k::rateType
  tmp::uType
  atmp::uNoUnitsType
  fsalfirst::rateType
end

struct MidpointConstantCache <: OrdinaryDiffEqConstantCache end

function alg_cache(alg::Midpoint,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tmp = similar(u); atmp = similar(u, uEltypeNoUnits)
  k = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  MidpointCache(u,uprev,k,tmp,atmp,fsalfirst)
end

alg_cache(alg::Midpoint,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = MidpointConstantCache()

@cache struct RK4Cache{uType,rateType,uNoUnitsType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  fsalfirst::rateType
  k₂::rateType
  k₃::rateType
  k₄::rateType
  k::rateType
  tmp::uType
  atmp::uNoUnitsType
end

struct RK4ConstantCache <: OrdinaryDiffEqConstantCache end

function alg_cache(alg::RK4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  k₁ = zero(rate_prototype)
  k₂ = zero(rate_prototype)
  k₃ = zero(rate_prototype)
  k₄ = zero(rate_prototype)
  k  = zero(rate_prototype)
  tmp = similar(u); atmp = similar(u, uEltypeNoUnits)
  RK4Cache(u,uprev,k₁,k₂,k₃,k₄,k,tmp,atmp)
end

alg_cache(alg::RK4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = RK4ConstantCache()

@cache struct CarpenterKennedy2N54Cache{uType,rateType,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k::rateType
  tmp::uType
  fsalfirst::rateType
  tab::TabType
end

struct CarpenterKennedy2N54ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache
  A2::T
  A3::T
  A4::T
  A5::T
  B1::T
  B2::T
  B3::T
  B4::T
  B5::T
  c2::T2
  c3::T2
  c4::T2
  c5::T2

  function CarpenterKennedy2N54ConstantCache(::Type{T}, ::Type{T2}) where {T,T2}
    A2 = T(-567301805773/1357537059087)
    A3 = T(-2404267990393/2016746695238)
    A4 = T(-3550918686646/2091501179385)
    A5 = T(-1275806237668/842570457699)
    B1 = T(1432997174477/9575080441755)
    B2 = T(5161836677717/13612068292357)
    B3 = T(1720146321549/2090206949498)
    B4 = T(3134564353537/4481467310338)
    B5 = T(2277821191437/14882151754819)
    c2 = T2(1432997174477/9575080441755)
    c3 = T2(2526269341429/6820363962896)
    c4 = T2(2006345519317/3224310063776)
    c5 = T2(2802321613138/2924317926251)
    new{T,T2}(A2, A3, A4, A5, B1, B2, B3, B4, B5, c2, c3, c4, c5)
  end
end

function alg_cache(alg::CarpenterKennedy2N54,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tmp = similar(u)
  k = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  tab = CarpenterKennedy2N54ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  CarpenterKennedy2N54Cache(u,uprev,k,tmp,fsalfirst,tab)
end

function alg_cache(alg::CarpenterKennedy2N54,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  CarpenterKennedy2N54ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
end

@cache struct BS3Cache{uType,rateType,uNoUnitsType,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  fsalfirst::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  utilde::uType
  tmp::uType
  atmp::uNoUnitsType
  tab::TabType
end

function alg_cache(alg::BS3,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tab = BS3ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  k1 = zero(rate_prototype)
  k2 = zero(rate_prototype)
  k3 = zero(rate_prototype)
  k4 = zero(rate_prototype)
  utilde = similar(u)
  atmp = similar(u,uEltypeNoUnits)
  tmp = similar(u)
  BS3Cache(u,uprev,k1,k2,k3,k4,utilde,tmp,atmp,tab)
end

alg_cache(alg::BS3,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = BS3ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))

@cache struct OwrenZen3Cache{uType,rateType,uNoUnitsType,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k1::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  utilde::uType
  tmp::uType
  atmp::uNoUnitsType
  tab::TabType
end

function alg_cache(alg::OwrenZen3,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tab = OwrenZen3ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  k1 = zero(rate_prototype)
  k2 = zero(rate_prototype)
  k3 = zero(rate_prototype)
  k4 = zero(rate_prototype)
  utilde = similar(u)
  atmp = similar(u,uEltypeNoUnits)
  tmp = similar(u)
  OwrenZen3Cache(u,uprev,k1,k2,k3,k4,utilde,tmp,atmp,tab)
end

alg_cache(alg::OwrenZen3,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = OwrenZen3ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))

@cache struct OwrenZen4Cache{uType,rateType,uNoUnitsType,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k1::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  k5::rateType
  k6::rateType
  utilde::uType
  tmp::uType
  atmp::uNoUnitsType
  tab::TabType
end

function alg_cache(alg::OwrenZen4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tab = OwrenZen4ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  k1 = zero(rate_prototype)
  k2 = zero(rate_prototype)
  k3 = zero(rate_prototype)
  k4 = zero(rate_prototype)
  k5 = zero(rate_prototype)
  k6 = zero(rate_prototype)
  utilde = similar(u)
  atmp = similar(u,uEltypeNoUnits)
  tmp = similar(u)
  OwrenZen4Cache(u,uprev,k1,k2,k3,k4,k5,k6,utilde,tmp,atmp,tab)
end

alg_cache(alg::OwrenZen4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = OwrenZen4ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))

@cache struct OwrenZen5Cache{uType,rateType,uNoUnitsType,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k1::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  k5::rateType
  k6::rateType
  k7::rateType
  k8::rateType
  utilde::uType
  tmp::uType
  atmp::uNoUnitsType
  tab::TabType
end

function alg_cache(alg::OwrenZen5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tab = OwrenZen5ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  k1 = zero(rate_prototype)
  k2 = zero(rate_prototype)
  k3 = zero(rate_prototype)
  k4 = zero(rate_prototype)
  k5 = zero(rate_prototype)
  k6 = zero(rate_prototype)
  k7 = zero(rate_prototype)
  k8 = zero(rate_prototype)
  utilde = similar(u)
  atmp = similar(u,uEltypeNoUnits)
  tmp = similar(u)
  OwrenZen5Cache(u,uprev,k1,k2,k3,k4,k5,k6,k7,k8,utilde,tmp,atmp,tab)
end

alg_cache(alg::OwrenZen5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = OwrenZen5ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))

@cache struct BS5Cache{uType,rateType,uNoUnitsType,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k1::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  k5::rateType
  k6::rateType
  k7::rateType
  k8::rateType
  utilde::uType
  tmp::uType
  atmp::uNoUnitsType
  tab::TabType
end

function alg_cache(alg::BS5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tab = BS5ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  k1 = zero(rate_prototype)
  k2 = zero(rate_prototype)
  k3 = zero(rate_prototype)
  k4 = zero(rate_prototype)
  k5 = zero(rate_prototype)
  k6 = zero(rate_prototype)
  k7 = zero(rate_prototype)
  k8 = zero(rate_prototype)
  utilde = similar(u)
  atmp = similar(u,uEltypeNoUnits)
  tmp = similar(u)
  BS5Cache(u,uprev,k1,k2,k3,k4,k5,k6,k7,k8,utilde,tmp,atmp,tab)
end

alg_cache(alg::BS5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = BS5ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))

@cache struct Tsit5Cache{uType,rateType,uNoUnitsType,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k1::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  k5::rateType
  k6::rateType
  k7::rateType
  utilde::uType
  tmp::uType
  atmp::uNoUnitsType
  tab::TabType
end

function alg_cache(alg::Tsit5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tab = Tsit5ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  k1 = zero(rate_prototype)
  k2 = zero(rate_prototype)
  k3 = zero(rate_prototype)
  k4 = zero(rate_prototype)
  k5 = zero(rate_prototype)
  k6 = zero(rate_prototype)
  k7 = zero(rate_prototype)
  utilde = similar(u)
  atmp = similar(u,uEltypeNoUnits)
  tmp = similar(u)
  Tsit5Cache(u,uprev,k1,k2,k3,k4,k5,k6,k7,utilde,tmp,atmp,tab)
end

alg_cache(alg::Tsit5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = Tsit5ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))

@cache struct DP5Cache{uType,rateType,uNoUnitsType,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k1::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  k5::rateType
  k6::rateType
  k7::rateType
  dense_tmp3::rateType
  dense_tmp4::rateType
  update::rateType
  bspl::rateType
  utilde::uType
  tmp::uType
  atmp::uNoUnitsType
  tab::TabType
end

function alg_cache(alg::DP5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  k1 = zero(rate_prototype)
  k2 = zero(rate_prototype)
  k3 = zero(rate_prototype)
  k4 = zero(rate_prototype)
  k5 = zero(rate_prototype)
  k6 = k2
  k7 = zero(rate_prototype) # This is FSAL'd to k1
  dense_tmp3 = k2
  dense_tmp4 = k5
  bspl = k3

  tmp = similar(u) # has to be separate for FSAL
  utilde = tmp

  if eltype(u) != uEltypeNoUnits || calck
    update = zero(rate_prototype)
    atmp = similar(u,uEltypeNoUnits)
  else
    update = k7
    atmp = k3
  end

  tab = DP5ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  cache = DP5Cache(u,uprev,k1,k2,k3,k4,k5,k6,k7,dense_tmp3,dense_tmp4,update,bspl,utilde,tmp,atmp,tab)
  cache
end

alg_cache(alg::DP5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = DP5ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))

@cache struct DP5ThreadedCache{uType,rateType,uNoUnitsType,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k1::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  k5::rateType
  k6::rateType
  k7::rateType
  dense_tmp3::rateType
  dense_tmp4::rateType
  update::rateType
  bspl::rateType
  utilde::uType
  tmp::uType
  atmp::uNoUnitsType
  tab::TabType
end

function alg_cache(alg::DP5Threaded,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  k1 = zero(rate_prototype)
  k2 = zero(rate_prototype)
  k3 = zero(rate_prototype)
  k4 = zero(rate_prototype)
  k5 = zero(rate_prototype)
  k6 = zero(rate_prototype)
  k7 = zero(rate_prototype)
  dense_tmp3 = zero(rate_prototype)
  dense_tmp4 = zero(rate_prototype)
  update = zero(rate_prototype)
  bspl = zero(rate_prototype)
  utilde = similar(u)
  tmp = similar(u); atmp = similar(u,uEltypeNoUnits)
  tab = DP5ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  DP5ThreadedCache(u,uprev,k1,k2,k3,k4,k5,k6,k7,dense_tmp3,dense_tmp4,update,bspl,utilde,tmp,atmp,tab)
end

@cache struct Anas5Cache{uType,rateType,uNoUnitsType,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k1::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  k5::rateType
  k6::rateType
  k7::rateType
  utilde::uType
  tmp::uType
  atmp::uNoUnitsType
  tab::TabType
end

function alg_cache(alg::Anas5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tab = Anas5ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  k1 = zero(rate_prototype)
  k2 = zero(rate_prototype)
  k3 = zero(rate_prototype)
  k4 = zero(rate_prototype)
  k5 = zero(rate_prototype)
  k6 = zero(rate_prototype)
  k7 = zero(rate_prototype)
  utilde = similar(u)
  atmp = similar(u,uEltypeNoUnits)
  tmp = similar(u)
  Anas5Cache(u,uprev,k1,k2,k3,k4,k5,k6,k7,utilde,tmp,atmp,tab)
end

alg_cache(alg::Anas5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = Anas5ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))


text:
### A Pluto.jl notebook ###
# v0.17.3

using Markdown
using InteractiveUtils

# ╔═╡ 581b00ef-e590-4b33-910c-f4f9525a9626
# hideall

using Genie

# ╔═╡ f3799f4a-4926-11ec-2bf8-a5ce535d651e
md"""
# Developing Genie Web Services

Starting up ad-hoc web servers at the REPL and writing small scripts to wrap micro-services works great, but production apps tend to become complex very quickly. They also have more stringent requirements, like managing dependencies, compressing assets, reloading code, logging, environments, or structuring the codebase in a way which promotes efficient workflows when working in teams.

Genie apps provide all these features, from dependency management and versioning (using Julia's `Pkg`, since a Genie app is a Julia project), to a powerful asset pipeline (using industry vetted tools like Yarn and Webpack), automatic code reloading in development (provided by `Revise.jl`), and a clear resource-oriented MVC layout.

Genie enables a modular approach towards app building, allowing to add more components as the need arises. You can start with the web service template (which includes dependencies management, logging, environments, and routing), and grow it by sequentially adding DB persistence  (through the SearchLight ORM), high performance HTML view templates with embedded Julia (via Flax), asset pipeline and compilation, and more.

## Setting up a Genie Web Service project

Genie packs handy generator features and templates which help bootstrapping and setting up various parts of an application. For bootstrapping a new app we need to invoke one of the functions in the `newapp` family:
"""

# ╔═╡ f1e11f9d-a5ce-4d6b-9e26-ce39800cac20
md"""

#### Adding Genie package

```julia
julia> using Pkg

julia> Pkg.add("Genie")
```

##### Importing Genie

```julia
julia> using Genie

julia> Genie.newapp_webservice("MyGenieApp")
```

"""

# ╔═╡ 4226f596-ce42-4953-966f-ea1afe32b698
# hideall

Genie.newapp_webservice("MyGenieApp")

# ╔═╡ 2ebe66be-24b5-4765-8a58-0ccac149799b
md"""
If you follow the log messages in the REPL you will see that the command will trigger a flurry of actions in order to set up the new project:

- it creates a new folder, `MyGenieApp/`, which will hosts the files of the app and whose name corresponds to the name of the app,
- within the `MyGenieApp/` folder, it creates the files and folders needed by the app,
- changes the active directory to `MyGenieApp/` and creates a new Julia project within it (adding the `Project.toml` file),
- installs all the required dependencies for the new Genie app (using `Pkg` and the standard `Manifest.toml` file), and finally,
- starts the web server

---
**TIP**

Check out the inline help for `Genie.newapp`, `Genie.newapp_webservice`, `Genie.newapp_mvc`, and `Genie.newapp_fullstack` too see what options are available for bootstrapping applications. We'll go over the different configurations in upcoming sections.

---

## The file structure

Our newly created web service has this file structure:

```julia
├── Manifest.toml
├── Project.toml
├── bin
├── bootstrap.jl
├── config
├── genie.jl
├── log
├── public
├── routes.jl
└── src
```

These are the roles of each of the files and folders:

- `Manifest.toml` and `Project.toml` are used by Julia and `Pkg` to manage the app's dependencies.
- `bin/` includes scripts for starting up a Genie REPL or a Genie server.
- `bootstrap.jl`, `genie.jl`, as well as all the files within `src/` are used by Genie to load the application and _should not be user modified_.
- `config/` includes the per-environment configuration files.
- `log/` is used by Genie to store per-environment log files.
- `public/` is the document root, which includes static files exposed by the app on the network/internet.
- `routes.jl` is the dedicated file for registering Genie routes.

---
**HEADS UP**

After creating a new app you might need to change the file permissions to allow editing/saving the files such as `routes.jl`.

---

## Adding logic

You can now edit the `routes.jl` file to add some logic, at the bottom of the file:"""

# ╔═╡ 8b7366bc-beca-428e-8e03-ae02e5000e8f
md"""
```julia
route("/hello") do
  "Welcome to Genie!"
end
```
"""

# ╔═╡ 685dec36-b77a-470f-aa59-0e460aaf1188
md"""
If you now visit <http://127.0.0.1:8000/hello> you'll see a warm greeting.
"""

# ╔═╡ 76c662cc-74f9-4f55-9c55-9fdd14febc88
md"""
## Growing the app

Genie apps are just plain Julia projects. This means that `routes.jl` will behave like any other Julia script - you can reference extra packages, you can switch into `pkg>` mode to manage per project dependencies, include other files, etcetera.

If you have existing Julia code that you want to quickly load into a Genie app, you can add a `lib/` folder in the root of the app and place your Julia files there. When available, `lib/` and all its subfolders are automatically added to the `LOAD_PATH` by Genie, recursively.

If you need to add database support, you can always add the SearchLight ORM by using the dedicated generator, running `julia> Genie.Generator.db_support()` in the app's REPL.

However, if your app grows in complexity and you develop it from scratch, it is more efficient to take advantage of Genie's resource-oriented MVC structure.
"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
Genie = "c43c736e-a2d1-11e8-161f-af95117fbd1e"

[compat]
Genie = "~4.12.0"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[ArgParse]]
deps = ["Logging", "TextWrap"]
git-tree-sha1 = "3102bce13da501c9104df33549f511cd25264d7d"
uuid = "c7e460c6-2fb9-53a9-8c5b-16f535851c63"
version = "1.1.4"

[[ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[CSTParser]]
deps = ["Tokenize"]
git-tree-sha1 = "6cc1759204bed5a4e2a5c2f00901fd5d90bc7a62"
uuid = "00ebfdb7-1f24-5e51-bd34-a7502290713f"
version = "3.3.1"

[[CodeTracking]]
deps = ["InteractiveUtils", "UUIDs"]
git-tree-sha1 = "9aa8a5ebb6b5bf469a7e0e2b5202cf6f8c291104"
uuid = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2"
version = "1.0.6"

[[CommonMark]]
deps = ["Crayons", "JSON", "URIs"]
git-tree-sha1 = "4cd7063c9bdebdbd55ede1af70f3c2f48fab4215"
uuid = "a80b9123-70ca-4bc0-993e-6e3bcb318db6"
version = "0.8.6"

[[Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "44c37b4636bc54afac5c574d2d02b625349d6582"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.41.0"

[[Crayons]]
git-tree-sha1 = "249fe38abf76d48563e2f4556bebd215aa317e15"
uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
version = "4.1.1"

[[DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "3daef5523dd2e769dad2365274f760ff5f282c7d"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.11"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[EzXML]]
deps = ["Printf", "XML2_jll"]
git-tree-sha1 = "0fa3b52a04a4e210aeb1626def9c90df3ae65268"
uuid = "8f5d6c58-4d21-5cfd-889c-e3ad7ee6a615"
version = "1.1.0"

[[FilePathsBase]]
deps = ["Compat", "Dates", "Mmap", "Printf", "Test", "UUIDs"]
git-tree-sha1 = "04d13bfa8ef11720c24e4d840c0033d145537df7"
uuid = "48062228-2e41-5def-b9a4-89aafe57970f"
version = "0.9.17"

[[FileWatching]]
uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee"

[[GMP_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "781609d7-10c4-51f6-84f2-b8444358ff6d"

[[Genie]]
deps = ["ArgParse", "Dates", "Distributed", "EzXML", "FilePathsBase", "HTTP", "HttpCommon", "Inflector", "JSON3", "JuliaFormatter", "Logging", "Markdown", "MbedTLS", "Millboard", "Nettle", "OrderedCollections", "Pkg", "REPL", "Random", "Reexport", "Revise", "SHA", "Serialization", "Sockets", "UUIDs", "Unicode", "VersionCheck", "YAML"]
git-tree-sha1 = "226ea9a883fc5d2112344fa3ac6c28ee335d2cd3"
uuid = "c43c736e-a2d1-11e8-161f-af95117fbd1e"
version = "4.12.0"

[[HTTP]]
deps = ["Base64", "Dates", "IniFile", "Logging", "MbedTLS", "NetworkOptions", "Sockets", "URIs"]
git-tree-sha1 = "0fa77022fe4b511826b39c894c90daf5fce3334a"
uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3"
version = "0.9.17"

[[HttpCommon]]
deps = ["Dates", "Nullables", "Test", "URIParser"]
git-tree-sha1 = "46313284237aa6ca67a6bce6d6fbd323d19cff59"
uuid = "77172c1b-203f-54ac-aa54-3f1198fe9f90"
version = "0.5.0"

[[Inflector]]
deps = ["Unicode"]
git-tree-sha1 = "8555b54ddf27806b070ce1d1cf623e1feb13750c"
uuid = "6d011eab-0732-4556-8808-e463c76bf3b6"
version = "1.0.1"

[[IniFile]]
deps = ["Test"]
git-tree-sha1 = "098e4d2c533924c921f9f9847274f2ad89e018b8"
uuid = "83e8ac13-25f8-5344-8a64-a9f2b223428f"
version = "0.5.0"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "abc9885a7ca2052a736a600f7fa66209f96506e1"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.4.1"

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "3c837543ddb02250ef42f4738347454f95079d4e"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.3"

[[JSON3]]
deps = ["Dates", "Mmap", "Parsers", "StructTypes", "UUIDs"]
git-tree-sha1 = "7d58534ffb62cd947950b3aa9b993e63307a6125"
uuid = "0f8b85d8-7281-11e9-16c2-39a750bddbf1"
version = "1.9.2"

[[JuliaFormatter]]
deps = ["CSTParser", "CommonMark", "DataStructures", "Pkg", "Tokenize"]
git-tree-sha1 = "fcfaddc61f766211b2c835d3eceaf999b6ea9555"
uuid = "98e50ef6-434e-11e9-1051-2b60c6c9e899"
version = "0.22.4"

[[JuliaInterpreter]]
deps = ["CodeTracking", "InteractiveUtils", "Random", "UUIDs"]
git-tree-sha1 = "b55aae9a2bf436fc797d9c253a900913e0e90178"
uuid = "aa1ae85d-cabe-5617-a682-6adf51b2e16a"
version = "0.9.3"

[[LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[Libiconv_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "42b62845d70a619f063a7da093d995ec8e15e778"
uuid = "94ce4f54-9a6c-5748-9c1c-f9c7231a4531"
version = "1.16.1+1"

[[LinearAlgebra]]
deps = ["Libdl"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[LoweredCodeUtils]]
deps = ["JuliaInterpreter"]
git-tree-sha1 = "6b0440822974cab904c8b14d79743565140567f6"
uuid = "6f1432cf-f94c-5a45-995e-cdbf5db27b0b"
version = "2.2.1"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[MbedTLS]]
deps = ["Dates", "MbedTLS_jll", "Random", "Sockets"]
git-tree-sha1 = "1c38e51c3d08ef2278062ebceade0e46cefc96fe"
uuid = "739be429-bea8-5141-9913-cc70e7f3736d"
version = "1.0.3"

[[MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[Millboard]]
git-tree-sha1 = "ea6a5b7e56e76d8051023faaa11d91d1d881dac3"
uuid = "39ec1447-df44-5f4c-beaa-866f30b4d3b2"
version = "0.2.5"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[Nettle]]
deps = ["Libdl", "Nettle_jll"]
git-tree-sha1 = "a68340b9edfd98d0ed96aee8137cb716ea3b6dea"
uuid = "49dea1ee-f6fa-5aa6-9a11-8816cee7d4b9"
version = "0.5.1"

[[Nettle_jll]]
deps = ["Artifacts", "GMP_jll", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "eca63e3847dad608cfa6a3329b95ef674c7160b4"
uuid = "4c82536e-c426-54e4-b420-14f461c4ed8b"
version = "3.7.2+0"

[[NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[Nullables]]
git-tree-sha1 = "8f87854cc8f3685a60689d8edecaa29d2251979b"
uuid = "4d1e1d77-625e-5b40-9113-a560ec7a8ecd"
version = "1.0.0"

[[OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "13468f237353112a01b2d6b32f3d0f80219944aa"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.2.2"

[[Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[Preferences]]
deps = ["TOML"]
git-tree-sha1 = "2cf929d64681236a2e074ffafb8d568733d2e6af"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.2.3"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[ProgressMeter]]
deps = ["Distributed", "Printf"]
git-tree-sha1 = "afadeba63d90ff223a6a48d2009434ecee2ec9e8"
uuid = "92933f4c-e287-5a05-a399-4b506db050ca"
version = "1.7.1"

[[REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.3.0"

[[Revise]]
deps = ["CodeTracking", "Distributed", "FileWatching", "JuliaInterpreter", "LibGit2", "LoweredCodeUtils", "OrderedCollections", "Pkg", "REPL", "Requires", "UUIDs", "Unicode"]
git-tree-sha1 = "2f9d4d6679b5f0394c52731db3794166f49d5131"
uuid = "295af30f-e4ad-537b-8983-00126c2a3abe"
version = "3.3.1"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[Scratch]]
deps = ["Dates"]
git-tree-sha1 = "0b4b7f1393cff97c33891da2a0bf69c6ed241fda"
uuid = "6c6a2e73-6563-6170-7368-637461726353"
version = "1.1.0"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[StringEncodings]]
deps = ["Libiconv_jll"]
git-tree-sha1 = "50ccd5ddb00d19392577902f0079267a72c5ab04"
uuid = "69024149-9ee7-55f6-a4c4-859efe599b68"
version = "0.3.5"

[[StructTypes]]
deps = ["Dates", "UUIDs"]
git-tree-sha1 = "d24a825a95a6d98c385001212dc9020d609f2d4f"
uuid = "856f2bd8-1eba-4b0a-8007-ebc267875bd4"
version = "1.8.1"

[[TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[TextWrap]]
git-tree-sha1 = "9250ef9b01b66667380cf3275b3f7488d0e25faf"
uuid = "b718987f-49a8-5099-9789-dcd902bef87d"
version = "1.0.1"

[[Tokenize]]
git-tree-sha1 = "0952c9cee34988092d73a5708780b3917166a0dd"
uuid = "0796e94c-ce3b-5d07-9a54-7f471281c624"
version = "0.5.21"

[[URIParser]]
deps = ["Unicode"]
git-tree-sha1 = "53a9f49546b8d2dd2e688d216421d050c9a31d0d"
uuid = "30578b45-9adc-5946-b283-645ec420af67"
version = "0.4.1"

[[URIs]]
git-tree-sha1 = "97bbe755a53fe859669cd907f2d96aee8d2c1355"
uuid = "5c2747f8-b7ea-4ff2-ba2e-563bfd36b1d4"
version = "1.3.0"

[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[UrlDownload]]
deps = ["HTTP", "ProgressMeter"]
git-tree-sha1 = "05f86730c7a53c9da603bd506a4fc9ad0851171c"
uuid = "856ac37a-3032-4c1c-9122-f86d88358c8b"
version = "1.0.0"

[[VersionCheck]]
deps = ["Dates", "JSON3", "Logging", "Pkg", "Random", "Scratch", "UrlDownload"]
git-tree-sha1 = "89ef2431dd59344ebaf052d0737205854ded0c62"
uuid = "a637dc6b-bca1-447e-a4fa-35264c9d0580"
version = "0.2.0"

[[XML2_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "Zlib_jll"]
git-tree-sha1 = "1acf5bdf07aa0907e0a37d3718bb88d4b687b74a"
uuid = "02c8fc9c-b97f-50b9-bbe4-9be30ff0a78a"
version = "2.9.12+0"

[[YAML]]
deps = ["Base64", "Dates", "Printf", "StringEncodings"]
git-tree-sha1 = "3c6e8b9f5cdaaa21340f841653942e1a6b6561e5"
uuid = "ddb6d928-2868-570f-bddf-ab3f9cf99eb6"
version = "0.4.7"

[[Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╟─f3799f4a-4926-11ec-2bf8-a5ce535d651e
# ╟─f1e11f9d-a5ce-4d6b-9e26-ce39800cac20
# ╠═581b00ef-e590-4b33-910c-f4f9525a9626
# ╠═4226f596-ce42-4953-966f-ea1afe32b698
# ╟─2ebe66be-24b5-4765-8a58-0ccac149799b
# ╟─8b7366bc-beca-428e-8e03-ae02e5000e8f
# ╟─685dec36-b77a-470f-aa59-0e460aaf1188
# ╟─76c662cc-74f9-4f55-9c55-9fdd14febc88
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
@testset "Longest Path" begin
    gfa_result = read_from_gfa("data/gfa_sample_1.gfa");
    res = find_graph_component(gfa_result);
    longest = find_longest_path(res.graph[1], Cbc.Optimizer, is_weighted = true);
    fasta = get_fasta(longest)
    longest_p_test = [7,2,3,4,5,6]
    longest_lp_test = ["4\t-","2\t+","3\t+","4\t+","5\t+","6\t+"]
    fasta_test = "TCACCCTGTATTAACTCCATCTTTGAGAAACATTTAATAATGTAATGTGTTTGTCATACAGGGTGAATACAGATGCAC"
    @test longest.path == longest_p_test
    @test longest.label_path == longest_lp_test
    @test fasta == fasta_test
    gfa_result_2 = read_from_gfa("data/gfa_sample_3.gfa");
    res_2 = find_graph_component(gfa_result_2);
    longest_2 = find_longest_path(res_2.graph[1], Cbc.Optimizer, is_weighted = false);
    fasta_2 = get_fasta(longest_2)
    longest_p_test_2 = [7,1,2,3,4,5,6]
    longest_lp_test_2 = ["4\t-","1\t+","2\t+","3\t+","4\t+","5\t+","6\t+"]
    fasta_test_2 = "TCACCCTGTAATTTAACTCCATCTTTGAGAAACATTTAATAATGTAATGTGTTTGTCATACAGGGTGAATACAGATGCAC"
    @test longest_2.path == longest_p_test_2
    @test longest_2.label_path == longest_lp_test_2
    @test fasta_2 == fasta_test_2
    gfa_result_3 = read_from_gfa("data/gfa_sample_4.gfa");
    res_3 = find_graph_component(gfa_result_3);
    longest_3 = find_longest_path(res_3.graph[1], Cbc.Optimizer, is_weighted = false);
    fasta_3 = get_fasta(longest_3)
    @test longest_3.path == longest_p_test_2
    @test longest_3.label_path == longest_lp_test_2
    @test fasta_3 == fasta_test_2
end


text:
module Metadata

if !isdefined(Core, :String)
    typealias String UTF8String
end

using FlatBuffers

@enum(Type_, BOOL = 0, INT8 = 1, INT16 = 2, INT32 = 3, INT64 = 4,
  UINT8 = 5, UINT16 = 6, UINT32 = 7, UINT64 = 8,
  FLOAT = 9, DOUBLE = 10,  UTF8 = 11,  BINARY = 12,
  CATEGORY = 13, TIMESTAMP = 14, DATE = 15, TIME = 16)

@enum(Encoding, PLAIN = 0, DICTIONARY = 1)

@enum(TimeUnit, SECOND = 0, MILLISECOND = 1, MICROSECOND = 2, NANOSECOND = 3)

# FlatBuffers.enumsizeof(::Type{TimeUnit}) = UInt8

type PrimitiveArray
    type_::Type_
    encoding::Encoding
    offset::Int64
    length::Int64
    null_count::Int64
    total_bytes::Int64
end

type CategoryMetadata
    levels::PrimitiveArray
    ordered::Bool
end

@default CategoryMetadata ordered=false

type TimestampMetadata
    unit::TimeUnit
    timezone::String
end

type DateMetadata
end

type TimeMetadata
    unit::TimeUnit
end

@union TypeMetadata Union{Void,CategoryMetadata,TimestampMetadata,DateMetadata,TimeMetadata}

type Column
    name::String
    values::PrimitiveArray
    metadata_type::Int8
    metadata::TypeMetadata
    user_metadata::String
end

function Column(name::String, values::PrimitiveArray, metadata::TypeMetadata=nothing, user_metadata::String="")
    return Column(name, values, FlatBuffers.typeorder(TypeMetadata, typeof(metadata)), metadata, user_metadata)
end

type CTable
    description::String
    num_rows::Int64
    columns::Vector{Column}
    version::Int32
    metadata::String
end

end # module

# wesm/feather/cpp/src/metadata_generated.h
# wesm/feather/cpp/src/types.h
const Type_2julia = Dict{Metadata.Type_,DataType}(
    Metadata.BOOL      => Arrow.Bool,
    Metadata.INT8      => Int8,
    Metadata.INT16     => Int16,
    Metadata.INT32     => Int32,
    Metadata.INT64     => Int64,
    Metadata.UINT8     => UInt8,
    Metadata.UINT16    => UInt16,
    Metadata.UINT32    => UInt32,
    Metadata.UINT64    => UInt64,
    Metadata.FLOAT     => Float32,
    Metadata.DOUBLE    => Float64,
    Metadata.UTF8      => WeakRefString{UInt8},
    Metadata.BINARY    => Vector{UInt8},
    Metadata.CATEGORY  => Int64,
    Metadata.TIMESTAMP => Int64,
    Metadata.DATE      => Int64,
    Metadata.TIME      => Int64
)

const julia2Type_ = Dict{DataType,Metadata.Type_}(
    Bool    => Metadata.BOOL,
    Int8    => Metadata.INT8,
    Int16   => Metadata.INT16,
    Int32   => Metadata.INT32,
    Int64   => Metadata.INT64,
    UInt8   => Metadata.UINT8,
    UInt16  => Metadata.UINT16,
    UInt32  => Metadata.UINT32,
    UInt64  => Metadata.UINT64,
    Float32 => Metadata.FLOAT,
    Float64 => Metadata.DOUBLE,
    String  => Metadata.UTF8,
    Vector{UInt8}   => Metadata.BINARY,
    DateTime   => Metadata.INT64,
    Date   => Metadata.INT32,
    WeakRefString{UInt8} => Metadata.UTF8
)

const NON_PRIMITIVE_TYPES = Set([Metadata.UTF8, Metadata.BINARY])
isprimitive(x::Metadata.Type_) = x in NON_PRIMITIVE_TYPES ? false : true

const TimeUnit2julia = Dict{Metadata.TimeUnit,DataType}(
    Metadata.SECOND => Arrow.Second,
    Metadata.MILLISECOND => Arrow.Millisecond,
    Metadata.MICROSECOND => Arrow.Microsecond,
    Metadata.NANOSECOND => Arrow.Nanosecond
)
const julia2TimeUnit = Dict{DataType,Metadata.TimeUnit}([(v, k) for (k,v) in TimeUnit2julia])


text:
using Base.Test
using POMDPToolbox

using LightDarkPOMDPs

pomdp = LightDark2D(init_dist=SymmetricNormal2([2.0, 2.0], 5.0))

policy = FunctionPolicy(s->-0.1*mean(s))

fnew = BootstrapFilter(pomdp, 100, rng=MersenneTwister(42))
ro = RolloutSimulator(rng=MersenneTwister(1), max_steps=10)
simulate(ro, pomdp, policy, fnew)

policy = FunctionPolicy(s->Vec2([0.0, 0.0]))
fnew = BasicParticleFilter(pomdp, LowVarianceResampler(1000), 1000, MersenneTwister(3))
is = Vec2([5.2, 0.0])
hr = HistoryRecorder(initial_state=is, max_steps=10, rng=MersenneTwister(3))
hist = simulate(hr, pomdp, policy, fnew)
@test norm(mean(last(belief_hist(hist)))-is) <= 0.5


text:
mutable struct Liability :> crAccount
end

max_stars_repo_name: adamwillisXanax/accountingJulia

max_issues_repo_name: adamwillisXanax/accountingJulia

max_forks_repo_name: adamwillisXanax/accountingJulia


text:
@testset "Wrapper" begin
    include("wrapper/exprs.jl")
    include("wrapper/naming_conventions.jl")
    include("wrapper/wrap.jl")
    include("wrapper/docs.jl")
end


text:
using Test, SparseArrays
import LuxurySparse: SparseMatrixCOO, allocated_coo, isdense

coo1 = SparseMatrixCOO([1,4,2,3,3,3], [1,1,2,4,3,4], [0.1, 0.2, 0.4im, 0.5, 0.3, 0.5im], 4, 4)

@testset "basic" begin
    @test eltype(coo1) == ComplexF64
    @test size(coo1, 1) == size(coo1, 2) == 4

    coo2 = allocated_coo(ComplexF64, 4, 4, 6)
    @test eltype(coo2) == ComplexF64
    @test nnz(coo2) == 6
    @test size(coo2) == (4, 4)
    @test size(coo2, 1) == size(coo2, 2) == 4

    @test coo1 == copy(coo1)
    @test coo1 == copyto!(coo2, coo1)

    @test Matrix(coo1) == [0.1 0 0 0; 0 0.4im 0 0; 0 0 0.3 0.5+0.5im; 0.2 0 0 0]
    @test sparse(coo1) == sparse(findnz(coo1)...)
    @test isdense(coo1) == false

    A = sprand(50,50, 0.2) |> Matrix
    @test SparseMatrixCOO(A) == A
end

@testset "sparse" begin
    p1 = copy(coo1)
    @test nonzeros(p1) == p1.vs
    p1.vs[2] = 0
    @test dropzeros!(p1) == [0.1 0 0 0; 0 0.4im 0 0; 0 0 0.3 0.5+0.5im; 0 0 0 0]
end

@testset "basicmath" begin
    @test coo1*2 == Matrix(coo1)*2
    @test 2*coo1 == 2*Matrix(coo1)
    @test coo1/2 == Matrix(coo1)/2
    @test coo1*2 isa SparseMatrixCOO
    @test coo1/2 isa SparseMatrixCOO
    @test 2*coo1 isa SparseMatrixCOO
end


text:
module WCSLIB
using Compat

export pscard,
       pvcard,
       wcsprm

export wcsbchk,
       #wcsbdx,
       wcsbth,
       wcscopy,
       wcserr_enable,
       wcserr_prt,
       wcsfree,
       wcshdo,
       #wcsidx,
       wcsini,
       #wcsmix,
       wcsmodify,
       wcsnps,
       wcsnpv,
       wcsp2s,
       wcsperr,
       wcspih,
       wcsprt,
       wcss2p,
       wcsset,
       #wcssptr,
       #wcssub,
       wcstab

import Base: convert

include("../deps/deps.jl")
include("libwcs_common.jl")
include("libwcs_h.jl")

if VERSION < v"0.3-"
    setfield!(a,b,c) = setfield(a,b,c)
    pscard(i,m,value) = pscard(convert(Cint,i),
                               convert(Cint,m),
                               convert(Array_72_Uint8,value))
    pvcard(i,m,value) = pvcard(convert(Cint,i),
                               convert(Cint,m),
                               convert(Cdouble,value))
end

function unsafe_store_vec!{T}(p::Ptr{T}, v::Vector{T})
    for i = 1:length(v)
        unsafe_store!(p, v[i], i)
    end
end

function convert(::Type{Array_72_Uint8}, s::ASCIIString)
    @assert length(s) < 72
    v = zeros(UInt8, 72)
    for i = 1:length(s)
        v[i] = s[i]
    end
    Array_72_Uint8(v...)
end

macro check_type(k, v, t)
    :(typeof($(esc(v))) <: $t || error($(esc(k))," must have type ",$t))
end

macro check_prop(k, prop, v, op, n)
    en = esc(n)
    :(($op)($prop($(esc(v))), $en) || error($(esc(k))," must have ",$prop," ",$op," ",$en))
end

function wcsmodify(w::wcsprm; kvs...)
    @assert w.flag != -1
    w.flag = 0
    naxis = w.naxis
    for (k,v) in kvs

        # double[naxis]
        if k in (:cdelt,:crder,:crota,:crpix,:crval,:csyer)
            @check_type k v Vector{Float64}
            @check_prop k length v (==) naxis
            unsafe_store_vec!(getfield(w,k), v)

        # char[72,naxis]
        elseif k in (:cname,:ctype,:cunit)
            @check_type k v Vector{ASCIIString}
            @check_prop k length v (==) naxis
            p = convert(Ptr{Array_72_Uint8}, getfield(w,k))
            x = Array_72_Uint8[convert(Array_72_Uint8,s) for s in v]
            unsafe_store_vec!(p, x)

        # pvcard[]
        elseif k === :pv
            npv = length(v)
            @check_type k v Vector{pvcard}
            @check_prop k length v (<=) w.npvmax
            unsafe_store_vec!(w.pv, v)
            w.npv = npv

        # pscard[]
        elseif k === :ps
            nps = length(v)
            @check_type k v Vector{pscard}
            @check_prop k length v (<=) w.npsmax
            unsafe_store_vec!(w.ps, v)
            w.nps = nps

        # double[naxis,naxis]
        elseif k in (:cd,:pc)
            @check_type k v Matrix{Float64}
            @check_prop k size v (==) (naxis,naxis)
            unsafe_store_vec!(getfield(w,k), vec(v'))

        # double
        elseif k in (:equinox,:latpole,:lonpole,:mjdavg,:mjdobs,
                     :restfrq,:restwav,:velangl,:velosys,:zsource)
            setfield!(w, k, convert(Float64,v))

        # int
        elseif k in (:colnum,)
            setfield!(w, k, convert(Cint,v))

        # char[72]
        elseif k in (:dateavg,:dateobs,:radesys,:specsys,:ssysobs,:ssyssrc,
                     :wcsname)
            setfield!(w, k, convert(Array_72_Uint8,v))

        # double[3]
        elseif k in (:obsgeo,)
            @check_type k v Vector{Float64}
            @check_prop k length v (==) 3
            w.obsgeo = Array_3_Cdouble(v[1], v[2], v[3])

        # char[4], but only uses first
        elseif k in (:alt,)
            @check_type k v Char
            (('A' <= v <= 'Z') || v == ' ') || error("alt must be 'A'-'Z' or ' '")
            x = 0x0
            w.alt = Array_4_Uint8(@compat(UInt8(v)), x, x, x)

        else
            error("unrecognized keyword argument \"$k\"")
        end
    end
end

function wcsprm(naxis::Integer; kvs...)
    w = wcsprm()
    stat = wcsini(1, naxis, w)
    @assert stat == 0
    finalizer(w, wcsfree)
    wcsmodify(w; kvs...)
    w
end

function wcscopy(alloc::Integer, src::wcsprm, dst::wcsprm)
    p = convert(Ptr{Cint}, C_NULL)
    wcssub(alloc, src, p, p, dst)
end

macro same_size(a, b)
    :(size($(esc(a))) == size($(esc(b))) || error($(string(a))," must have the same dimensions as ",$(string(b))))
end

function wcsp2s(wcs::wcsprm, pixcrd::Matrix{Float64};
                imgcrd::Matrix{Float64}=similar(pixcrd),
                phi::Matrix{Float64}=similar(pixcrd),
                theta::Matrix{Float64}=similar(pixcrd),
                world::Matrix{Float64}=similar(pixcrd),
                stat::Matrix{Cint}=similar(pixcrd,Cint))
    @same_size imgcrd pixcrd
    @same_size phi pixcrd
    @same_size theta pixcrd
    @same_size world pixcrd
    @same_size stat pixcrd
    (nelem, ncoord) = size(pixcrd)
    nelem >= wcs.naxis || error("number of rows must be greater than or equal to naxis")
    wcsp2s(wcs, ncoord, nelem, pointer(pixcrd), pointer(imgcrd),
           pointer(phi), pointer(theta), pointer(world), pointer(stat))
    world
end

function wcss2p(wcs::wcsprm, world::Matrix{Float64};
                phi::Matrix{Float64}=similar(world),
                theta::Matrix{Float64}=similar(world),
                imgcrd::Matrix{Float64}=similar(world),
                pixcrd::Matrix{Float64}=similar(world),
                stat::Matrix{Cint}=similar(world,Cint))
    @same_size phi world
    @same_size theta world
    @same_size imgcrd world
    @same_size pixcrd world
    @same_size stat world
    (nelem, ncoord) = size(world)
    nelem >= wcs.naxis || error("number of rows must be greater than or equal to naxis")
    wcss2p(wcs, ncoord, nelem, pointer(world), pointer(phi), pointer(theta),
           pointer(imgcrd), pointer(pixcrd), pointer(stat))
    pixcrd
end

function wcspih(header::ASCIIString; nkeyrec::Integer=div(length(header),80),
                                     relax::Integer=0, ctrl::Integer=0)
    @assert ctrl >= 0 # < 0 modifies the header
    nreject = Cint[0]
    nwcs = Cint[0]
    wcs = Ptr{wcsprm}[0]
    stat = wcspih(pointer(header), nkeyrec, relax, ctrl,
                  pointer(nreject), pointer(nwcs), pointer(wcs))
    @assert stat == 0
    p = wcs[1]
    a = wcsprm[unsafe_load(p,i) for i = 1:nwcs[1]]
    Libc.free(p)
    (a, @compat(Int(nreject[1])))
end

function wcsbth(header::ASCIIString; nkeyrec::Integer=div(length(header),80),
                                     relax::Integer=0, ctrl::Integer=0,
                                     keysel::Integer=0)
    @assert ctrl >= 0 # < 0 modifies the header
    nreject = Cint[0]
    nwcs = Cint[0]
    wcs = Ptr{wcsprm}[0]
    stat = wcsbth(pointer(header), nkeyrec, relax, ctrl, keysel, C_NULL,
                  pointer(nreject), pointer(nwcs), pointer(wcs))
    @assert stat == 0
    p = wcs[1]
    a = wcsprm[unsafe_load(p,i) for i = 1:nwcs[1]]
    Libc.free(p)
    (a, @compat(Int(nreject[1])))
end

function wcshdo(w::wcsprm; relax::Integer=0)
    nkeyrec = Cint[0]
    header = Ptr{UInt8}[0]
    stat = wcshdo(relax, w, pointer(nkeyrec), pointer(header))
    @assert stat == 0
    p = header[1]
    s = bytestring(p)
    Libc.free(p)
    s
end

end # module


text: abstract type PathIntegral end


text:
function allocate(datakind::Symbol, dtype::DataType, dims...)
    if datakind == :local
        return zeros(dtype, dims...)
    elseif datakind == :shared
        return SharedArray{dtype}(dims...)
    else
        ErrorException("Unknown kind: " * string(datakind)) |> throw
    end
end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
# sample from the Muller potential using BBK
using Printf
using Random
using LinearAlgebra
using ForwardDiff
using BasicMD
using Statistics

include("potentials.jl")

β = 0.1;
γ = 0.1;
M = [1.5, 0.75];
q₀ = [-0.75, 1.5];
p₀ = [0.0, 0.0];
Δt = 1e-2;
n_iters = 10^6;

V = x->Muller(x);
cfg = ForwardDiff.GradientConfig(V, q₀);
gradV! = (gradV, x)-> ForwardDiff.gradient!(gradV, V, x, cfg);

# second moment and energy obsrevables
f₁ = x-> x[1][1]^2;
f₂ = x-> V(x[1]);
observables = (f₁, f₂);

sampler = BBK(gradV!, β, γ, M, Δt)

Random.seed!(100);
Xvals = sample_trajectory([q₀, p₀], sampler, options=MDOptions(n_iters=n_iters,n_save_iters=n_save_iters));
f₁_estimate = mean(f₁.(Xvals));
f₂_estimate = mean(f₂.(Xvals));
@printf("f₁ estimate with %d samples: %g\n",length(Xvals), f₁_estimate);
@printf("f₂ estimate with %d samples: %g\n", length(Xvals), f₂_estimate);
Random.seed!(100);
observable_samples = sample_observables([q₀, p₀], sampler, observables, options=MDOptions(n_iters=n_iters,n_save_iters=n_save_iters));
@printf("f₁ estimate with %d samples: %g\n",length(observable_samples), mean(observable_samples[1,:]));
@printf("f₂ estimate with %d samples: %g\n", length(observable_samples), mean(observable_samples[2,:]));

max_stars_repo_name: gideonsimpson/JuBasicMD


text:
module NNNCVXF

using LinearAlgebra
using Random
using Statistics
using Printf
using PyPlot

using  Flux
import Flux.Optimise.update!
import Flux: gradient
import Flux.Losses.crossentropy
using  SetIntersectionProjection
using  InvertibleNetworks
using  Distributions
import Images
using CUDA
using Parameters

include("DataSetUtils.jl")
include("LossFunctions.jl")
include("NetworkOutputSampling.jl")
include("Train.jl")
include("PlotUtils.jl")
include("AugmentDataLabel.jl")

export TrOpts, Log

@with_kw  mutable struct TrOpts
           eval_every      ::Int  = 5
           alpha           ::Union{AbstractFloat,Vector{AbstractFloat},Array{Any}} = []
           batchsize       ::Int  = 1
           use_gpu         ::Bool = true
           lossf           ::Any  = []
           lossg           ::Any  = []
           active_channels ::Array{Int}  = [1]
           flip_dims       ::Array{Int}  = []
           permute_dims    ::Array{Int}  = []
           maxiter         ::Int         = 10
           opt             ::Any         = Flux.ADAM(1f-3)
           rand_grad_perc_zero  ::Float32     = 0f0
end

@with_kw mutable struct Log
  train     = []
  val       = []
  dc2_train = []
  dc2_val   = []
  IoU_train = []
  IoU_val   = []
end

end #end module