text:
#----------------------------------------
# OpenGL    GLSL      GLSL
# version   version   string
#----------------------------------------
#  2.0       110       "#version 110"
#  2.1       120
#  3.0       130
#  3.1       140
#  3.2       150       "#version 150"
#  3.3       330
#  4.0       400
#  4.1       410       "#version 410 core"
#  4.2       420
#  4.3       430
#  ES 2.0    100       "#version 100"
#  ES 3.0    300       "#version 300 es"
#----------------------------------------

const IMGUI_BACKEND_RENDERER_NAME = "imgui_impl_opengl3"
const GLFW_BACKEND_PLATFORM_NAME = "julia_imgui_glfw"

"""
    init(ctx::Context)
Initialize GLFW backend. Do initialize a IMGUI context before calling this function.
"""
function init(ctx::Context)
    ctx.Time = 0.0

    # setup back-end capabilities flags
    io::Ptr{ImGuiIO} = igGetIO()
    io.BackendFlags = unsafe_load(io.BackendFlags) | ImGuiBackendFlags_HasMouseCursors
    io.BackendFlags = unsafe_load(io.BackendFlags) | ImGuiBackendFlags_HasSetMousePos
    io.BackendFlags = unsafe_load(io.BackendFlags) | ImGuiBackendFlags_PlatformHasViewports
    # TODO: pending glfw 3.4
    # io.BackendFlags = unsafe_load(io.BackendFlags) | ImGuiBackendFlags_HasMouseHoveredViewport
    io.BackendPlatformName = pointer(GLFW_BACKEND_PLATFORM_NAME)

    io.BackendRendererName = pointer(IMGUI_BACKEND_RENDERER_NAME)
    io.BackendFlags = unsafe_load(io.BackendFlags) | ImGuiBackendFlags_RendererHasVtxOffset  # version ≥ 320
    io.BackendFlags = unsafe_load(io.BackendFlags) | ImGuiBackendFlags_RendererHasViewports

    # store the contextual object reference to IMGUI
    io.UserData = pointer_from_objref(ctx)

    # keyboard mapping
    c_set!(io.KeyMap, ImGuiKey_Tab, GLFW.KEY_TAB)
    c_set!(io.KeyMap, ImGuiKey_LeftArrow, GLFW.KEY_LEFT)
    c_set!(io.KeyMap, ImGuiKey_RightArrow, GLFW.KEY_RIGHT)
    c_set!(io.KeyMap, ImGuiKey_UpArrow, GLFW.KEY_UP)
    c_set!(io.KeyMap, ImGuiKey_DownArrow, GLFW.KEY_DOWN)
    c_set!(io.KeyMap, ImGuiKey_PageUp, GLFW.KEY_PAGE_UP)
    c_set!(io.KeyMap, ImGuiKey_PageDown, GLFW.KEY_PAGE_DOWN)
    c_set!(io.KeyMap, ImGuiKey_Home, GLFW.KEY_HOME)
    c_set!(io.KeyMap, ImGuiKey_End, GLFW.KEY_END)
    c_set!(io.KeyMap, ImGuiKey_Insert, GLFW.KEY_INSERT)
    c_set!(io.KeyMap, ImGuiKey_Delete, GLFW.KEY_DELETE)
    c_set!(io.KeyMap, ImGuiKey_Backspace, GLFW.KEY_BACKSPACE)
    c_set!(io.KeyMap, ImGuiKey_Space, GLFW.KEY_SPACE)
    c_set!(io.KeyMap, ImGuiKey_Enter, GLFW.KEY_ENTER)
    c_set!(io.KeyMap, ImGuiKey_Escape, GLFW.KEY_ESCAPE)
    c_set!(io.KeyMap, ImGuiKey_KeyPadEnter, GLFW.KEY_KP_ENTER)
    c_set!(io.KeyMap, ImGuiKey_A, GLFW.KEY_A)
    c_set!(io.KeyMap, ImGuiKey_C, GLFW.KEY_C)
    c_set!(io.KeyMap, ImGuiKey_V, GLFW.KEY_V)
    c_set!(io.KeyMap, ImGuiKey_X, GLFW.KEY_X)
    c_set!(io.KeyMap, ImGuiKey_Y, GLFW.KEY_Y)
    c_set!(io.KeyMap, ImGuiKey_Z, GLFW.KEY_Z)

    # set clipboard
    io.GetClipboardTextFn = GLFW_GET_CLIPBOARD_TEXT_FUNCPTR[]
    io.SetClipboardTextFn = GLFW_SET_CLIPBOARD_TEXT_FUNCPTR[]
    io.ClipboardUserData = Ptr{Cvoid}(ctx.Window.handle)

    # create mouse cursors
    ctx.MouseCursors[ImGuiMouseCursor_Arrow+1] = GLFW.CreateStandardCursor(GLFW.ARROW_CURSOR)
    ctx.MouseCursors[ImGuiMouseCursor_TextInput+1] = GLFW.CreateStandardCursor(GLFW.IBEAM_CURSOR)
    ctx.MouseCursors[ImGuiMouseCursor_ResizeNS+1] = GLFW.CreateStandardCursor(GLFW.VRESIZE_CURSOR)
    ctx.MouseCursors[ImGuiMouseCursor_ResizeEW+1] = GLFW.CreateStandardCursor(GLFW.HRESIZE_CURSOR)
    ctx.MouseCursors[ImGuiMouseCursor_Hand+1] = GLFW.CreateStandardCursor(GLFW.HAND_CURSOR)

    # prepare for GLFW 3.4+
    # ctx.MouseCursors[ImGuiMouseCursor_ResizeAll+1] = glfwCreateStandardCursor(GLFW_RESIZE_ALL_CURSOR)
    # ctx.MouseCursors[ImGuiMouseCursor_ResizeNESW+1] = glfwCreateStandardCursor(GLFW_RESIZE_NESW_CURSOR)
    # ctx.MouseCursors[ImGuiMouseCursor_ResizeNWSE+1] = glfwCreateStandardCursor(GLFW_RESIZE_NWSE_CURSOR)
    # ctx.MouseCursors[ImGuiMouseCursor_NotAllowed+1] = glfwCreateStandardCursor(GLFW_NOT_ALLOWED_CURSOR)
    ctx.MouseCursors[ImGuiMouseCursor_ResizeAll+1] = GLFW.CreateStandardCursor(GLFW.ARROW_CURSOR)
    ctx.MouseCursors[ImGuiMouseCursor_ResizeNESW+1] = GLFW.CreateStandardCursor(GLFW.ARROW_CURSOR)
    ctx.MouseCursors[ImGuiMouseCursor_ResizeNWSE+1] = GLFW.CreateStandardCursor(GLFW.ARROW_CURSOR)
    ctx.MouseCursors[ImGuiMouseCursor_NotAllowed+1] = GLFW.CreateStandardCursor(GLFW.ARROW_CURSOR)

    # chain GLFW callbacks
    ctx.PrevUserCallbackMousebutton = C_NULL
    ctx.PrevUserCallbackScroll = C_NULL
    ctx.PrevUserCallbackKey = C_NULL
    ctx.PrevUserCallbackChar = C_NULL
    ctx.PrevUserCallbackMonitor = C_NULL
    if ctx.InstalledCallbacks
        ctx.PrevUserCallbackMousebutton = glfwSetMouseButtonCallback(ctx.Window, @cfunction(ImGui_ImplGlfw_MouseButtonCallback, Cvoid, (GLFW.Window, GLFW.MouseButton, GLFW.Action, Cint)))
        ctx.PrevUserCallbackScroll = glfwSetScrollCallback(ctx.Window, @cfunction(ImGui_ImplGlfw_ScrollCallback, Cvoid, (GLFW.Window, Cdouble, Cdouble)))
        ctx.PrevUserCallbackKey = glfwSetKeyCallback(ctx.Window, @cfunction(ImGui_ImplGlfw_KeyCallback, Cvoid, (GLFW.Window, GLFW.Key, Cint, GLFW.Action, Cint)))
        ctx.PrevUserCallbackChar = glfwSetCharCallback(ctx.Window, @cfunction(ImGui_ImplGlfw_CharCallback, Cvoid, (GLFW.Window, Cuint)))
        ctx.PrevUserCallbackMonitor = glfwSetMonitorCallback(@cfunction(ImGui_ImplGlfw_MonitorCallback, Cvoid, (Ptr{Cvoid}, Cint)))
    end

    # update monitors the first time
    ImGui_ImplGlfw_UpdateMonitors(ctx)
    glfwSetMonitorCallback(@cfunction(ImGui_ImplGlfw_MonitorCallback, Cvoid, (Ptr{Cvoid}, Cint)))

    # the mouse update function expect PlatformHandle to be filled for the main viewport
    main_viewport::Ptr{ImGuiViewport} = igGetMainViewport()
    main_viewport.PlatformHandle = Ptr{Cvoid}(ctx.Window.handle)
    if Sys.iswindows()
        main_viewport.PlatformHandleRaw = ccall((:glfwGetWin32Window, GLFW.libglfw), Ptr{Cvoid}, (GLFW.Window,), ctx.Window)
    end

    if unsafe_load(io.ConfigFlags) & ImGuiConfigFlags_ViewportsEnable == ImGuiConfigFlags_ViewportsEnable
        ImGui_ImplGlfw_InitPlatformInterface(ctx)
        ImGui_ImplOpenGL3_InitPlatformInterface()
    end

    return true
end

"""
    shutdown((ctx::Context)
Clean up resources.
"""
function shutdown(ctx::Context)
    ImGui_ImplOpenGL3_ShutdownPlatformInterface()
    ImGui_ImplOpenGL3_DestroyDeviceObjects(ctx)

    ImGui_ImplGlfw_ShutdownPlatformInterface(ctx)

    if ctx.InstalledCallbacks
        ccall((:glfwSetMouseButtonCallback, libglfw), Ptr{Cvoid}, (GLFW.Window, Ptr{Cvoid}), ctx.Window, ctx.PrevUserCallbackMousebutton)
        ccall((:glfwSetScrollCallback, libglfw), Ptr{Cvoid}, (GLFW.Window, Ptr{Cvoid}), ctx.Window, ctx.PrevUserCallbackScroll)
        ccall((:glfwSetKeyCallback, libglfw), Ptr{Cvoid}, (GLFW.Window, Ptr{Cvoid}), ctx.Window, ctx.PrevUserCallbackKey)
        ccall((:glfwSetCharCallback, libglfw), Ptr{Cvoid}, (GLFW.Window, Ptr{Cvoid}), ctx.Window, ctx.PrevUserCallbackChar)
        ctx.InstalledCallbacks = false
    end

    for i = 1:length(ctx.MouseCursors)
        GLFW.DestroyCursor(ctx.MouseCursors[i])
        ctx.MouseCursors[i] = C_NULL
    end

    return true
end

function ImGui_ImplGlfw_UpdateMousePosAndButtons(ctx::Context)
    # update buttons
    io::Ptr{ImGuiIO} = igGetIO()
    for n = 0:length(ctx.MouseJustPressed)-1
        # if a mouse press event came, always pass it as "mouse held this frame",
        # so we don't miss click-release events that are shorter than 1 frame.
        is_down = ctx.MouseJustPressed[n+1] || GLFW.GetMouseButton(ctx.Window, n)
        c_set!(io.MouseDown, n, is_down)
        ctx.MouseJustPressed[n+1] = false
    end

    # update mouse position
    mouse_pos_backup = unsafe_load(io.MousePos)
    FLT_MAX = igGET_FLT_MAX()
    io.MousePos = ImVec2(-FLT_MAX, -FLT_MAX)
    io.MouseHoveredViewport = 0
    platform_io::Ptr{ImGuiPlatformIO} = igGetPlatformIO()
    vp = unsafe_load(platform_io.Viewports)
    viewport_ptrs = unsafe_wrap(Vector{Ptr{ImGuiViewport}}, vp.Data, vp.Size)
    for viewport in viewport_ptrs
        window = GLFW.Window(unsafe_load(viewport.PlatformHandle))
        @assert window != C_NULL

        if GLFW.GetWindowAttrib(window, GLFW.FOCUSED) != 0
            if unsafe_load(io.WantSetMousePos)
                x = mouse_pos_backup.x - unsafe_load(viewport.Pos.x)
                y = mouse_pos_backup.y - unsafe_load(viewport.Pos.y)
                GLFW.SetCursorPos(window, Cdouble(x), Cdouble(y))
            else
                mouse_x, mouse_y = GLFW.GetCursorPos(GLFW.Window(window))
                if unsafe_load(io.ConfigFlags) & ImGuiConfigFlags_ViewportsEnable == ImGuiConfigFlags_ViewportsEnable
                    # Multi-viewport mode: mouse position in OS absolute coordinates (io.MousePos is (0,0) when the mouse is on the upper-left of the primary monitor)
                    window_x, window_y = GLFW.GetWindowPos(GLFW.Window(window))
                    io.MousePos = ImVec2(Cfloat(mouse_x + window_x), Cfloat(mouse_y + window_y))
                else
                    # Single viewport mode: mouse position in client window coordinates (io.MousePos is (0,0) when the mouse is on the upper-left corner of the app window)
                    io.MousePos = ImVec2(Cfloat(mouse_x), Cfloat(mouse_y))
                end
            end

        end
        for n = 0:length(ctx.MouseJustPressed)-1
            c_set!(io.MouseDown, n, c_get(io.MouseDown, n) || GLFW.GetMouseButton(window, n))
        end
    end
    # TODO: pending glfw 3.4 GLFW_HAS_MOUSE_PASSTHROUGH
    # TODO: _WIN32
end

function ImGui_ImplGlfw_UpdateMouseCursor(ctx::Context)
    io::Ptr{ImGuiIO} = igGetIO()
    if (unsafe_load(io.ConfigFlags) & ImGuiConfigFlags_NoMouseCursorChange == ImGuiConfigFlags_NoMouseCursorChange) ||
        GLFW.GetInputMode(ctx.Window, GLFW.CURSOR) == GLFW.CURSOR_DISABLED
        return nothing
    end

    imgui_cursor = igGetMouseCursor()
    platform_io::Ptr{ImGuiPlatformIO} = igGetPlatformIO()
    vp = unsafe_load(platform_io.Viewports)
    viewport_ptrs = unsafe_wrap(Vector{Ptr{ImGuiViewport}}, vp.Data, vp.Size)
    for viewport in viewport_ptrs
        window = unsafe_load(viewport.PlatformHandle)
        @assert window != C_NULL
        if imgui_cursor == ImGuiMouseCursor_None || unsafe_load(io.MouseDrawCursor)
            # hide OS mouse cursor if imgui is drawing it or if it wants no cursor
            GLFW.SetInputMode(window, GLFW.CURSOR, GLFW.CURSOR_HIDDEN)
        else
            # show OS mouse cursor
            cursor = ctx.MouseCursors[imgui_cursor+1]
            GLFW.SetCursor(window, cursor != C_NULL ? ctx.MouseCursors[imgui_cursor+1] : ctx.MouseCursors[ImGuiMouseCursor_Arrow+1])
            GLFW.SetInputMode(window, GLFW.CURSOR, GLFW.CURSOR_NORMAL)
        end
    end

    return nothing
end

function ImGui_ImplGlfw_UpdateMonitors(ctx::Context)
    platform_io::Ptr{ImGuiPlatformIO} = igGetPlatformIO()
    monitors_count = Cint(0)
	ptr = @c glfwGetMonitors(&monitors_count)
	glfw_monitors = unsafe_wrap(Array, ptr, monitors_count)
    monitors_ptr::Ptr{ImGuiPlatformMonitor} = Libc.malloc(monitors_count * sizeof(ImGuiPlatformMonitor))
    for i = 1:monitors_count
        glfw_monitor = glfw_monitors[i]
        mptr::Ptr{ImGuiPlatformMonitor} = monitors_ptr + (i-1) * sizeof(ImGuiPlatformMonitor)

        x, y = Cint(0), Cint(0)
        @c glfwGetMonitorPos(glfw_monitor, &x, &y)
        vid_mode = unsafe_load(glfwGetVideoMode(glfw_monitor))
        mptr.MainPos = ImVec2(x, y)
        mptr.MainSize = ImVec2(vid_mode.width, vid_mode.height)
        mptr.WorkPos = ImVec2(x, y)
        mptr.WorkSize = ImVec2(vid_mode.width, vid_mode.height)

        w, h = Cint(0), Cint(0)
        @c glfwGetMonitorWorkarea(glfw_monitors[i], &x, &y, &w, &h)
        if w > 0 && h > 0
            mptr.WorkPos = ImVec2(Cfloat(x), Cfloat(y))
            mptr.WorkSize = ImVec2(Cfloat(w), Cfloat(h))
        end

        x_scale, y_scale = Cfloat(0), Cfloat(0)
        @c glfwGetMonitorContentScale(glfw_monitor, &x_scale, &y_scale)
        mptr.DpiScale = x_scale
    end

    platform_io.Monitors = ImVector_ImGuiPlatformMonitor(monitors_count, monitors_count, monitors_ptr)
    ctx.WantUpdateMonitors = false

    return nothing
end

function new_frame(ctx::Context)
    if ctx.ShaderHandle != C_NULL
        ImGui_ImplOpenGL3_CreateDeviceObjects(ctx)
    end

    io::Ptr{ImGuiIO} = igGetIO()
    @assert ImFontAtlas_IsBuilt(unsafe_load(io.Fonts)) "Font atlas not built! It is generally built by the renderer back-end. Missing call to renderer _NewFrame() function? e.g. ImGui_ImplOpenGL3_NewFrame()."

    # setup display size (every frame to accommodate for window resizing)
    w, h = Cint(0), Cint(0)
    display_w, display_h = Cint(0), Cint(0)
    @c glfwGetWindowSize(ctx.Window, &w, &h)
    @c glfwGetFramebufferSize(ctx.Window, &display_w, &display_h)
    io.DisplaySize = ImVec2(Cfloat(w), Cfloat(h))
    if w > 0 && h > 0
        w_scale = Cfloat(display_w / w)
        h_scale = Cfloat(display_h / h)
        io.DisplayFramebufferScale = ImVec2(w_scale, h_scale)
    end

    ctx.WantUpdateMonitors && ImGui_ImplGlfw_UpdateMonitors(ctx)

    # setup time step
    current_time = glfwGetTime()
    io.DeltaTime = ctx.Time > 0.0 ? Cfloat(current_time - ctx.Time) : Cfloat(1.0/60.0)
    ctx.Time = current_time

    ImGui_ImplGlfw_UpdateMousePosAndButtons(ctx)
    ImGui_ImplGlfw_UpdateMouseCursor(ctx)

    # TODO: Gamepad navigation mapping
    # ImGui_ImplGlfw_UpdateGamepads()
end

render(ctx::Context) = ImGui_ImplOpenGL3_RenderDrawData(ctx, igGetDrawData())


text:
"""Calculates line of sight distance between two points on the earth (earth radius and conversion to miles hardcoded below)"""
function haversine(x, y)
    @inbounds begin
        # longitudes
        Δλ = deg2rad(y[1] - x[1])

        # latitudes
        φ₁ = deg2rad(x[2])
        φ₂ = deg2rad(y[2])
    end

    Δφ = φ₂ - φ₁

    # haversine formula
    a = sin(Δφ/2)^2 + cos(φ₁)*cos(φ₂)*sin(Δλ/2)^2

    # distance on the sphere
    distance_in_meters = round(2 * 6.3781e6 * asin( min(√a, one(a)) )) # take care of floating point errors
    distance_in_miles = round(distance_in_meters / 1_609, 2)
    distance_in_miles
end

ext: jl


text:
#######################################################################
function marginaline(posterior::FullNormal, σ2)
#######################################################################

    # dimension of problem is D

    D = round(Int, length(posterior)/2)

    # marginals, conditionals

    μ, Σ = mean(posterior), cov(posterior)

    a_idx, b_idx = collect(1:D), collect(D+1:2D)

    μ_a = μ[a_idx]

    μ_b = μ[b_idx]

    Σ_aa = Σ[a_idx, a_idx]

    Σ_ab = Σ[a_idx, b_idx]

    Σ_ba = Σ[b_idx, a_idx]

    Σ_bb = Σ[b_idx, b_idx]

    Σ_a_given_b = Σ_aa - Σ_ab*(Σ_bb\Σ_ba)

    σ2Diag = Diagonal(σ2)

    function linemeansigma(t)

        local mu = μ_a * t + μ_b

        local K = t*(Σ_ab/Σ_bb) + I

        local Sigma = σ2Diag + t * Σ_a_given_b * t  +  K*Σ_bb*K'

        Sigma = (Sigma + Sigma')/2

        mu, Sigma

    end

end

#######################################################################
function marginaline(posterior::DiagNormal, σ2)
#######################################################################

    # dimension of problem is D

    D = round(Int, length(posterior)/2)

    # marginals, conditionals

    μ, Σ = mean(posterior), cov(posterior)

    a_idx, b_idx = collect(1:D), collect(D+1:2D)

    μ_a = μ[a_idx]

    μ_b = μ[b_idx]

    Σ_aa = Diagonal(Σ[a_idx, a_idx])

    # Σ_ab = Σ[a_idx, b_idx]      # These are zero because of diagonal
    #                             # covariance matrix
    # Σ_ba = Σ[b_idx, a_idx]

    Σ_bb = Diagonal(Σ[b_idx, b_idx])

    Σ_a_given_b = Σ_aa # - Σ_ab*(Σ_bb\Σ_ba) # The last term is zero, same reason as above

    σ2Diag = Diagonal(σ2)

    function linemeansigma(t)

        local mu = μ_a * t + μ_b

        # local K = t*(Σ_ab/Σ_bb) + I # The first term is zero, hence K is idenity

        local Sigma = σ2Diag + t * Σ_a_given_b * t  +  Σ_bb# K*Σ_bb*K'

        Sigma = (Sigma + Sigma')/2

        mu, Sigma

    end

end


text:
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the compressible Euler equations
gamma = 1.4
equations = CompressibleEulerEquations2D(gamma)

# Same as the default Trixi example, but v2 is a non-symmetric perturbation
function initial_condition_khi_nonsymmetric_perturbation(x, t, equations::CompressibleEulerEquations2D)
    # change discontinuity to tanh
    # typical resolution 128^2, 256^2
    # domain size is [-1,+1]^2
    slope = 15
     # discontinuous function with value 2 for -.5 <= x <= .5 and 0 elsewhere
    B = tanh(slope * x[2] + 7.5) - tanh(slope * x[2] - 7.5)
    rho = 0.5 + 0.75 * B
    v1 = 0.5 * (B - 1)
    # scale by (1 + .01 * sin(pi*x) * sin(pi*y)) to break symmetry
    v2 = 0.1 * sin(2 * pi * x[1]) * (1.0 + 1e-2 * sin(pi * x[1]) * sin(pi * x[2]))
    p = 1.0
    return prim2cons(SVector(rho, v1, v2, p), equations)
end
initial_condition = initial_condition_khi_nonsymmetric_perturbation

# surface_flux = FluxPlusDissipation(flux_ranocha, DissipationLocalLaxFriedrichs())
surface_flux = flux_lax_friedrichs
volume_flux  = flux_ranocha
polydeg = 3
basis = LobattoLegendreBasis(polydeg)
indicator_sc = IndicatorHennemannGassner(equations, basis,
                                         alpha_max=0.0025,
                                         alpha_min=0.0001,
                                         alpha_smooth=true,
                                         variable=density_pressure)
volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
                                                 volume_flux_dg=volume_flux,
                                                 volume_flux_fv=surface_flux)
# volume_integral = VolumeIntegralFluxDifferencing(flux_ranocha)
solver = DGSEM(basis, surface_flux, volume_integral)

coordinates_min = (-1.0, -1.0)
coordinates_max = ( 1.0,  1.0)
cells_per_dimension = 32
mesh = TreeMesh(coordinates_min, coordinates_max,
                initial_refinement_level=Int(log2(cells_per_dimension)),
                n_cells_max=100_000)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)

###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 3.0)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 500
analysis_callback = AnalysisCallback(semi, interval=analysis_interval)

alive_callback = AliveCallback(alive_interval=100, analysis_interval=analysis_interval)

stepsize_callback = StepsizeCallback(cfl=0.9)

callbacks = CallbackSet(summary_callback,
                        analysis_callback,
                        alive_callback,
                        stepsize_callback)

###############################################################################
# run the simulation

tsave = LinRange(tspan..., 2)
# use_positivity_limiting = true
# if use_positivity_limiting==true
    limiter! = PositivityPreservingLimiterZhangShu(thresholds=(5.0e-6, 5.0e-6),
                                                variables=(Trixi.density, pressure))
    sol = solve(ode, SSPRK43(limiter!), saveat=tsave, save_everystep=false, callback=callbacks);
# else
#     sol = solve(ode, RDPK3SpFSAL49(), abstol = 1.0e-7, reltol = 1.0e-7,
#                 save_everystep = false, saveat = tsave, callback = callbacks)
# end

summary_callback() # print the timer summary


text:
const QueryKeys = AbstractSelection
const TermKeys = Union{Symbol,OtherSelection,ElseSelection}
const Renaming = AbstractRenaming
const Trans = AbstractTransformation
const MaybeRenaming = Union{Nothing,AbstractRenaming}
const MaybeTrans = Union{Nothing,AbstractTransformation}
const MaybeComp = Union{Nothing,Composition}

struct SelectionQuery{S<:QueryKeys, R<:MaybeRenaming, T<:MaybeTrans} <: AbstractMetaSelection
    s::S
    r::R
    t::T
end
SelectionQuery(s::SelectionQuery) = s
SelectionQuery(s::S) where {S<:QueryKeys} = SelectionQuery(s, nothing, nothing)

SelectionQuery(s::S, r::R) where {S<:QueryKeys, R<:MaybeRenaming} =
    SelectionQuery(x, r, nothing)
SelectionQuery(p::Pair{S,R}) where {S<:QueryKeys, R<:MaybeRenaming} =
    SelectionQuery(first(p), last(p), nothing)

SelectionQuery(s::S, t::T) where {S<:QueryKeys, T<:MaybeTrans} =
    SelectionQuery(x, nothing, t)
SelectionQuery(p::Pair{S,T}) where {S<:QueryKeys, T<:MaybeTrans} =
    SelectionQuery(first(p), nothing, last(p))

SelectionQuery(s::S, t::T, r::R) where {S<:QueryKeys, R<:Renaming, T<:Trans} =
    SelectionQuery(s, r, t)
SelectionQuery(p::Pair{S,Pair{R,T}}) where {S<:QueryKeys, R<:MaybeRenaming, T<:MaybeTrans} =
    SelectionQuery(first(p), first(last(p)), last(last(p)))
SelectionQuery(p::Pair{S,Pair{T,R}}) where {S<:QueryKeys, R<:MaybeRenaming, T<:MaybeTrans} =
    SelectionQuery(first(p), last(last(p)), first(last(p)))

struct SelectionTerm{S<:TermKeys, R<:MaybeComp, T<:MaybeComp} <: AbstractMetaSelection
    s::S
    r::R
    t::T
end
colname(s::SelectionTerm) = s.s
keyfunc(s::SelectionTerm) = s.r
valfunc(s::SelectionTerm) = s.t
Base.iterate(s::SelectionTerm) = (colname(s), (keyfunc(s), valfunc(s)))
Base.iterate(s::SelectionTerm, state) = (first(state), Base.tail(state))
Base.iterate(s::SelectionTerm, ::Type{Tuple}) = nothing

struct SelectionResult{S} <: AbstractVector{S}
    s::Vector{S}
    SelectionResult(s::AbstractVector{S}) where {S<:SelectionTerm} = new{S}(collect(s))
end
Base.size(s::SelectionResult) = size(s.s)
Base.getindex(s::SelectionResult, i) = getindex(s.s, i)
Base.iterate(x::SelectionResult) = iterate(x.s)
Base.iterate(x::SelectionResult, state) = iterate(x.s, state)
function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, s::SelectionResult)
    print(io, length(s.s), "-element SelectionResult\n")
    foreach(x->println(io, " ", x), s.s)
end

function SelectionResult(xs::Vector{S}, r::RenamingSymbols, t::MaybeTrans) where {
        S <: Union{Symbol,OtherSelection,ElseSelection}
    }
    if length(xs) == length(r.s)
        SelectionResult([
            SelectionTerm(x, Composition(RenamingSymbol(s)), Composition(t))
            for (x, s)
            in zip(xs, r.s)
        ])
    else
        @warn("Renaming array had different length ($(length(r.s))) than target selections ($(length(xs))), renaming skipped.")
        SelectionResult([SelectionTerm(x, nothing, Composition(t)) for x in xs])
    end
end

function SelectionResult(xs::Vector{S}, r::MaybeRenaming, t::MaybeTrans) where {
        S <: Union{Symbol,OtherSelection,ElseSelection}
    }
    SelectionResult([SelectionTerm(x, Composition(r), Composition(t)) for x in xs])
end
function SelectionResult(xs::Vector{S}, r::Union{Composition, MaybeRenaming}, t::Union{Composition, MaybeTrans}) where {
        S <: Union{Symbol,OtherSelection,ElseSelection}
    }
    SelectionResult([SelectionTerm(x, Composition(r), Composition(t)) for x in xs])
end

function SelectionResult(x::S, r::MaybeRenaming, t::MaybeTrans) where {
        S <: Union{Symbol,OtherSelection,ElseSelection}
    }
    SelectionResult([SelectionTerm(x, Composition(r), Composition(t))])
end

function SelectionResult(
        xs::AbstractVector{S},
        rs::AbstractVector{<:MaybeComp},
        ts::AbstractVector{<:MaybeComp}) where {
        S <: Union{Symbol,OtherSelection,ElseSelection}
    }
    SelectionResult([SelectionTerm(x, r, t) for (x, r, t) in zip(xs, rs, ts)])
end

function SelectionResult(xs::SelectionResult, rn::MaybeRenaming, tn::MaybeTrans)
    SelectionResult([SelectionTerm(x, extend(r, rn), extend(t, tn)) for (x, r, t) in xs])
end

colnames(x::SelectionResult) = map(colname, x)
keyfuncs(x::SelectionResult) = map(keyfunc, x)
valfuncs(x::SelectionResult) = map(valfunc, x)

Base.show(io::IO, s::SelectionQuery) = (print(io, "SelectionQuery("); _print_srt(io, s); print(io, ")"))
Base.show(io::IO, s::SelectionTerm) = (print(io, "SelectionTerm("); _print_srt(io, s); print(io, ")"))
function _print_srt(io::IO, s)
    print(io, s.s isa Symbol ? ":" : "", s.s)
    isnothing(s.r) || print(io, " => ", s.r)
    isnothing(s.t) || print(io, " => ", s.t)
end


text:
const BOLT = b"\x60\x60\xB0\x17"
const BOLT_VERSION = 1
# Versions are coded as UInt32, so 4*4 elements for 4 versions
const RAW_BOLT_VERSIONS = pack([0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0])[3:end]

const CLIENT = Dict(
    "INIT"          => 0x01,
    "ACK_FAILURE"   => 0x0e,
    "RESET"         => 0x0f,
    "RUN"           => 0x10,
    "DISCARD_ALL"   => 0x2f,
    "PULL_ALL"      => 0x3f,
)

const SERVER = Dict(
    "SUCESS"    => 0x70,
    "RECORD"    => 0x71,
    "IGNORED"   => 0x7e,
    "FAILURE"   => 0x7f,
)

const DEFAULT_USER_AGENT = "boltkit/0.0" # Unicode ?
const MAX_CHUNK_SIZE = 65535
const SOCKET_TIMEOUT = 5


text:
## Look at the NCAA data using our simple rating system.

## Read the data
using CSV
using DataFrames
content = CSV.read("ncaa-2018.csv")

## Convert teams into numbers and also build our game table
m = size(content,1)
teams_dict = Dict{String,Int}()
teams = Vector{String}()
index = 1
data = Array{Int64}(undef, m, 4)
for i = 1 : m
    global index
    team = content[:Team][i]
    opponent = content[:Opponent][i]
    teamscore = content[Symbol("Team Score")][i]
    oppscore = content[Symbol("Opponent Score")][i]
    if !haskey(teams_dict, team)
        push!(teams, team)
        teams_dict[team] = index
        index = index + 1
    end

    if !haskey(teams_dict, opponent)
        push!(teams, opponent)
        teams_dict[opponent] = index
        index = index + 1
    end

    data[i, 1] = teams_dict[team]
    data[i, 2] = teams_dict[opponent]
    data[i, 3] = teamscore
    data[i, 4] = oppscore
end

## Build the least-squares model
B = zeros(m, length(teams))
p = zeros(m)
for i=1:m
    # form the model
    ti = data[i,1]
    tj = data[i,2]
    p[i] = data[i,3] - data[i,4]
    B[i,ti] = 1
    B[i,tj] = -1
end
## Solve the model
# This solves as a least squares problem.
s = B\p

##
rank = sortperm(s,rev=true)
##
teams[rank]
##
## Numerical Computing in Data Science top 10 in 2018 are ..
using Printf
for i=1:10
    @printf("%20s - %2i - %.3f\n ",
        teams[rank[i]], i, s[rank[i]])
end


text:
baremodule EventTracker

export @recordinterval, @recordpoint

using EventTrackerBase:
    @recordinterval,
    @recordpoint,
    EVENT_INTERVAL,
    EVENT_POINT,
    Event,
    EventKind,
    Location,
    append_events!,
    append_locations!,
    clear

function summary end
function summary_dataframe end

function event_table end
function event_dataframe end
function location_table end
function location_dataframe end

function stacks end

module Implementations

using BangBang.Extras: modify!!
using DataFrames: DataFrame, combine, groupby, select!, Not, leftjoin
using RecipesBase: RecipesBase, @recipe, @series
using Statistics: mean
using StructArrays: StructVector
using Tables: Tables

using EventTrackerBase:
    EVENT_INTERVAL,
    EVENT_POINT,
    Event,
    EventID,
    Location,
    LocationID,
    TimeNS,
    append_events!,
    append_locations!

using ..EventTracker: EventTracker

include("utils.jl")
include("analysis.jl")
include("stacks.jl")
include("plots.jl")

end

Implementations.define_docstrings()

end  # module


text:
import StatsBase: sample, Weights

"""
$(TYPEDEF)

`DaviesGenerator` defines a Davies generator.

# Fields

$(FIELDS)
"""
struct DaviesGenerator <: AbstractLiouvillian
    "System bath coupling operators"
    coupling::AbstractCouplings
    "Spectrum density"
    γ::Any
    "Lambshift spectral density"
    S::Any
end

function (D::DaviesGenerator)(du, ρ, ω_ba, v, s::Real)
    γm = D.γ.(ω_ba)
    sm = D.S.(ω_ba)
    for op in D.coupling(s)
        A = v' * op * v
        davies_update!(du, ρ, A, γm, sm)
    end
end

function (D::DaviesGenerator)(du, ρ, ω_ba, s::Real)
    γm = D.γ.(ω_ba)
    sm = D.S.(ω_ba)
    for op in D.coupling(s)
        davies_update!(du, ρ, op, γm, sm)
    end
end

function davies_update!(du, u, A, γ, S)
    A2 = abs2.(A)
    γA = γ .* A2
    Γ = sum(γA, dims=1)
    dim = size(du, 1)
    for a = 1:dim
        for b = 1:a - 1
            @inbounds du[a, a] += γA[a, b] * u[b, b] - γA[b, a] * u[a, a]
            @inbounds du[a, b] +=
                -0.5 * (Γ[a] + Γ[b]) * u[a, b] +
                γ[1, 1] * A[a, a] * A[b, b] * u[a, b]
        end
        for b = a + 1:dim
            @inbounds du[a, a] += γA[a, b] * u[b, b] - γA[b, a] * u[a, a]
            @inbounds du[a, b] +=
                -0.5 * (Γ[a] + Γ[b]) * u[a, b] +
                γ[1, 1] * A[a, a] * A[b, b] * u[a, b]
        end
    end
    H_ls = Diagonal(sum(S .* A2, dims=1)[1, :])
    axpy!(-1.0im, H_ls * u - u * H_ls, du)
end

function update_cache!(cache, D::DaviesGenerator, ω_ba, v, s::Real)
    len = size(ω_ba, 1)
    γm = D.γ.(ω_ba)
    sm = D.S.(ω_ba)
    for op in D.coupling(s)
        A2 = abs2.(v' * op * v)
        for b = 1:len
            for a = 1:len
                @inbounds cache[b, b] -=
                    A2[a, b] * (0.5 * γm[a, b] + 1.0im * sm[a, b])
            end
        end
    end
end

"""
$(TYPEDEF)

Defines correlated Davies generator

# Fields

$(FIELDS)
"""
struct CorrelatedDaviesGenerator <: AbstractLiouvillian
    "System bath coupling operators"
    coupling::AbstractCouplings
    "Spectrum density"
    γ::Any
    "Lambshift spectral density"
    S::Any
    "Indices to iterate"
    inds::Any
end

function (D::CorrelatedDaviesGenerator)(du, ρ, ω_ba, s::Real)
    for (α, β) in D.inds
        γm = D.γ[α,β].(ω_ba)
        sm = D.S[α,β].(ω_ba)
        Aα = D.coupling[α](s)
        Aβ = D.coupling[β](s)
        correlated_davies_update!(du, ρ, Aα, Aβ, γm, sm)
    end
end

function correlated_davies_update!(du, u, Aα, Aβ, γ, S)
    A2 = transpose(Aα) .* Aβ
    γA = γ .* A2
    Γ = sum(γA, dims=1)
    dim = size(du, 1)
    for a = 1:dim
        for b = 1:a - 1
            @inbounds du[a, a] += γA[a, b] * u[b, b] - γA[b, a] * u[a, a]
            @inbounds du[a, b] +=
                -0.5 * (Γ[a] + Γ[b]) * u[a, b] +
                γ[1, 1] * Aα[a, a] * Aβ[b, b] * u[a, b]
        end
        for b = a + 1:dim
            @inbounds du[a, a] += γA[a, b] * u[b, b] - γA[b, a] * u[a, a]
            @inbounds du[a, b] +=
                -0.5 * (Γ[a] + Γ[b]) * u[a, b] +
                γ[1, 1] * Aα[a, a] * Aβ[b, b] * u[a, b]
        end
    end
    H_ls = Diagonal(sum(S .* A2, dims=1)[1, :])
    axpy!(-1.0im, H_ls * u - u * H_ls, du)
end

"""
$(TYPEDEF)

Defines the one-sided AME Liouvillian operator.

# Fields

$(FIELDS)
"""
struct OneSidedAMELiouvillian <: AbstractLiouvillian
    "System bath coupling operators"
    coupling::AbstractCouplings
    "Spectrum density"
    γ::Any
    "Lambshift spectral density"
    S::Any
    "Indices to iterate"
    inds::Any
end

function (A::OneSidedAMELiouvillian)(dρ, ρ, ω_ba, v, s::Real)
    for (α, β) in A.inds
        γm = A.γ[α,β].(ω_ba)
        sm = A.S[α,β].(ω_ba)
        Aα = v' * A.coupling[α](s) * v
        Λ = (0.5 * γm + 1.0im * sm) .* Aα
        Aβ =  v' * A.coupling[β](s) * v
        𝐊₂ = Aβ * Λ * ρ - Λ * ρ * Aβ
        𝐊₂ = 𝐊₂ + 𝐊₂'
        axpy!(-1.0, 𝐊₂, dρ)
    end
end

function (A::OneSidedAMELiouvillian)(du, u, ω_ba, s::Real)
    for (α, β) in A.inds
        γm = A.γ[α,β].(ω_ba)
        sm = A.S[α,β].(ω_ba)
        Aα = A.coupling[α](s)
        Aβ = A.coupling[β](s)
        Λ = (0.5 * γm + 1.0im * sm) .* Aα
        𝐊₂ = Aβ * Λ * u - Λ * u * Aβ
        𝐊₂ = 𝐊₂ + 𝐊₂'
        axpy!(-1.0, 𝐊₂, du)
    end
end


text:
# WHERE clause.

mutable struct WhereClause <: AbstractSQLClause
    over::Union{SQLClause, Nothing}
    condition::SQLClause

    WhereClause(;
                over = nothing,
                condition) =
        new(over, condition)
end

WhereClause(condition; over = nothing) =
    WhereClause(over = over, condition = condition)

"""
    WHERE(; over = nothing, condition)
    WHERE(condition; over = nothing)

A `WHERE` clause.

# Examples

```jldoctest
julia> c = FROM(:location) |>
           WHERE(OP("=", :zip, "60614")) |>
           SELECT(:location_id);

julia> print(render(c))
SELECT "location_id"
FROM "location"
WHERE ("zip" = '60614')
```
"""
WHERE(args...; kws...) =
    WhereClause(args...; kws...) |> SQLClause

dissect(scr::Symbol, ::typeof(WHERE), pats::Vector{Any}) =
    dissect(scr, WhereClause, pats)

function PrettyPrinting.quoteof(c::WhereClause, ctx::QuoteContext)
    ex = Expr(:call, nameof(WHERE), quoteof(c.condition, ctx))
    if c.over !== nothing
        ex = Expr(:call, :|>, quoteof(c.over, ctx), ex)
    end
    ex
end

rebase(c::WhereClause, c′) =
    WhereClause(over = rebase(c.over, c′), condition = c.condition)


text:
##############################################################################
##
## group combines multiple refs
## Missings have a ref of 0
##
##############################################################################

mutable struct GroupedArray{N} <: AbstractArray{UInt32, N}
    refs::Array{UInt32, N}   # refs must be between 0 and n. 0 means missing
    n::Int                   # Number of potential values (= maximum(refs))
end

size(g::GroupedArray) = size(g.refs)
@propagate_inbounds getindex(g::GroupedArray, i::Int) = getindex(g.refs, i)
@propagate_inbounds getindex(g::GroupedArray, I) = getindex(g.refs, I)

group(xs::GroupedArray) = xs

function group(xs::AbstractArray)
	refs = Array{UInt32}(undef, size(xs))
	invpool = Dict{eltype(xs), UInt32}()
	n = 0
	z = UInt32(0)
	@inbounds for i in eachindex(xs)
		x = xs[i]
		if x === missing
			refs[i] = 0
		else
			lbl = get(invpool, x, z)
			if !iszero(lbl)
				refs[i] = lbl
			else
				n += 1
				refs[i] = n
				invpool[x] = n
			end
		end
	end
	return GroupedArray{ndims(xs)}(refs, n)
end

function group(args...)
    g1 = deepcopy(group(args[1]))
    for j in 2:length(args)
        gj = group(args[j])
        size(g1) == size(gj) || throw(DimensionMismatch(
            "cannot match array of size $(size(g1)) with array of size $(size(gj))"))
        combine!(g1, gj)
    end
    factorize!(g1)
end

function combine!(g1::GroupedArray, g2::GroupedArray)
    @inbounds for i in eachindex(g1.refs, g2.refs)
        # if previous one is missing or this one is missing, set to missing
        g1.refs[i] = (g1.refs[i] == 0 || g2.refs[i] == 0) ? 0 : g1.refs[i] + (g2.refs[i] - 1) * g1.n
    end
    g1.n = g1.n * g2.n
    return g1
end

# An in-place version of group() that relabels the refs
function factorize!(g::GroupedArray{N}) where {N}
    refs = g.refs
    invpool = Dict{UInt32, UInt32}()
    n = 0
	z = UInt32(0)
    @inbounds for i in eachindex(refs)
        x = refs[i]
        if !iszero(x)
            lbl = get(invpool, x, z)
            if !iszero(lbl)
                refs[i] = lbl
            else
                n += 1
                refs[i] = n
                invpool[x] = n
            end
        end
    end
    return GroupedArray{N}(refs, n)
end

##############################################################################
##
## Cluster Covariance
##
##############################################################################

struct ClusterCovariance <: CovarianceEstimator
    clusternames::Tuple{Symbol, Vararg{Symbol}} # names are not allowed to be empty
    clusters::Union{Tuple{Vararg{GroupedArray}}, Nothing}
end

"""
    cluster(names::Symbol...)

Estimate variance-covariance matrix with a cluster-robust estimator.

# Arguments
- `names::Symbol...`: column names of variables that represent the clusters.
"""
cluster(ns::Symbol...) = ClusterCovariance((ns...,), nothing)
cluster(ns::NTuple{N, Symbol}, gs::NTuple{N, GroupedArray}) where N =
    ClusterCovariance(ns, gs)

"""
    names(vce::ClusterCovariance)

Return column names of variables used to form clusters for `vce`.
"""
names(v::ClusterCovariance) = v.clusternames

length(v::ClusterCovariance) = length(names(v))

show(io::IO, ::ClusterCovariance) = print(io, "Cluster-robust covariance estimator")

function show(io::IO, ::MIME"text/plain", v::ClusterCovariance)
    print(io, length(v), "-way cluster-robust covariance estimator:")
    for n in names(v)
        print(io, "\n  ", n)
    end
end

function Vcov.completecases(table, v::ClusterCovariance)
    Tables.istable(table) || throw(ArgumentError("completecases requires a table input"))
    out = trues(length(Tables.rows(table)))
    columns = Tables.columns(table)
    for name in names(v)
        out .&= .!ismissing.(Tables.getcolumn(columns, name))
    end
    return out
end

function materialize(table, v::ClusterCovariance)
    Tables.istable(table) || throw(ArgumentError("materialize requires a table input"))
    ns = names(v)
    return cluster(ns, map(n->group(Tables.getcolumn(table, n)), ns))
end

"""
    nclusters(vce::ClusterCovariance)

Return the number of clusters for each dimension/way of clustering.
"""
function nclusters(v::ClusterCovariance)
    NamedTuple{names(v)}(map(x -> x.n, v.clusters))
end

function df_FStat(x::RegressionModel, v::ClusterCovariance, ::Bool)
    minimum(nclusters(v)) - 1
end

function S_hat(x::RegressionModel, v::ClusterCovariance)
    # Cameron, Gelbach, & Miller (2011): section 2.3
    dim = size(modelmatrix(x), 2) * size(residuals(x), 2)
    S = zeros(dim, dim)
    for c in combinations(1:length(v))
        g = group((v.clusters[i] for i in c)...)
        S += (-1)^(length(c) - 1) * helper_cluster(modelmatrix(x), residuals(x), g)
    end
    # scale total vcov estimate by ((N-1)/(N-K)) * (G/(G-1))
    # another option would be to adjust each matrix given by helper_cluster by number of its categories
    # both methods are presented in Cameron, Gelbach and Miller (2011), section 2.3
    # I choose the first option following reghdfe
    G = minimum(nclusters(v))
    rmul!(S, (size(modelmatrix(x), 1) - 1) / dof_residual(x) * G / (G - 1))
end

# res is a Vector in OLS, Matrix in IV
function helper_cluster(X::Matrix, res::Union{Vector, Matrix}, g::GroupedArray)
    X2 = zeros(eltype(X), g.n, size(X, 2) * size(res, 2))
    idx = 0
    for k in 1:size(res, 2)
        for j in 1:size(X, 2)
            idx += 1
            @inbounds @simd for i in 1:size(X, 1)
                X2[g.refs[i], idx] += X[i, j] * res[i, k]
            end
        end
    end
    return Symmetric(X2' * X2)
end

function StatsBase.vcov(x::RegressionModel, v::ClusterCovariance)
    xtx = inv(crossmodelmatrix(x))
    pinvertible(Symmetric(xtx * S_hat(x, v) * xtx))
end


text:
#using Base.Threads

function learn_basic()
    a = zeros(10)
    Threads.@threads for i = 1:10
        a[i] = Threads.threadid()
    end
    @show a

    i = Threads.Atomic{Int}(0)
    ids = zeros(4)
    old_is = zeros(4)
    Threads.@threads for id in 1:4
        old_is[id] = Threads.atomic_add!(i,id) # return old value
        ids[id] = id
    end
    @show old_is
    @show ids
end

function learn_write_read_file()
    n = 5000000
    done = Threads.Atomic{Bool}(false)

    writefile = function()
        open("output.txt", "w") do io
            for i in 1:n
                #sleep(0.1)
                println(io, "Hello! Number $i")
                #flush(io)
            end
        end

        Threads.atomic_or!(done, true)
        println("Done writing")
    end

    readfile = function()
        while !isfile("output.txt")
            sleep(0.1)
        end

        pos = 0 # position of from last file opening
        first_line = ""
        final_line = ""

        while !done[]
            sleep(0.2)
            open("output.txt", "r") do io
                skip(io, pos)

                first_line = readline(io)
                for line in eachline(io)
                    final_line = line
                end

                pos = position(io)
                #println("EOF: $(eof(io))") # EOF: true
            end
            # When first_line is blank meaning no update on the file between
            # opening.
            println("First line: $first_line")
            println("Final line: $final_line")
        end

        # do one more time here?
    end

    write_task = Threads.@spawn writefile()
    read_task = Threads.@spawn readfile()

    wait(write_task) # use fetch() to get return value
    wait(read_task)
end

function learn_error()
    # task = Threads.@spawn (function()
    #     error("Message not printed")
    # end)()
    # wait(task) # error printed because of wait()

    task = Threads.@spawn (function()
        try
            error("Message will not be printed without try/catch")
        catch e
            @error e.msg
            @error stacktrace(catch_backtrace())
            rethrow(e)
        end
    end)()

    println("Will not reach here")
end

function learn_Task()
end
function hello()
    println("First")
    sleep(5)
    println("Second")
end

#learn_basic()
#learn_write_read_file()
learn_error()

nothing


text:
using Test
using Mathematica: weval, buildexpr, getexpr
using MathLink: @W_str

@testset "buildexpr" begin
    @test buildexpr(1) == 1
    @test buildexpr("s") == "s"
    @test buildexpr(:x) == W"x"
    @test buildexpr([1,2,3]) == W"List"(1,2,3)
    @test buildexpr([1 2; 3 4]) == W"List"(W"List"(1, 2), W"List"(3, 4))
    @test buildexpr(1//2) == W"Rational"(1, 2)
    @test buildexpr(:x=>1) == W"Rule"(W"x", 1)
    @test buildexpr(:f(:x, (1, 2))) == buildexpr(:f(:x, [1, 2])) == W"f"(W"x", W"List"(1, 2))
    @test buildexpr(Dict(:a => 1, :b => 2)) == W"Association"(W"Rule"(W"a", 1), W"Rule"(W"b", 2))
end

@testset "getexpr" begin
    @test getexpr(W"List"(1, 2, 3)) == [1, 2, 3]
    @test getexpr(W"List"(W"List"(1, 2), W"List"(3, 4))) == [1 2; 3 4]
    @test getexpr(W"List"(W"List"(1, 2, 3), W"List"(4, 5, 6))) == [1 2 3; 4 5 6]
    arr = reshape(1:12, 2, 3, 2)
    @test getexpr(buildexpr(arr)) == arr
    @test getexpr(W"Rule"(1, W"x")) == (1=>:x)
    @test getexpr(W"Association"(W"Rule"(W"a", 1), W"Rule"(W"b", 2))) == Dict(:a => 1, :b => 2)
end

@testset "operators" begin
    @testset "call" begin
        @test :x(1, :x) == W"x"(1, W"x")
    end
    @testset "plus" begin
        @test :x + 1 == W"Plus"(W"x", 1)
        @test :x + 2 + 3.2 == W"Plus"(W"x", 2, 3.2)
        @test 1 + :Sin(:x) + 3.2 == W"Plus"(1, W"Sin"(W"x"), 3.2)
    end
    @testset "minus" begin
        @test -:x == W"Times"(-1, W"x")
        @test 1 - :x == W"Plus"(1, W"Times"(-1, W"x"))
        @test :x - 1 == W"Plus"(W"x", -1)
        @test :x - :y - 1 == W"Plus"(W"x", W"Times"(-1, W"y"), -1)
        @test 1 - (:x + 2) == W"Plus"(1, W"Times"(-1, W"Plus"(W"x", 2)))
    end
    @testset "times" begin
        @test 2 * :x * :y == W"Times"(2, W"x", W"y")
    end
    @testset "divide" begin
        @test 1 / :x == W"Times"(1, W"Power"(W"x", -1))
        @test :x / 3 == W"Times"(W"x", W"Power"(3, -1))
    end
    @testset "power" begin
        @test :a ^ :b ^ :c == W"Power"(W"a", W"Power"(W"b", W"c"))
    end
end

@testset "functions" begin
    @test (1 |> :Simplify) == W"Simplify"(1)
end


text:
function create_splittable_from_params(params::cluster_parameters, α::Float32)
    params_l = deepcopy(params)
    params_l.distribution = sample_distribution(params.posterior_hyperparams)
    params_r = deepcopy(params)
    params_r.distribution = sample_distribution(params.posterior_hyperparams)
    lr_weights = rand(Dirichlet(Float64.([α / 2, α / 2])))
    return splittable_cluster_params(
        params,
        params_l,
        params_r,
        lr_weights,
        false,
        ones(burnout_period + 5) * -Inf,
    )
end

function merge_clusters_to_splittable(
    cpl::cluster_parameters,
    cpr::cluster_parameters,
    α::Float32,
)
    suff_stats = aggregate_suff_stats(cpl.suff_statistics, cpr.suff_statistics)
    posterior_hyperparams = calc_posterior(cpl.hyperparams, suff_stats)
    lr_weights = rand(Dirichlet(Float64.([
        cpl.suff_statistics.N + (α / 2),
        cpr.suff_statistics.N + (α / 2),
    ])))
    cp = cluster_parameters(
        cpl.hyperparams,
        cpl.distribution,
        suff_stats,
        posterior_hyperparams,
    )
    return splittable_cluster_params(
        cp,
        cpl,
        cpr,
        lr_weights,
        false,
        ones(burnout_period + 5) * -Inf,
    )
end

function should_merge!(
    should_merge::AbstractArray{Float32,1},
    cpl::cluster_parameters,
    cpr::cluster_parameters,
    α::Float32,
    final::Bool,
)
    new_suff = aggregate_suff_stats(cpl.suff_statistics, cpr.suff_statistics)
    cp = cluster_parameters(
        cpl.hyperparams,
        cpl.distribution,
        new_suff,
        cpl.posterior_hyperparams,
    )
    cp.posterior_hyperparams = calc_posterior(cp.hyperparams, cp.suff_statistics)
    log_likihood_l = log_marginal_likelihood(
        cpl.hyperparams,
        cpl.posterior_hyperparams,
        cpl.suff_statistics,
    )
    log_likihood_r = log_marginal_likelihood(
        cpr.hyperparams,
        cpr.posterior_hyperparams,
        cpr.suff_statistics,
    )
    log_likihood = log_marginal_likelihood(
        cp.hyperparams,
        cp.posterior_hyperparams,
        cp.suff_statistics,
    )
    log_HR = (
        -log(α) + logabsgamma(α)[1] - 2 * logabsgamma(0.5 * α)[1] +
        logabsgamma(cp.suff_statistics.N)[1] -
        logabsgamma(cp.suff_statistics.N + α)[1] +
        logabsgamma(cpl.suff_statistics.N + 0.5 * α)[1] -
        logabsgamma(cpl.suff_statistics.N)[1] - logabsgamma(cpr.suff_statistics.N)[1] +
        logabsgamma(cpr.suff_statistics.N + 0.5 * α)[1] +
        log_likihood - log_likihood_l - log_likihood_r
    )
    # log_HR = -(log(α) +
    #     lgamma(cpl.suff_statistics.N) + log_likihood_l +
    #     lgamma(cpr.suff_statistics.N) + log_likihood_r -
    #     (lgamma(cp.suff_statistics.N) + log_likihood))
    if (log_HR > log(rand())) || (final && log_HR > log(0.1))
        should_merge .= 1
    end
end

function sample_cluster_params(params::splittable_cluster_params, α::Float32, first::Bool)
    points_count = Vector{Float32}()
    params.cluster_params.distribution = sample_distribution(
        first ? params.cluster_params.hyperparams :
        params.cluster_params.posterior_hyperparams,
    )
    params.cluster_params_l.distribution = sample_distribution(
        first ? params.cluster_params_l.hyperparams :
        params.cluster_params_l.posterior_hyperparams,
    )
    params.cluster_params_r.distribution = sample_distribution(
        first ? params.cluster_params_r.hyperparams :
        params.cluster_params_r.posterior_hyperparams,
    )
    push!(points_count, params.cluster_params_l.suff_statistics.N)
    push!(points_count, params.cluster_params_r.suff_statistics.N)
    points_count .+= α / 2
    params.lr_weights = rand(Dirichlet(Float64.(points_count)))

    log_likihood_l = log_marginal_likelihood(
        params.cluster_params_l.hyperparams,
        params.cluster_params_l.posterior_hyperparams,
        params.cluster_params_l.suff_statistics,
    )
    log_likihood_r = log_marginal_likelihood(
        params.cluster_params_r.hyperparams,
        params.cluster_params_r.posterior_hyperparams,
        params.cluster_params_r.suff_statistics,
    )

    params.logsublikelihood_hist[1:burnout_period-1] =
        params.logsublikelihood_hist[2:burnout_period]
    params.logsublikelihood_hist[burnout_period] = log_likihood_l + log_likihood_r
    logsublikelihood_now = 0.0
    for i = 1:burnout_period
        logsublikelihood_now +=
            params.logsublikelihood_hist[i] * (1 / (burnout_period - 0.1))
    end
    if logsublikelihood_now != -Inf &&
       logsublikelihood_now - params.logsublikelihood_hist[burnout_period] < 1e-2 # propogate abs change to other versions?
        # println(params.logsublikelihood_hist)
        params.splittable = true
    end

    return params
end


text:
# https://github.com/joshday/OnlineStats.jl/blob/master/notebooks/animations.jl
#using UnicodePlots
using Flux
using Printf
using Plots
using Distributions: Uniform, Normal
using Statistics: mean
using Base.Iterators: partition
using Random: randperm, seed!
using SpikingNeuralNetworks
SNN = SpikingNeuralNetworks
SNN.@load_units
#using MLJ
using Markdown
using InteractiveUtils
using JLD
using Pkg
using PyCall
using Flux
using PyPlot
using Zygote
##
# add zygote gradient tracking.
##
# https://discourse.julialang.org/t/multilayer-perceptron-with-multidimensional-output-array-using-flux/69747/2
dir(x) = fieldnames(typeof(x))
singlify(x::Float64) = Float32(x)
predict(x, m) = m.W*x .+ m.b

filename = string("../JLD/PopulationScatter_adexp.jld")

trace = load(filename , "trace")
evo_population = [t.metadata["pop"] for t in trace]
fitpop = [t.metadata["fitpop"] for t in trace]
singlify(x::Float64) = Float32(x)
x = []#evo_population
y = []#evo_loss
refitpop = []
rex = []
for (i,l) in enumerate(evo_population)
	if i>1
		l[1] = singlify.(l[1])
		append!(x,l[1])
		append!(refitpop,fitpop[i])

	end
end

y = refitpop
long_error=[]
for f in refitpop
	for i in f
		append!(long_error,i[1])
	end
end
y = singlify.(y)
data = Flux.Data.DataLoader((x,y), batchsize=4)#,shuffle=true);
#data = Flux.Data.DataLoader((x,y), batchsize=4,shuffle=true); # automatic batching convenience
X1, Y1 = data#first(data);
@show size(X1) size(Y1)
#X2, Y2 = data
#@show size(X1) size(Y1)

Nx = size(X1, 1)
Ny = size(Y1, 1)
model = Chain(Dense(Nx, Nx, sigmoid), Dense(Nx, Ny, sigmoid), Dense(Ny, Ny, identity));
#gs = gradient(() -> sum(model(X1)), params(model))
# Loss/cost function
loss(x, y) = Flux.mse(model(x), y)
loss(X1,Y1)
# Parameters of the model
para = Flux.params(model);
trnloss = [];
opt = ADAM()

for epoch = 1:500
	@show(size(data))
	@show(size(X1))

    Flux.train!(loss,para,data,opt)
    append!(trnloss, loss(X1, Y1))
	@show(trnloss)
end

amod = model(last(keys(opt.state)))
amod.W
predict(x, m) = m.W*x .+ m.b

gs[amod.W]

mach = machine(model, X1, Y1)
evaluate!(mach, resampling=cv, measure=l2, verbosity=0)

ALLEN_DURATION = 2000 * ms
ALLEN_DELAY = 1000 * ms

function nloss(E,ngt_spikes,ground_spikes)
    spikes = get_spikes(E)
    spikes = [s/1000.0 for s in spikes]

	#opt_vec = [i[1] for i in opt_vec]
	#s_a = signal(opt_vec, length(opt_vec)/last(vmgtt))
	#s_b = signal(vmgtv, length(vmgtt)/last(vmgtt))

	maxt = findmax(sort!(unique(vcat(spikes,ground_spikes))))[1]
    if size(spikes)[1]>1
        t, S = SPIKE_distance_profile(spikes, ground_spikes;t0=0,tf = maxt)
        spkdistance = SpikeSynchrony.trapezoid_integral(t, S)/(t[end]-t[1]) # == SPIKE_distance(y1, y2)
    else
        spkdistance = 10.0
    end
	if length(spikes)>1

		custom_raster2(spikes,ground_spikes)
		custom_raster(spikes,ground_spikes)

		#savefig(crp, "aligned_VMs_adexp.png")
		#display(crp)
	end
	spkdistance*=spkdistance

    delta = abs(size(spikes)[1] - ngt_spikes)
    return spkdistance+delta

end

function nloss(param)
    current = current_search(cell_type,param,ngt_spikes)
	if cell_type=="IZHI"
        pp = SNN.IZParameter(;a = param[1], b = param[2], c = param[3], d = param[4])
        E = SNN.IZ(;N = 1, param = pp)
    end
    if cell_type=="ADEXP"
        adparam = SNN.ADEXParameter(;a = param[1],
            b = param[2],
            cm = param[3],
            v0 = param[4],
            τ_m = param[5],
            τ_w = param[6],
            θ = param[7],
            delta_T = param[8],
            v_reset = param[9],
            spike_delta = param[10])

        E = SNN.AD(;N = 1, param=adparam)
    end

    SNN.monitor(E, [:v])
    SNN.monitor(E, [:fire])
    ALLEN_DURATION = 2000 * ms
    ALLEN_DELAY = 1000 * ms

    E.I = [current*nA]
	#@show(simdur)
    SNN.sim!([E]; dt =1*ms, delay=ALLEN_DELAY,stimulus_duration=ALLEN_DURATION,simulation_duration = simdur)
    vecplot(E, :v) |> display
    error = loss(E,ngt_spikes,ground_spikes)
	@show(error)

    error
end


text:
@struct(type NameStr
 ntype::Cshort # VARIABLE TYPE: 1=NUMERIC, 2=CHAR
 nhfun::Cshort # HASH OF NNAME (always 0)
 nlng::Cshort # LENGTH OF VARIABLE IN OBSERVATION
 nvar0::Cshort # VARNUM
 nname::ASCIIString(8); # NAME OF VARIABLE
 nlabel::ASCIIString(40) # LABEL OF VARIABLE
 nform::ASCIIString(8) # NAME OF FORMAT
 nfl::Cshort # FORMAT FIELD LENGTH OR 0
 nfd::Cshort # FORMAT NUMBER OF DECIMALS
 nfj::Cshort # 0=LEFT JUSTIFICATION, 1=RIGHT JUST
 nfill::Array{Cchar, 1}(2) # (UNUSED, FOR ALIGNMENT AND FUTURE)
 niform::ASCIIString(8) # NAME OF INPUT FORMAT
 nifl::Cshort # INFORMAT LENGTH ATTRIBUTE
 nifd::Cshort # INFORMAT NUMBER OF DECIMALS
 npos::Int32 # POSITION OF VALUE IN OBSERVATION
 rest::ASCIIString(48) # remaining fields are irrelevant. there are an additional
end,
align_packed,
:BigEndian)
function show(io::IO, ns::NameStr)
    print(io, "Name: ", strip(ns.nname)", ")
    print(io, "Label: ", strip(ns.nlabel)", ")
    print(io, ns.ntype == 1? "Numeric, " : "Character, ")
    print(io, "Pos: ", ns.npos, ", ")
    print(io, "Len: ", ns.nlng)
end

const NAMESTR_SIZE = calcsize(NameStr)


text:
@userplot vfpoil

@recipe function f(h::vfpoil;ϵs=0.0006, Tss=80, ∇Ts=1, GORs=false, dis=2.99, Qns=10, Dns=50, tols=0.05, sgs=0.6)

        Depth,Thp,Bsw,OilPVT,GasPVT,WaterPVT, Qmax=h.args

DepthRange, Qrange, Pwf, l, P=OilOutflowSen(Depth,Thp,Bsw,OilPVT,GasPVT,WaterPVT, Qmax;ϵ=ϵs, Ts=Tss, ∇T=∇Ts, GOR=GORs, di=dis, Qn=Qns, Dn=Dns, tol=tols, sg=sgs)

    xlabel --> "Rate [bbl/d]"
    ylabel --> "Pwf [psi]"
    xformatter --> :plain
    label --> l
    @series begin
        seriestype := :path
        Qrange, Pwf'
    end
    end

@userplot vfpgas

@recipe function f(h::vfpgas; dis=2.99, Tss=80, ∇Ts=1, ϵs=0.0006, Qns=10, Dns=50, tols=0.05, sgs=0.6, Incs=0)

    Depth,Thp,GasPVT, Qmax=h.args

DepthRange, Qrange, Pwf, leg, P = GasOutflow(Depth,Thp,GasPVT, Qmax; di=dis, Ts=Tss, ∇T=∇Ts, ϵ=ϵs, Qn=Qns, Dn=Dns, tol=tols, sg=sgs, Inc=Incs)

xlabel --> "Rate [Mscfd]"
ylabel --> "Pwf [psi]"
xformatter --> :plain
label --> leg
@series begin
    seriestype := :path
    Qrange, Pwf'
end
end


text:
using RecipesPipeline
using Test

makie_test_dir = joinpath(@__DIR__, "test_makie")
mkpath(makie_test_dir)

@testset "RecipesPipeline.jl" begin
    @testset "Makie integration" begin cd(makie_test_dir) do
        include("makie.jl")
    end end
end


text:
importall Bukdu

type PageController <: ApplicationController
end

type UserController <: ApplicationController
end

index(::PageController) = "index page"
show(c::PageController) = parse(Int, c[:params]["id"])
delete(::PageController) = nothing

index(::UserController) = "index user"
show(c::UserController) = c[:params]

Router() do
    resources("/pages", PageController) do
        resources("/users", UserController)
    end
end

Logger.set_level(false)
using Base.Test
conn = (Router)(get, "/pages")
@test conn.status == 200
@test conn.resp_body == "index page"

conn = (Router)(get, "/pages/1")
@test conn.status == 200
@test conn.resp_body == 1
@test conn.params["id"] == "1"

conn = (Router)(get, "/pages/1/users")
@test conn.status == 200
@test conn.resp_body == "index user"
@test conn.params["page_id"] == "1"
@test !haskey(conn.params, "id")

conn = (Router)(get, "/pages/1/users/2")
@test conn.status == 200
@test conn.resp_body == Assoc(page_id="1", id="2")
@test conn.params["page_id"] == "1"
@test conn.params["id"] == "2"

conn = (Router)(delete, "/pages/1")
@test conn.status == 200

@test_throws NoRouteError (Router)(get, "/unknown/1")

reset(Router)

@test_throws NoRouteError (Router)(get, "/pages")

Router() do
    resources("/pages", PageController, except=[delete], singleton=true)
    resources("/users", UserController, only=[index], singleton=true)
end

conn = (Router)(get, "/pages")
@test conn.status == 200

@test_throws NoRouteError (Router)(delete, "/pages/1")

conn = (Router)(get, "/users")
@test conn.status == 200

@test_throws NoRouteError (Router)(delete, "/users/1")


text:
using DashBootstrapComponents, DashHtmlComponents

button = html_div([
    dbc_button("Click me", id = "example-button", className = "me-2", n_clicks = 0),
    html_span(id = "example-output", style = Dict("verticalAlign" => "middle")),
]);

callback!(
    app,
    Output("example-output", "children"),
    Input("example-button", "n_clicks"),
) do n
    return n > 0 ? "Clicked $n times." : "Not clicked."
end;

max_stars_repo_name: glsdown/dash-bootstrap-components

max_issues_repo_name: glsdown/dash-bootstrap-components

max_forks_repo_name: glsdown/dash-bootstrap-components


text:
module _PhasEn
export PhasEn
using Plots
using Random: randperm
    """
        Phas = PhasEn(Sig)

    Returns the phase entropy (`Phas`) estimate of the data sequence (`Sig`)
    using the default parameters:
    angular partitions = 4, time delay = 1, logarithm = natural,

        Phas = PhasEn(Sig::AbstractArray{T,1} where T<:Real; K::Int=4, tau::Int=1, Logx::Real=exp(1), Norm::Bool=true, Plotx::Bool=false)

    Returns the phase entropy (`Phas`) estimate of the data sequence (`Sig`)
    using the specified 'keyword' arguments:

    # Arguments:
    `K`     - Angular partitions (coarse graining), an integer > 1  \n
                *Note: Division of partitions begins along the positive x-axis. As this point is somewhat arbitrary, it is
                 recommended to use even-numbered (preferably multiples of 4) partitions for sake of symmetry. \n
    `tau`   - Time Delay, a positive integer    \n
    `Logx`  - Logarithm base, a positive scalar \n
    `Norm`  - Normalisation of `Phas` value:    \n
              [false]  no normalisation
              [true]   normalises w.r.t. the number of partitions Log(`K`)
    `Plotx` - When `Plotx` == true, returns Poincaré plot (default: false)  \n

    # See also `SampEn`, `ApEn`, `GridEn`, `MSEn`, `SlopEn`, `CoSiEn`, `BubbEn`

    # References:
        [1] Ashish Rohila and Ambalika Sharma,
            "Phase entropy: a new complexity measure for heart rate
            variability."
            Physiological measurement
            40.10 (2019): 105006.

    """
    function PhasEn(Sig::AbstractArray{T,1} where T<:Real; K::Int=4, tau::Int=1,
        Logx::Real=exp(1), Norm::Bool=true, Plotx::Bool=false)

    N = size(Sig,1)
    (N > 10) ? nothing : error("Sig:   must be a numeric vector")
    (K > 1) ? nothing :  error("K:     must be an integer > 1")
    (tau >0) ? nothing : error("tau:   must be an integer > 0")
    (Logx>0) ? nothing : error("Logx:  must be a positive number > 0")

    Yn = Sig[1+2*tau:end] .- Sig[tau+1:end-tau]
    Xn = Sig[tau+1:end-tau] .- Sig[1:end-2*tau]
    Theta_r = atan.(Yn./Xn)
    Theta_r[(Yn.<0) .& (Xn.<0)] .+= pi
    Theta_r[(Yn.<0) .& (Xn.>0)] .+= 2*pi
    Theta_r[(Yn.>0) .& (Xn.<0)] .+= pi

    Limx = ceil.(maximum(abs.(vcat(Yn, Xn))))
    Angs = range(0,2*pi,length=K+1)
    Tx = zeros(Int, K, length(Theta_r))
    Si = zeros(K)
    for n = 1:K
        Temp = (Theta_r .> Angs[n]) .& (Theta_r .< Angs[n+1]);
        Tx[n,Temp] .= 1
        Si[n] = sum(Theta_r[Temp])
    end
    Si = Si[Si.!=0]
    Phas = -sum((Si./sum(Si)).*(log.(Logx, Si./sum(Si))))
    if Norm
        Phas = Phas/(log(Logx, K))
    end

    if Plotx
        Ys = sin.(Angs)*Limx*sqrt(2);
        Xs = cos.(Angs)*Limx*sqrt(2);
        Cols = hcat(zeros(K), repeat(randperm(K)/K,outer=(1,2)))
        Tx = Bool.(Tx)
        xx = plot()
        for n = 1:K
            plot!(Xn[Tx[n,:]], Yn[Tx[n,:]], seriestype =:scatter,
                markersize = 2,
                markercolor=RGB(Cols[n,:]...,),
                markerstrokecolor=RGB(Cols[n,:]...,))
        end
        #plot!(hcat(zeros(K+1), Xs), hcat(Ys,zeros(K+1)),c=:magenta,
        plot!(hcat(Xs,zeros(K+1))', hcat(Ys,zeros(K+1))',c=:magenta,
            xlim = (-Limx, Limx), ylim = (-Limx, Limx),
            size = (400, 400), legend=false,
            xticks = [-Limx, 0, Limx],  yticks = [-Limx, 0 ,Limx],
            grid = false)
        xlabel!("X ₙ₊ₜ - X ₙ")
        ylabel!("X ₙ₊₂ₜ - X ₙ₊ₜ")
        display(xx)
    end

    return Phas
    end

end

"""
Copyright 2021 Matthew W. Flood, EntropyHub

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

     http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.

For Terms of Use see https://github.com/MattWillFlood/EntropyHub

"""


text:
using CSV
using CairoMakie
using Colors

color_lst = [colorant"#332288" colorant"#555555" colorant"#117733" colorant"#882255"]

f = CSV.File("SamplingData.csv")

dist1 = f.dist1
dist2 = f.dist2
dist3 = f.dist3

noto_sans = "../resources/NotoSans-Regular.ttf"

tickfontsize = 28
labelfontsize = 34

f = Figure(resolution = (750,450), font=noto_sans)

ax = Axis(  f[1,1],
            xticklabelsize=tickfontsize,
            yticklabelsize=tickfontsize,
            yticklabelpad=2,
            xlabel="Sample Count", ylabel = "E[Min Dist]",
            xlabelsize=labelfontsize,
            ylabelsize=labelfontsize,
            yscale = log10,
            xscale = log10,
            yminorticksvisible = true,
            yminorgridvisible = true,
            yminorticks = IntervalsBetween(10),
            xminorticksvisible = true,
            xminorgridvisible = true,
            xminorticks = IntervalsBetween(10))

l1 = lines!(ax, 1:length(dist1), dist1, color = color_lst[1], linewidth=3)
l2 = lines!(ax, 1:length(dist2), dist2, color = color_lst[2], linewidth=3)
l3 = lines!(ax, 1:length(dist3), dist3, color = color_lst[3], linewidth=3)

l1elem = LineElement(color=color_lst[1], linewidth=10)
l2elem = LineElement(color=color_lst[2], linewidth=10)
l3elem = LineElement(color=color_lst[3], linewidth=10)

xlims!(ax,(1,length(dist1)))

Legend(f[1,2], [l3elem,l2elem,l1elem], ["3D", "2D", "1D"], labelsize=labelfontsize)

save("SamplingPlot.pdf",f)


text:
using JuliaParser

include(joinpath(dirname(@__FILE__),"../src/interactiveutil.jl"))

function dotest(input)
lines = split(input,'\n')
code = String(lines[1])
ts = Lexer.TokenStream{Lexer.SourceLocToken}(code)
res = Parser.parse(ts)
textrep = create_textrep(res.loc, code)
if strip(input) != strip(textrep)
    println("GOT:")
    println(textrep)
    error()
end
end

"""
f()
x
xxx
""" |> dotest

"""
t...
x
xxxx
""" |> dotest

"""
for i = 1:10; end
        x xx
    x   xxxx
    xxxxxxxx  xxx
xxxxxxxxxxxxxxxxx
""" |> dotest


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct OperationStatusResponse <: SwaggerModel
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    status::Any # spec type: Union{ Nothing, String } # spec name: status
    startTime::Any # spec type: Union{ Nothing, DateTime } # spec name: startTime
    endTime::Any # spec type: Union{ Nothing, DateTime } # spec name: endTime
    error::Any # spec type: Union{ Nothing, ApiError } # spec name: error

    function OperationStatusResponse(;name=nothing, status=nothing, startTime=nothing, endTime=nothing, error=nothing)
        o = new()
        validate_property(OperationStatusResponse, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(OperationStatusResponse, Symbol("status"), status)
        setfield!(o, Symbol("status"), status)
        validate_property(OperationStatusResponse, Symbol("startTime"), startTime)
        setfield!(o, Symbol("startTime"), startTime)
        validate_property(OperationStatusResponse, Symbol("endTime"), endTime)
        setfield!(o, Symbol("endTime"), endTime)
        validate_property(OperationStatusResponse, Symbol("error"), error)
        setfield!(o, Symbol("error"), error)
        o
    end
end # type OperationStatusResponse

const _property_map_OperationStatusResponse = Dict{Symbol,Symbol}(Symbol("name")=>Symbol("name"), Symbol("status")=>Symbol("status"), Symbol("startTime")=>Symbol("startTime"), Symbol("endTime")=>Symbol("endTime"), Symbol("error")=>Symbol("error"))
const _property_types_OperationStatusResponse = Dict{Symbol,String}(Symbol("name")=>"String", Symbol("status")=>"String", Symbol("startTime")=>"DateTime", Symbol("endTime")=>"DateTime", Symbol("error")=>"ApiError")
Base.propertynames(::Type{ OperationStatusResponse }) = collect(keys(_property_map_OperationStatusResponse))
Swagger.property_type(::Type{ OperationStatusResponse }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_OperationStatusResponse[name]))}
Swagger.field_name(::Type{ OperationStatusResponse }, property_name::Symbol) =  _property_map_OperationStatusResponse[property_name]

function check_required(o::OperationStatusResponse)
    true
end

function validate_property(::Type{ OperationStatusResponse }, name::Symbol, val)
end


text:
module Mergers

using Flux
import Base: show

using ..AbstractMergers

export Merger, inputname

"""
    Merger{F,O}

An element in a `FeedbackChain` in which the forward stream and feedback stream
are combined according to an operation `op`.

# Fields
- `splitname::String`: name of the `Splitter` node from which the feedback is taken
- `fb::F`: feedback branch
- `op::O`: operation to combine forward and feedback branches

# Details
`fb` typically takes the form of a Flux operation or chain. When a `FeedbackChain`
encounters a `Merger`, it will look up the state `s` of the `Splitter` given by `forkname`
from the previous timestep, apply `fb` to it and combine it with the forward input `x`
according  to `op(x, fb(s))`
"""
struct Merger{F,O} <: AbstractMerger
    splitname::String
    fb::F
    op::O
end # struct Merger

function (m::Merger)(x, y)
    m.op(x, m.fb(y[m.splitname]))
end

# These overloads ensure that a Merger behaves as Flux expects, e.g.,
# when moving to gpu or collecting parameters.
Flux.@functor Merger
#children(m::Merger) = (m.fb, m.op)
#mapchildren(f, m::Merger) = Merger(m.splitname, f(m.fb), m.op)

"""
    inputname(m::Merger)

Return the name of the `Splitter` from which `m` gets its input.
"""
inputname(m::Merger) = m.splitname

function show(io::IO, m::Merger)
    print(io, "Merger(\"", m.splitname, "\", ", m.fb, ", ", m.op, ")")
end # function show
end # module Mergers


text:
using PrettyPrint
using Test
struct Account
    username :: String
    password :: String
end

struct S1
  i :: Int
  f :: Float64
end

struct S2
  s  :: Vector{String}
  s1 :: S1
end

@testset "PrettyPrint.jl" begin

    data = S2(
        ["114514", "as we can"],
        S1(42, 9.96)
    )

    @test pformat(data) == """S2(
  s=[
    "114514",
    "as we can",
  ],
  s1=S1(
    i=42,
    f=9.96,
  ),
)"""

    @testset "extension" begin
    @test pformat([Account("van", "gd"), Account("thautwarm", "996icu")]) == """[
  Account(
    username="van",
    password="gd",
  ),
  Account(
    username="thautwarm",
    password="996icu",
  ),
]"""
    PrettyPrint.pprint_impl(io, account::Account, indent::Int, newline::Bool) = print(io, "Account($(account.username))")
    @test pformat([Account("van", "gd"), Account("thautwarm", "996icu")]) == """[
  Account(van),
  Account(thautwarm),
]"""
    @test pformat(Set([Account("van", "gd")])) == """{
  Account(van),
}"""
    @test pformat((Account("van", "gd"), Account("thautwarm", "996icu"))) == """(
  Account(van),
  Account(thautwarm),
)"""
    end

end


text:
## Example script for linear-strength vortex panel method
using LinearAlgebra
using Base.Iterators
using Seaborn
using AeroMDAO

## Airfoil
alpha_u = [0.2, 0.3, 0.2, 0.15, 0.2]
alpha_l = [-0.2, -0.1, -0.1, -0.001]
dzs     = (0., 0.)
# airfoil = (Foil ∘ kulfan_CST)(alpha_u, alpha_l, dzs, 0.0, 60);
airfoil = Foil(naca4((0,0,1,2), 81; sharp_trailing_edge = false))
V, α    = 1., 1.
ρ       = 1.225
uniform = Uniform2D(V, α)
num_pans = 80
panels  = paneller(airfoil, num_pans);

## Vortex panel
A = @time vortex_influence_matrix(panels)
b = @time neumann_boundary_condition(panels, velocity(uniform))
γs = A \ b

##
cl = 2 * sum((fwdsum(γs) / 2) .* panel_length.(panels)) / uniform.mag

## Velocities and pressure coefficient
# panel_vels = [ velocity(uniform) .+ sum(vortex_velocity.(γs[1:end-1], γs[2:end], panels, Ref(panel_i))) for panel_i in panels ]

qts = γs
cps = @. 1 - qts^2 / uniform.mag^2

## Pressure coefficient
# upper, lower = get_surface_values(panels, cps)
# upper, lower = [ upper; lower[1] ], [ lower; upper[1] ]
# plot(first.(upper), last.(upper), label = "Upper")
# plot(first.(lower), last.(lower), label = "Lower")

plot(first.(panel_points(panels)), cps)
ylim(maximum(cps), minimum(cps))
xlabel("(x/c)")
ylabel("Cp")
legend()
show()

## Plotting
x_domain, y_domain = (-1, 2), (-1, 1)
grid_size = 50
x_dom, y_dom = range(x_domain..., length = grid_size), range(y_domain..., length = grid_size)
grid = product(x_dom, y_dom)
pts = panel_points(panels);

vortex_vels = [ velocity(uniform) .+ sum(vortex_velocity.(γs[1:end-1], γs[2:end], panels, x, y)) for (x, y) in grid ]

speeds = @. sqrt(first(vortex_vels)^2 + last(vortex_vels)^2)
cps = @. 1 - speeds^2 / uniform.mag^2

##
contourf(first.(grid), last.(grid), cps)
CB1 = colorbar(label = "Pressure Coefficient (Cp)")
# quiver(first.(grid), last.(grid), first.(vortex_vels), last.(vortex_vels), speeds)
streamplot(first.(grid)', last.(grid)', first.(vortex_vels)', last.(vortex_vels)', color = speeds', cmap = "coolwarm", density = 2)
CB2 = colorbar(orientation="horizontal", label = "Relative Airspeed (U/U∞)")
fill(first.(pts), last.(pts), color = "black", zorder = 3)
tight_layout()
show()

max_forks_repo_name: GodotMisogi/AeroMDAO


text:
using Test
using TextAnalysis
using TextModels

println("Running tests:")

include("crf.jl")
include("ner.jl")
include("pos.jl")
include("averagePerceptronTagger.jl")
include("ulmfit.jl")
include("sentiment.jl")