text:
"""
This file initializes the ClimateERA module by setting and determining the
ECMWF reanalysis parameters to be analyzed and the regions upon which the data
are to be extracted from.  Functionalities include:
    - Setting up of reanalysis module type
    - Setting up of reanalysis parameters to be analyzed
    - Setting up of time steps upon which data are to be downloaded
    - Setting up of region of analysis based on ClimateEasy

"""

# ClimateERA Module Setup

function eramoduledisp(init::Dict)

    if init["actionID"] == 1; len = 4; elseif init["actionID"] == 2; len = 6; end
    @debug "$(Dates.now()) - There are $(len) types of modules that ClimateERA can $(init["action"])."

    @debug "$(Dates.now()) - 1) Dry Surface Modules    (e.g. Surface Winds)"
    @debug "$(Dates.now()) - 2) Dry Pressure Modules   (e.g. Winds at Pressure Height)"
    @debug "$(Dates.now()) - 3) Moist Surface Modules  (e.g. Rainfall, Total Column Water)"
    @debug "$(Dates.now()) - 4) Moist Pressure Modules (e.g. Humidity at Pressure Height)"

    if init["actionID"] == 2
        @debug "$(Dates.now()) - 5) Calc Surface Modules   (e.g. PI)"
        @debug "$(Dates.now()) - 6) Calc Pressure Modules  (e.g. Eddy Kinetic Energy, Psi)"
    end

    return len

end

# ClimateERA Parameter Setup

function eraparametercopy(;overwrite::Bool=false)

    jfol = joinpath(DEPOT_PATH[1],"files/ClimateERA/"); mkpath(jfol);
    ftem = joinpath(@__DIR__,"../extra/epartemplate.txt")
    fpar = joinpath(jfol,"eraparameters.txt")

    if !overwrite
        if !isfile(fpar)

            @debug "$(Dates.now()) - Unable to find eraparameters.txt, copying data from epartemplate.txt ..."

            open(fpar,"w") do io
                open(ftem) do f
                    for line in readlines(f)
                        write(io,"$line\n")
                    end
                end
            end

        end
    else

        if isfile(fpar)
            @warn "$(Dates.now()) - Overwriting eraparameters.txt in $jfol ..."
            rm(fpar,force=true)
        end

        open(fpar,"w") do io
            open(ftem) do f
                for line in readlines(f)
                    write(io,"$line\n")
                end
            end
        end

    end

    return fpar

end

function eraparameterload()

    @debug "$(Dates.now()) - Loading information on parameters used in ERA reanalysis."
    return readdlm(eraparametercopy(),',',comments=true);

end

function eraparameterload(init::Dict)

    @debug "$(Dates.now()) - Loading information on parameters used in ERA reanalysis."
    allparams = readdlm(eraparametercopy(),',',comments=true);

    @debug "$(Dates.now()) - Filtering out for parameters in the $(init["modulename"]) module."
    parmods = allparams[:,1]; return allparams[(parmods.==init["moduletype"]),:];

end

function eraparameterdisp(parlist::AbstractArray,init::Dict)
    @info "$(Dates.now()) - The following variables are offered in the $(init["modulename"]) module:"
    for ii = 1 : size(parlist,1); @info "$(Dates.now()) - $(ii)) $(parlist[ii,6])" end
end

function eraparameteradd(fadd::AbstractString)

    if !isfile(fadd); error("$(Dates.now()) - The file $(fadd) does not exist."); end
    ainfo = readdlm(fadd,',',comments=true); aparID = ainfo[:,2]; nadd = length(aparID);

    for iadd = 1 : nadd
        eraparameteradd(modID=ainfo[iadd,1],parID=ainfo[iadd,2],ncID=ainfo[iadd,3],
                        era5=ainfo[iadd,4],erai=ainfo[iadd,5],
                        full=ainfo[iadd,6],unit=ainfo[iadd,7],throw=false);
    end

end

function eraparameteradd(;
    modID::AbstractString, parID::AbstractString, ncID::AbstractString,
    era5::AbstractString, erai::AbstractString,
    full::AbstractString, unit::AbstractString,
    throw::Bool=true
)
    fpar = eraparametercopy(); pinfo = eraparameterload(); eparID = pinfo[:,2];

    if sum(eparID.==parID) > 0

        if throw
            error("$(Dates.now()) - Parameter ID already exists.  Please choose a new parID.")
        else
            @info "$(Dates.now()) - $(parID) has already been added to eraparameters.txt"
        end

    else

        open(fpar,"a") do io
            writedlm(io,[modID parID ncID era5 erai full unit],',')
        end

    end

end

function erapressure(emod::Dict)

    if occursin("pre",emod["moduleID"])
        @debug "$(Dates.now()) - A pressure module was selected, and therefore all available pressure levels will be saved into the parameter Dictionary."
        emod["levels"] = erapressureload();
    else
        @debug "$(Dates.now()) - A surface module was selected, and therefore we will save 'sfc' into the parameter level Dictionary."
        emod["levels"] = ["sfc"];
    end

    return

end

function erapressureload()
    return [1,2,3,5,7,10,20,30,50,70,100,125,150,175,200,
            225,250,300,350,400,450,500,550,600,650,700,750,
            775,800,825,850,875,900,925,950,975,1000]
end

# ClimateERA Region Setup

function eraregionload(gregID::AbstractString,init::Dict)

    @info "$(Dates.now()) - Loading available GeoRegions ..."
    geo = GeoRegion(gregID); greggrid = [geo.N,geo.S,geo.E,geo.W]
    if gregID == "GLB"; regglbe = true; else; regglbe = false; end

    @info "$(Dates.now()) - Storing GeoRegion properties and information for the $(geo.name) region ..."
    return Dict("region"=>gregID,"grid"=>greggrid,"name"=>geo.name,"isglobe"=>regglbe)

end

function eraregionstep(gregID::AbstractString,step::Real=0)

    @debug "$(Dates.now()) - Determining spacing between grid points in the GeoRegion ..."
    if step == 0
        if gregID == "GLB";
              step = 1.0;
        else; step = 0.25;
        end
    else
        if !checkegrid(step)
            error("$(Dates.now()) - The grid resolution specified is not valid.")
        end
    end

    return step

end

function eraregionvec(ereg::Dict,emod::Dict,step::Real)

    step = eraregionstep(ereg["region"],step); ereg["step"] = step; N,S,E,W = ereg["grid"]

    @info "$(Dates.now()) - Creating longitude and latitude vectors for the GeoRegion ..."
    lon = convert(Array,W:step:E); if (mod(E,360) == mod(W,360)) && (E!=W); pop!(lon); end
    lat = convert(Array,N:-step:S); nlon = size(lon,1); nlat = size(lat,1);
    ereg["lon"] = lon; ereg["lat"] = lat; ereg["size"] = [nlon,nlat];
    ereg["fol"] = "$(ereg["region"])x$(@sprintf("%.2f",ereg["step"]))"

    return ereg

end

# Initialization

function eramodule(moduleID::AbstractString,init::Dict)

    init["moduletype"] = moduleID; len = eramoduledisp(init);
    if init["actionID"] == 1 && (sum(["csfc","cpre"] .== moduleID) != 0)
        error("$(Dates.now()) - Module ID \"$(moduleID)\" not defined for action \"$(init["action"])\".  Call queryemod(modID=$(moduleID)) for more details.")
    end

    if     moduleID == "dsfc"; init["modulename"] = "dry surface";
    elseif moduleID == "dpre"; init["modulename"] = "dry pressure";
    elseif moduleID == "msfc"; init["modulename"] = "moist surface";
    elseif moduleID == "mpre"; init["modulename"] = "moist pressure";
    elseif moduleID == "csfc"; init["modulename"] = "calc surface";
    elseif moduleID == "cpre"; init["modulename"] = "calc pressure";
    end

    if init["actionID"] == 1 && init["datasetID"] == 1
        if occursin("sfc",moduleID)
            init["moduleprint"] = "reanalysis-era5-single-levels";
        elseif occursin("pre",moduleID)
            init["moduleprint"] = "reanalysis-era5-pressure-levels";
        end
    end

    if occursin("pre",moduleID)
        @debug "$(Dates.now()) - A pressure module was selected, and therefore all available pressure levels will be saved into the parameter Dictionary."
        init["levels"] = erapressureload();
    else
        @debug "$(Dates.now()) - A surface module was selected, and therefore we will save 'sfc' into the parameter level Dictionary."
        init["levels"] = ["sfc"];
    end

    return deepcopy(init)

end

function eraparameter(parameterID::AbstractString,emod::Dict)

    parlist = eraparameterload(emod);

    if sum(parlist[:,2] .== parameterID) == 0
        error("$(Dates.now()) - Invalid parameter choice for $(emod["moduletype"]).  Call queryepar(modID=$(emod["moduletype"]),parID=$(parameterID)) for more information.")
    else
        ID = (parlist[:,2] .== parameterID);
    end

    parinfo = parlist[ID,:];
    @info "$(Dates.now()) - ClimateERA will $(emod["action"]) $(parinfo[6]) data."

    if occursin("sfc",emod["moduletype"])
        return Dict("ID"  =>parinfo[2],"IDnc"=>parinfo[3],
                    "era5"=>parinfo[4],"erai"=>parinfo[5],
                    "name"=>parinfo[6],"unit"=>parinfo[7],"level"=>"sfc");
    else
        return Dict("ID"  =>parinfo[2],"IDnc"=>parinfo[3],
                    "era5"=>parinfo[4],"erai"=>parinfo[5],
                    "name"=>parinfo[6],"unit"=>parinfo[7],"level"=>undef);
    end

end

function eratime(timeID::Integer,emod::Dict)
    if timeID == 0;

        if emod["datasetID"] == 1
              fin = Dates.year(Dates.now())-1;
        else; fin = 2018
        end

        return Dict("Begin"=>1979,"End"=>fin);
        @info "$(Dates.now()) - User has chosen to $(emod["action"]) $(emod["dataset"]) datasets from 1979 to $(fin)."

    else
        return Dict("Begin"=>timeID,"End"=>timeID)
        @info "$(Dates.now()) - User has chosen to $(emod["action"]) $(emod["dataset"]) datasets in $(timeID)."
    end
end

function eratime(timeID::Vector,emod::Dict)
    beg = minimum(timeID); fin = maximum(timeID)
    return Dict("Begin"=>beg,"End"=>fin)
    @info "$(Dates.now()) - User has chosen to $(emod["action"]) $(emod["dataset"]) datasets from $(beg) to $(fin)."
end

function eraregion(gregID::AbstractString,emod::Dict,gres::Real)
    return eraregionvec(eraregionload(gregID,emod),emod,gres)
end

function erainitialize(
    init::Dict;
    modID::AbstractString, parID::AbstractString,
    regID::AbstractString="GLB", timeID::Union{Integer,Vector}=0,
    gres::Real=0
)

    emod = eramodule(modID,init); epar = eraparameter(parID,emod);
    ereg = eraregion(regID,emod,gres); etime = eratime(timeID,emod);

    return emod,epar,ereg,etime

end


text:
using Trixi, OrdinaryDiffEq

dg = DGMulti(polydeg = 3, element_type = Quad(), approximation_type = SBP(),
             surface_integral = SurfaceIntegralWeakForm(FluxLaxFriedrichs()),
             volume_integral = VolumeIntegralFluxDifferencing(flux_ranocha))

equations = CompressibleEulerEquations2D(1.4)
initial_condition = initial_condition_khi

vertex_coordinates, EToV = StartUpDG.uniform_mesh(dg.basis.elementType, 32)
mesh = VertexMappedMesh(vertex_coordinates, EToV, dg, is_periodic=(true, true))
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, dg)

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

summary_callback = SummaryCallback()
alive_callback = AliveCallback(alive_interval=10)
callbacks = CallbackSet(summary_callback,
                        alive_callback)

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

sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt= estimate_dt(mesh, dg), save_everystep=false, callback=callbacks);

summary_callback() # print the timer summary

ext: jl


text:
using Oceananigans.BuoyancyModels
using Oceananigans.Coriolis
using Oceananigans.Operators
using Oceananigans.Operators: ∂xᶠᶜᵃ, ∂yᶜᶠᵃ
using Oceananigans.StokesDrift
using Oceananigans.TurbulenceClosures: ∂ⱼ_τ₁ⱼ, ∂ⱼ_τ₂ⱼ, ∇_dot_qᶜ
using Oceananigans.Advection: div_Uc

using Oceananigans.TurbulenceClosures: shear_production, buoyancy_flux, dissipation, TKETracerIndex
import Oceananigans.TurbulenceClosures: hydrostatic_turbulent_kinetic_energy_tendency

"""
Return the tendency for the horizontal velocity in the x-direction, or the east-west
direction, ``u``, at grid point `i, j, k` for a HydrostaticFreeSurfaceModel.

The tendency for ``u`` is called ``G_u`` and defined via

    ``∂_t u = G_u - ∂_x ϕ_n``

where ϕ_n is the part of the barotropic pressure that's treated
implicitly during time-stepping.
"""
@inline function hydrostatic_free_surface_u_velocity_tendency(i, j, k, grid,
                                                              advection,
                                                              coriolis,
                                                              closure,
                                                              velocities,
                                                              free_surface,
                                                              tracers,
                                                              buoyancy,
                                                              diffusivities,
                                                              hydrostatic_pressure_anomaly,
                                                              auxiliary_fields,
                                                              forcings,
                                                              clock)

    model_fields = merge(hydrostatic_prognostic_fields(velocities, free_surface, tracers), auxiliary_fields)

    return ( - U_dot_∇u(i, j, k, grid, advection, velocities)
             - explicit_barotropic_pressure_x_gradient(i, j, k, grid, free_surface)
             - x_f_cross_U(i, j, k, grid, coriolis, velocities)
             - ∂xᶠᶜᵃ(i, j, k, grid, hydrostatic_pressure_anomaly)
             - ∂ⱼ_τ₁ⱼ(i, j, k, grid, closure, clock, velocities, diffusivities, tracers, buoyancy)
             + forcings.u(i, j, k, grid, clock, hydrostatic_prognostic_fields(velocities, free_surface, tracers)))
end

"""
Return the tendency for the horizontal velocity in the y-direction, or the east-west
direction, ``v``, at grid point `i, j, k` for a HydrostaticFreeSurfaceModel.

The tendency for ``v`` is called ``G_v`` and defined via

    ``∂_t v = G_v - ∂_y ϕ_n``

where ϕ_n is the part of the barotropic pressure that's treated
implicitly during time-stepping.
"""
@inline function hydrostatic_free_surface_v_velocity_tendency(i, j, k, grid,
                                                              advection,
                                                              coriolis,
                                                              closure,
                                                              velocities,
                                                              free_surface,
                                                              tracers,
                                                              buoyancy,
                                                              diffusivities,
                                                              hydrostatic_pressure_anomaly,
                                                              auxiliary_fields,
                                                              forcings,
                                                              clock)

    model_fields = merge(hydrostatic_prognostic_fields(velocities, free_surface, tracers), auxiliary_fields)

    return ( - U_dot_∇v(i, j, k, grid, advection, velocities)
             - explicit_barotropic_pressure_y_gradient(i, j, k, grid, free_surface)
             - y_f_cross_U(i, j, k, grid, coriolis, velocities)
             - ∂yᶜᶠᵃ(i, j, k, grid, hydrostatic_pressure_anomaly)
             - ∂ⱼ_τ₂ⱼ(i, j, k, grid, closure, clock, velocities, diffusivities, tracers, buoyancy)
             + forcings.v(i, j, k, grid, clock, hydrostatic_prognostic_fields(velocities, free_surface, tracers)))
end

"""
Return the tendency for a tracer field with index `tracer_index`
at grid point `i, j, k`.

The tendency is called ``G_c`` and defined via

    ``∂_t c = G_c``

where `c = C[tracer_index]`.
"""
@inline function hydrostatic_free_surface_tracer_tendency(i, j, k, grid,
                                                          val_tracer_index::Val{tracer_index},
                                                          advection,
                                                          closure,
                                                          buoyancy,
                                                          velocities,
                                                          free_surface,
                                                          tracers,
                                                          diffusivities,
                                                          auxiliary_fields,
                                                          forcing,
                                                          clock) where tracer_index

    @inbounds c = tracers[tracer_index]

    model_fields = merge(hydrostatic_prognostic_fields(velocities, free_surface, tracers), auxiliary_fields)

    return ( - div_Uc(i, j, k, grid, advection, velocities, c)
             - ∇_dot_qᶜ(i, j, k, grid, closure, c, val_tracer_index, clock, diffusivities, tracers, buoyancy, velocities)
             + forcing(i, j, k, grid, clock, hydrostatic_prognostic_fields(velocities, free_surface, tracers)))
end

"""
Return the tendency for an explicit free surface at horizontal grid point `i, j`.

The tendency is called ``G_η`` and defined via

    ``∂_t η = G_η``
"""
@inline function free_surface_tendency(i, j, grid,
                                       velocities,
                                       free_surface,
                                       tracers,
                                       auxiliary_fields,
                                       forcings,
                                       clock)

    k_surface = grid.Nz + 1
    model_fields = merge(hydrostatic_prognostic_fields(velocities, free_surface, tracers), auxiliary_fields)

    return @inbounds (   velocities.w[i, j, k_surface]
                       + forcings.η(i, j, k_surface, grid, clock, model_fields))
end

@inline function hydrostatic_turbulent_kinetic_energy_tendency(i, j, k, grid,
                                                               val_tracer_index::Val{tracer_index},
                                                               advection,
                                                               closure,
                                                               buoyancy,
                                                               velocities,
                                                               free_surface,
                                                               tracers,
                                                               diffusivities,
                                                               auxiliary_fields,
                                                               forcing,
                                                               clock) where tracer_index

    tke_index = TKETracerIndex(tracer_index)
    @inbounds e = tracers[tracer_index]

    model_fields = merge(hydrostatic_prognostic_fields(velocities, free_surface, tracers), auxiliary_fields)

    return ( - div_Uc(i, j, k, grid, advection, velocities, e)
             - ∇_dot_qᶜ(i, j, k, grid, closure, e, tke_index, clock, diffusivities, tracers, buoyancy, velocities)
             + shear_production(i, j, k, grid, closure, clock, velocities, tracers, buoyancy, diffusivities)
             + buoyancy_flux(i, j, k, grid, closure, velocities, tracers, buoyancy)
             - dissipation(i, j, k, grid, closure, tracers, buoyancy)
             + forcing(i, j, k, grid, clock, hydrostatic_prognostic_fields(velocities, free_surface, tracers)))
end


text:
#=
TODO
    * move all gl_ methods to GLPlot
    * integrate GLPlot UI
    * clean up corner cases
    * find a cleaner way for extracting properties
    * polar plots
    * labes and axis
    * fix units in all visuals (e.g dotted lines, marker scale, surfaces)
=#

slice_arg(img::Matrix{C}, idx::Int) where {C<:Colorant} = img
is_marker_supported(::GLVisualizeBackend, shape::GLVisualize.AllPrimitives) = true
is_marker_supported(::GLVisualizeBackend, shape::Union{Vector{Matrix{C}}, Matrix{C}}) where {C<:Colorant} = true
is_marker_supported(::GLVisualizeBackend, shape::Shape) = true
GL = Plots

# ---------------------------------------------------------------------------

# initialize the figure/window
# function _create_backend_figure(plt::Plot{GLVisualizeBackend})
#     # init a screen
#
#     GLPlot.init()
# end
const _glplot_deletes = []

function get_plot_screen(list::Vector, name, result = [])
    for elem in list
        get_plot_screen(elem, name, result)
    end
    return result
end
function get_plot_screen(screen, name, result = [])
    if screen.name == name
        push!(result, screen)
        return result
    end
    get_plot_screen(screen.children, name, result)
end

function create_window(plt::Plot{GLVisualizeBackend}, visible)
    name = Symbol("__Plots.jl")
    # make sure we have any screen open
    if isempty(GLVisualize.get_screens())
        # create a fresh, new screen
        parent_screen = GLVisualize.glscreen(
            "Plots",
            resolution = plt[:size],
            visible = visible
        )
        @async GLWindow.renderloop(parent_screen)
        GLVisualize.add_screen(parent_screen)
    end
    # now lets get ourselves a permanent Plotting screen
    plot_screens = get_plot_screen(GLVisualize.current_screen(), name)
    screen = if isempty(plot_screens) # no screen with `name`
        parent = GLVisualize.current_screen()
        screen = GLWindow.Screen(
            parent, area = map(GLWindow.zeroposition, parent.area),
            name = name
        )
        screen
    elseif length(plot_screens) == 1
        plot_screens[1]
    else
        # okay this is silly! Lets see if we can. There is an ID we could use
        # will not be fine for more than 255 screens though -.-.
        error("multiple Plot screens. Please don't use any screen with the name $name")
    end
    # Since we own this window, we can do deep cleansing
    empty!(screen)
    plt.o = screen
    GLWindow.set_visibility!(screen, visible)
    resize!(screen, plt[:size]...)
    screen
end
# ---------------------------------------------------------------------------

const _gl_marker_map = KW(
    :rect => '■',
    :star5 => '★',
    :diamond => '◆',
    :hexagon => '⬢',
    :cross => '✚',
    :xcross => '❌',
    :utriangle => '▲',
    :dtriangle => '▼',
    :ltriangle => '◀',
    :rtriangle => '▶',
    :pentagon => '⬟',
    :octagon => '⯄',
    :star4 => '✦',
    :star6 => '🟋',
    :star8 => '✷',
    :vline => '┃',
    :hline => '━',
    :+ => '+',
    :x => 'x',
    :circle => '●'
)

function gl_marker(shape)
    shape
end
function gl_marker(shape::Shape)
    points = Point2f0[GeometryTypes.Vec{2, Float32}(p) for p in zip(shape.x, shape.y)]
    bb = GeometryTypes.AABB(points)
    mini, maxi = minimum(bb), maximum(bb)
    w3 = maxi-mini
    origin, width = Point2f0(mini[1], mini[2]), Point2f0(w3[1], w3[2])
    map!(p -> ((p - origin) ./ width) - 0.5f0, points, points) # normalize and center
    GeometryTypes.GLNormalMesh(points)
end
# create a marker/shape type
function gl_marker(shape::Vector{Symbol})
    String(map(shape) do sym
        get(_gl_marker_map, sym, '●')
    end)
end

function gl_marker(shape::Symbol)
    if shape == :rect
        GeometryTypes.HyperRectangle(Vec2f0(0), Vec2f0(1))
    elseif shape == :circle || shape == :none
        GeometryTypes.HyperSphere(Point2f0(0), 1f0)
    elseif haskey(_gl_marker_map, shape)
        _gl_marker_map[shape]
    elseif haskey(_shapes, shape)
        gl_marker(_shapes[shape])
    else
        error("Shape $shape not supported by GLVisualize")
    end
end

function extract_limits(sp, plotattributes, kw_args)
    clims = sp[:clims]
    if is_2tuple(clims)
        if isfinite(clims[1]) && isfinite(clims[2])
          kw_args[:limits] = Vec2f0(clims)
        end
    end
    nothing
end

to_vec(::Type{T}, vec::T) where {T <: StaticArrays.StaticVector} = vec
to_vec(::Type{T}, s::Number) where {T <: StaticArrays.StaticVector} = T(s)

to_vec(::Type{T}, vec::StaticArrays.StaticVector{3}) where {T <: StaticArrays.StaticVector{2}} = T(vec[1], vec[2])
to_vec(::Type{T}, vec::StaticArrays.StaticVector{2}) where {T <: StaticArrays.StaticVector{3}} = T(vec[1], vec[2], 0)

to_vec(::Type{T}, vecs::AbstractVector) where {T <: StaticArrays.StaticVector} = map(x-> to_vec(T, x), vecs)

function extract_marker(plotattributes, kw_args)
    dim = Plots.is3d(plotattributes) ? 3 : 2
    scaling = dim == 3 ? 0.003 : 2
    if haskey(plotattributes, :markershape)
        shape = plotattributes[:markershape]
        shape = gl_marker(shape)
        if shape != :none
            kw_args[:primitive] = shape
        end
    end
    dim = isa(kw_args[:primitive], GLVisualize.Sprites) ? 2 : 3
    if haskey(plotattributes, :markersize)
        msize = plotattributes[:markersize]
        kw_args[:scale] = to_vec(GeometryTypes.Vec{dim, Float32}, msize .* scaling)
    end
    if haskey(plotattributes, :offset)
        kw_args[:offset] = plotattributes[:offset]
    end
    # get the color
    key = :markercolor
    haskey(plotattributes, key) || return
    c = gl_color(plotattributes[key])
    if isa(c, AbstractVector) && plotattributes[:marker_z] != nothing
        extract_colornorm(plotattributes, kw_args)
        kw_args[:color] = nothing
        kw_args[:color_map] = c
        kw_args[:intensity] = convert(Vector{Float32}, plotattributes[:marker_z])
    else
        kw_args[:color] = c
    end
    key = :markerstrokecolor
    haskey(plotattributes, key) || return
    c = gl_color(plotattributes[key])
    if c != nothing
        if !(isa(c, Colorant) || (isa(c, Vector) && eltype(c) <: Colorant))
            error("Stroke Color not supported: $c")
        end
        kw_args[:stroke_color] = c
        kw_args[:stroke_width] = Float32(plotattributes[:markerstrokewidth])
    end
end

function _extract_surface(plotattributes::Plots.Surface)
    plotattributes.surf
end
function _extract_surface(plotattributes::AbstractArray)
    plotattributes
end

# TODO when to transpose??
function extract_surface(plotattributes)
    map(_extract_surface, (plotattributes[:x], plotattributes[:y], plotattributes[:z]))
end
function topoints(::Type{P}, array) where P
    [P(x) for x in zip(array...)]
end
function extract_points(plotattributes)
    dim = is3d(plotattributes) ? 3 : 2
    array = if plotattributes[:seriestype] == :straightline
        straightline_data(plotattributes)
    elseif plotattributes[:seriestype] == :shape
        shape_data(plotattributes)
    else
        (plotattributes[:x], plotattributes[:y], plotattributes[:z])[1:dim]
    end
    topoints(Point{dim, Float32}, array)
end
function make_gradient(grad::Vector{C}) where C <: Colorant
    grad
end
function make_gradient(grad::ColorGradient)
    RGBA{Float32}[c for c in grad.colors]
end
make_gradient(c) = make_gradient(cgrad())

function extract_any_color(plotattributes, kw_args)
    if plotattributes[:marker_z] == nothing
        c = scalar_color(plotattributes, :fill)
        extract_c(plotattributes, kw_args, :fill)
        if isa(c, Colorant)
            kw_args[:color] = c
        else
            kw_args[:color] = nothing
            kw_args[:color_map] = make_gradient(c)
            clims = plotattributes[:subplot][:clims]
            if Plots.is_2tuple(clims)
                if isfinite(clims[1]) && isfinite(clims[2])
                    kw_args[:color_norm] = Vec2f0(clims)
                end
            elseif clims == :auto
                kw_args[:color_norm] = Vec2f0(ignorenan_extrema(plotattributes[:y]))
            end
        end
    else
        kw_args[:color] = nothing
        clims = plotattributes[:subplot][:clims]
        if Plots.is_2tuple(clims)
            if isfinite(clims[1]) && isfinite(clims[2])
                kw_args[:color_norm] = Vec2f0(clims)
            end
        elseif clims == :auto
            kw_args[:color_norm] = Vec2f0(ignorenan_extrema(plotattributes[:y]))
        else
            error("Unsupported limits: $clims")
        end
        kw_args[:intensity] = convert(Vector{Float32}, plotattributes[:marker_z])
        kw_args[:color_map] = gl_color_map(plotattributes, :marker)
    end
end

function extract_stroke(plotattributes, kw_args)
    extract_c(plotattributes, kw_args, :line)
    if haskey(plotattributes, :linewidth)
        kw_args[:thickness] = Float32(plotattributes[:linewidth] * 3)
    end
end

function extract_color(plotattributes, sym)
    plotattributes[Symbol("$(sym)color")]
end

gl_color(c::PlotUtils.ColorGradient) = c.colors
gl_color(c::Vector{T}) where {T<:Colorant} = c
gl_color(c::RGBA{Float32}) = c
gl_color(c::Colorant) = RGBA{Float32}(c)

function gl_color(tuple::Tuple)
    gl_color(tuple...)
end

# convert to RGBA
function gl_color(c, a)
    c = convertColor(c, a)
    RGBA{Float32}(c)
end
function scalar_color(plotattributes, sym)
    gl_color(extract_color(plotattributes, sym))
end

function gl_color_map(plotattributes, sym)
    colors = extract_color(plotattributes, sym)
    _gl_color_map(colors)
end
function _gl_color_map(colors::PlotUtils.ColorGradient)
    colors.colors
end
function _gl_color_map(c)
    Plots.default_gradient()
end

dist(a, b) = abs(a-b)
mindist(x, a, b) = NaNMath.min(dist(a, x), dist(b, x))

function gappy(x, ps)
    n = length(ps)
    x <= first(ps) && return first(ps) - x
    for j=1:(n-1)
        p0 = ps[j]
        p1 = ps[NaNMath.min(j+1, n)]
        if p0 <= x && p1 >= x
            return mindist(x, p0, p1) * (isodd(j) ? 1 : -1)
        end
    end
    return last(ps) - x
end
function ticks(points, resolution)
    Float16[gappy(x, points) for x = range(first(points), stop=last(points), length=resolution)]
end

function insert_pattern!(points, kw_args)
    tex = GLAbstraction.Texture(ticks(points, 100), x_repeat=:repeat)
    kw_args[:pattern] = tex
    kw_args[:pattern_length] = Float32(last(points))
end
function extract_linestyle(plotattributes, kw_args)
    haskey(plotattributes, :linestyle) || return
    ls = plotattributes[:linestyle]
    lw = plotattributes[:linewidth]
    kw_args[:thickness] = Float32(lw)
    if ls == :dash
        points = [0.0, lw, 2lw, 3lw, 4lw]
        insert_pattern!(points, kw_args)
    elseif ls == :dot
        tick, gap = lw/2, lw/4
        points = [0.0, tick, tick+gap, 2tick+gap, 2tick+2gap]
        insert_pattern!(points, kw_args)
    elseif ls == :dashdot
        dtick, dgap = lw, lw
        ptick, pgap = lw/2, lw/4
        points = [0.0, dtick, dtick+dgap, dtick+dgap+ptick, dtick+dgap+ptick+pgap]
        insert_pattern!(points, kw_args)
    elseif ls == :dashdotdot
        dtick, dgap = lw, lw
        ptick, pgap = lw/2, lw/4
        points = [0.0, dtick, dtick+dgap, dtick+dgap+ptick, dtick+dgap+ptick+pgap, dtick+dgap+ptick+pgap+ptick,  dtick+dgap+ptick+pgap+ptick+pgap]
        insert_pattern!(points, kw_args)
    end
    extract_c(plotattributes, kw_args, :line)
    nothing
end

function hover(to_hover::Vector, to_display, window)
    hover(to_hover[], to_display, window)
end

function get_cam(x)
    if isa(x, GLAbstraction.Context)
        return get_cam(x.children)
    elseif isa(x, Vector)
        return get_cam(first(x))
    elseif isa(x, GLAbstraction.RenderObject)
        return x[:preferred_camera]
    end
end

function hover(to_hover, to_display, window)
    if isa(to_hover, GLAbstraction.Context)
        return hover(to_hover.children, to_display, window)
    end
    area = map(window.inputs[:mouseposition]) do mp
        SimpleRectangle{Int}(round(Int, mp+10)..., 100, 70)
    end
    mh = GLWindow.mouse2id(window)
    popup = GLWindow.Screen(
        window,
        hidden = map(mh-> !(mh.id == to_hover.id), mh),
        area = area,
        stroke = (2f0, RGBA(0f0, 0f0, 0f0, 0.8f0))
    )
    cam = get!(popup.cameras, :perspective) do
        GLAbstraction.PerspectiveCamera(
            popup.inputs, Vec3f0(3), Vec3f0(0),
            keep = Signal(false),
            theta = Signal(Vec3f0(0)), trans = Signal(Vec3f0(0))
        )
    end

    map(enumerate(to_display)) do id
        i,d = id
        robj = visualize(d)
        viewit = Reactive.droprepeats(map(mh->mh.id == to_hover.id && mh.index == i, mh))
        camtype = get_cam(robj)
        Reactive.preserve(map(viewit) do vi
            if vi
                empty!(popup)
                if camtype == :perspective
                    cam.projectiontype.value = GLVisualize.PERSPECTIVE
                else
                    cam.projectiontype.value = GLVisualize.ORTHOGRAPHIC
                end
                GLVisualize._view(robj, popup, camera = cam)
                bb = GLAbstraction.boundingbox(robj).value
                mini = minimum(bb)
                w = GeometryTypes.widths(bb)
                wborder = w * 0.08f0 #8 percent border
                bb = GeometryTypes.AABB{Float32}(mini - wborder, w + 2 * wborder)
                GLAbstraction.center!(cam, bb)
            end
        end)
    end
    nothing
end

function extract_extrema(plotattributes, kw_args)
    xmin, xmax = ignorenan_extrema(plotattributes[:x]); ymin, ymax = ignorenan_extrema(plotattributes[:y])
    kw_args[:primitive] = GeometryTypes.SimpleRectangle{Float32}(xmin, ymin, xmax-xmin, ymax-ymin)
    nothing
end

function extract_font(font, kw_args)
    kw_args[:family] = font.family
    kw_args[:relative_scale] = pointsize(font)
    kw_args[:color] = gl_color(font.color)
end

function extract_colornorm(plotattributes, kw_args)
    clims = plotattributes[:subplot][:clims]
    if Plots.is_2tuple(clims)
        if isfinite(clims[1]) && isfinite(clims[2])
            kw_args[:color_norm] = Vec2f0(clims)
        end
    elseif clims == :auto
        z = if haskey(plotattributes, :marker_z) && plotattributes[:marker_z] != nothing
            plotattributes[:marker_z]
        elseif haskey(plotattributes, :line_z) && plotattributes[:line_z] != nothing
            plotattributes[:line_z]
        elseif isa(plotattributes[:z], Plots.Surface)
            plotattributes[:z].surf
        else
            plotattributes[:y]
        end
        kw_args[:color_norm] = Vec2f0(ignorenan_extrema(z))
        kw_args[:intensity] = map(Float32, collect(z))
    end
end

function extract_gradient(plotattributes, kw_args, sym)
    key = Symbol("$(sym)color")
    haskey(plotattributes, key) || return
    c = make_gradient(plotattributes[key])
    kw_args[:color] = nothing
    extract_colornorm(plotattributes, kw_args)
    kw_args[:color_map] = c
    return
end

function extract_c(plotattributes, kw_args, sym)
    key = Symbol("$(sym)color")
    haskey(plotattributes, key) || return
    c = gl_color(plotattributes[key])
    kw_args[:color] = nothing
    kw_args[:color_map] = nothing
    kw_args[:color_norm] = nothing
    if (
            isa(c, AbstractVector) &&
            ((haskey(plotattributes, :marker_z) && plotattributes[:marker_z] != nothing) ||
            (haskey(plotattributes, :line_z) && plotattributes[:line_z] != nothing))
        )
        extract_colornorm(plotattributes, kw_args)
        kw_args[:color_map] = c
    else
        kw_args[:color] = c
    end
    return
end

function extract_stroke(plotattributes, kw_args, sym)
    key = Symbol("$(sym)strokecolor")
    haskey(plotattributes, key) || return
    c = gl_color(plotattributes[key])
    if c != nothing
        if !isa(c, Colorant)
            error("Stroke Color not supported: $c")
        end
        kw_args[:stroke_color] = c
        kw_args[:stroke_width] = Float32(plotattributes[Symbol("$(sym)strokewidth")]) * 2
    end
    return
end

function draw_grid_lines(sp, grid_segs, thickness, style, model, color)

    kw_args = Dict{Symbol, Any}(
        :model => model
    )
    plotattributes = Dict(
        :linestyle => style,
        :linewidth => Float32(thickness),
        :linecolor => color
    )
    Plots.extract_linestyle(plotattributes, kw_args)
    GL.gl_lines(map(Point2f0, grid_segs.pts), kw_args)
end

function align_offset(startpos, lastpos, atlas, rscale, font, align)
    xscale, yscale = GLVisualize.glyph_scale!('X', rscale)
    xmove = (lastpos-startpos)[1] + xscale
    if isa(align, GeometryTypes.Vec)
        return -Vec2f0(xmove, yscale) .* align
    elseif align == :top
        return -Vec2f0(xmove/2f0, yscale)
    elseif align == :right
        return -Vec2f0(xmove, yscale/2f0)
    else
        error("Align $align not known")
    end
end

function alignment2num(x::Symbol)
    (x in (:hcenter, :vcenter)) && return 0.5
    (x in (:left, :bottom)) && return 0.0
    (x in (:right, :top)) && return 1.0
    0.0 # 0 default, or better to error?
end

function alignment2num(font::Plots.Font)
    Vec2f0(map(alignment2num, (font.halign, font.valign)))
end

pointsize(font) = font.pointsize * 2

function draw_ticks(
        axis, ticks, isx, isorigin, lims, m, text = "",
        positions = Point2f0[], offsets=Vec2f0[]
    )
    sz = pointsize(tickfont(axis))
    atlas = GLVisualize.get_texture_atlas()
    font = GLVisualize.defaultfont()

    flip = axis[:flip]; mirror = axis[:mirror]

    align = if isx
        mirror ? :bottom : :top
    else
        mirror ? :left : :right
    end
    axis_gap = Point2f0(isx ? 0 : sz / 2, isx ? sz / 2 : 0)
    for (cv, dv) in zip(ticks...)

        x, y = cv, lims[1]
        xy = if isorigin
            isx ? (x, 0) : (0, x)
        else
            isx ? (x, y) : (y, x)
        end
        _pos = m * GeometryTypes.Vec4f0(xy[1], xy[2], 0, 1)
        startpos = Point2f0(_pos[1], _pos[2]) - axis_gap
        str = string(dv)
        # need to tag a new UnicodeFun version for this... also the numbers become
        # so small that it looks terrible -.-
        # _str = split(string(dv), "^")
        # if length(_str) == 2
        #     _str[2] = UnicodeFun.to_superscript(_str[2])
        # end
        # str = join(_str, "")
        position = GLVisualize.calc_position(str, startpos, sz, font, atlas)
        offset = GLVisualize.calc_offset(str, sz, font, atlas)
        alignoff = align_offset(startpos, last(position), atlas, sz, font, align)
        map!(position, position) do pos
            pos .+ alignoff
        end
        append!(positions, position)
        append!(offsets, offset)
        text *= str

    end
    text, positions, offsets
end

function glvisualize_text(position, text, kw_args)
    text_align = alignment2num(text.font)
    startpos = Vec2f0(position)
    atlas = GLVisualize.get_texture_atlas()
    font = GLVisualize.defaultfont()
    rscale = kw_args[:relative_scale]

    position = GLVisualize.calc_position(text.str, startpos, rscale, font, atlas)
    offset = GLVisualize.calc_offset(text.str, rscale, font, atlas)
    alignoff = align_offset(startpos, last(position), atlas, rscale, font, text_align)

    map!(position, position) do pos
        pos .+ alignoff
    end
    kw_args[:position] = position
    kw_args[:offset] = offset
    kw_args[:scale_primitive] = true
    visualize(text.str, Style(:default), kw_args)
end

function text_model(font, pivot)
    pv = GeometryTypes.Vec3f0(pivot[1], pivot[2], 0)
    if font.rotation != 0.0
        rot = Float32(deg2rad(font.rotation))
        rotm = GLAbstraction.rotationmatrix_z(rot)
        return GLAbstraction.translationmatrix(pv)*rotm*GLAbstraction.translationmatrix(-pv)
    else
        eye(GeometryTypes.Mat4f0)
    end
end
function gl_draw_axes_2d(sp::Plots.Subplot{Plots.GLVisualizeBackend}, model, area)
    xticks, yticks, xspine_segs, yspine_segs, xtick_segs, ytick_segs, xgrid_segs, ygrid_segs, xminorgrid_segs, yminorgrid_segs, xborder_segs, yborder_segs = Plots.axis_drawing_info(sp)
    xaxis = sp[:xaxis]; yaxis = sp[:yaxis]

    xgc = Colors.color(Plots.gl_color(xaxis[:foreground_color_grid]))
    ygc = Colors.color(Plots.gl_color(yaxis[:foreground_color_grid]))
    axis_vis = []
    if xaxis[:grid]
        grid = draw_grid_lines(sp, xgrid_segs, xaxis[:gridlinewidth], xaxis[:gridstyle], model, RGBA(xgc, xaxis[:gridalpha]))
        push!(axis_vis, grid)
    end
    if yaxis[:grid]
        grid = draw_grid_lines(sp, ygrid_segs, yaxis[:gridlinewidth], yaxis[:gridstyle], model, RGBA(ygc, yaxis[:gridalpha]))
        push!(axis_vis, grid)
    end
    if xaxis[:minorgrid]
        minorgrid = draw_minorgrid_lines(sp, xminorgrid_segs, xaxis[:minorgridlinewidth], xaxis[:minorgridstyle], model, RGBA(xgc, xaxis[:minorgridalpha]))
        push!(axis_vis, minorgrid)
    end
    if yaxis[:minorgrid]
        minorgrid = draw_minorgrid_lines(sp, yminorgrid_segs, yaxis[:minorgridlinewidth], yaxis[:minorgridstyle], model, RGBA(ygc, yaxis[:minorgridalpha]))
        push!(axis_vis, minorgrid)
    end

    xac = Colors.color(Plots.gl_color(xaxis[:foreground_color_axis]))
    yac = Colors.color(Plots.gl_color(yaxis[:foreground_color_axis]))
    if alpha(xaxis[:foreground_color_axis]) > 0
        spine = draw_grid_lines(sp, xspine_segs, 1f0, :solid, model, RGBA(xac, 1.0f0))
        push!(axis_vis, spine)
    end
    if alpha(yaxis[:foreground_color_axis]) > 0
        spine = draw_grid_lines(sp, yspine_segs, 1f0, :solid, model, RGBA(yac, 1.0f0))
        push!(axis_vis, spine)
    end
    if sp[:framestyle] in (:zerolines, :grid)
        if alpha(xaxis[:foreground_color_grid]) > 0
            spine = draw_grid_lines(sp, xtick_segs, 1f0, :solid, model, RGBA(xgc, xaxis[:gridalpha]))
            push!(axis_vis, spine)
        end
        if alpha(yaxis[:foreground_color_grid]) > 0
            spine = draw_grid_lines(sp, ytick_segs, 1f0, :solid, model, RGBA(ygc, yaxis[:gridalpha]))
            push!(axis_vis, spine)
        end
    else
        if alpha(xaxis[:foreground_color_axis]) > 0
            spine = draw_grid_lines(sp, xtick_segs, 1f0, :solid, model, RGBA(xac, 1.0f0))
            push!(axis_vis, spine)
        end
        if alpha(yaxis[:foreground_color_axis]) > 0
            spine = draw_grid_lines(sp, ytick_segs, 1f0, :solid, model, RGBA(yac, 1.0f0))
            push!(axis_vis, spine)
        end
    end
    fcolor = Plots.gl_color(xaxis[:foreground_color_axis])

    xlim = Plots.axis_limits(xaxis)
    ylim = Plots.axis_limits(yaxis)

    if !(xaxis[:ticks] in (nothing, false, :none)) && !(sp[:framestyle] == :none) && xaxis[:showaxis]
        ticklabels = map(model) do m
            mirror = xaxis[:mirror]
            t, positions, offsets = draw_ticks(xaxis, xticks, true, sp[:framestyle] == :origin, ylim, m)
        end
        kw_args = Dict{Symbol, Any}(
            :position => map(x-> x[2], ticklabels),
            :offset => map(last, ticklabels),
            :color => fcolor,
            :relative_scale => pointsize(tickfont(xaxis)),
            :scale_primitive => false
        )
        push!(axis_vis, visualize(map(first, ticklabels), Style(:default), kw_args))
    end

    if !(yaxis[:ticks] in (nothing, false, :none)) && !(sp[:framestyle] == :none) && yaxis[:showaxis]
        ticklabels = map(model) do m
            mirror = yaxis[:mirror]
            t, positions, offsets = draw_ticks(yaxis, yticks, false, sp[:framestyle] == :origin, xlim, m)
        end
        kw_args = Dict{Symbol, Any}(
            :position => map(x-> x[2], ticklabels),
            :offset => map(last, ticklabels),
            :color => fcolor,
            :relative_scale => pointsize(tickfont(xaxis)),
            :scale_primitive => false
        )
        push!(axis_vis, visualize(map(first, ticklabels), Style(:default), kw_args))
    end

    xbc = Colors.color(Plots.gl_color(xaxis[:foreground_color_border]))
    ybc = Colors.color(Plots.gl_color(yaxis[:foreground_color_border]))
    intensity = sp[:framestyle] == :semi ? 0.5f0 : 1.0f0
    if sp[:framestyle] in (:box, :semi)
        xborder = draw_grid_lines(sp, xborder_segs, intensity, :solid, model, RGBA(xbc, intensity))
        yborder = draw_grid_lines(sp, yborder_segs, intensity, :solid, model, RGBA(ybc, intensity))
        push!(axis_vis, xborder, yborder)
    end

    area_w = GeometryTypes.widths(area)
    if sp[:title] != ""
        tf = titlefont(sp)
        font = Plots.Font(tf.family, tf.pointsize, :hcenter, :top, tf.rotation, tf.color)
        xy = Point2f0(area.w/2, area_w[2] + pointsize(tf)/2)
        kw = Dict(:model => text_model(font, xy), :scale_primitive => true)
        extract_font(font, kw)
        t = PlotText(sp[:title], font)
        push!(axis_vis, glvisualize_text(xy, t, kw))
    end
    if xaxis[:guide] != ""
        tf = guidefont(xaxis)
        xy = Point2f0(area.w/2, - pointsize(tf)/2)
        font = Plots.Font(tf.family, tf.pointsize, :hcenter, :bottom, tf.rotation, tf.color)
        kw = Dict(:model => text_model(font, xy), :scale_primitive => true)
        t = PlotText(xaxis[:guide], font)
        extract_font(font, kw)
        push!(axis_vis, glvisualize_text(xy, t, kw))
    end

    if yaxis[:guide] != ""
        tf = guidefont(yaxis)
        font = Plots.Font(tf.family, tf.pointsize, :hcenter, :top, 90f0, tf.color)
        xy = Point2f0(-pointsize(tf)/2, area.h/2)
        kw = Dict(:model => text_model(font, xy), :scale_primitive=>true)
        t = PlotText(yaxis[:guide], font)
        extract_font(font, kw)
        push!(axis_vis, glvisualize_text(xy, t, kw))
    end

    axis_vis
end

function gl_draw_axes_3d(sp, model)
    x = Plots.axis_limits(sp[:xaxis])
    y = Plots.axis_limits(sp[:yaxis])
    z = Plots.axis_limits(sp[:zaxis])

    min = Vec3f0(x[1], y[1], z[1])
    visualize(
        GeometryTypes.AABB{Float32}(min, Vec3f0(x[2], y[2], z[2])-min),
        :grid, model=model
    )
end

function gl_bar(plotattributes, kw_args)
    x, y = plotattributes[:x], plotattributes[:y]
    nx, ny = length(x), length(y)
    axis = plotattributes[:subplot][isvertical(plotattributes) ? :xaxis : :yaxis]
    cv = [discrete_value!(axis, xi)[1] for xi=x]
    x = if nx == ny
        cv
    elseif nx == ny + 1
        0.5diff(cv) + cv[1:end-1]
    else
        error("bar recipe: x must be same length as y (centers), or one more than y (edges).\n\t\tlength(x)=$(length(x)), length(y)=$(length(y))")
    end
    if haskey(kw_args, :stroke_width) # stroke is inside for bars
        #kw_args[:stroke_width] = -kw_args[:stroke_width]
    end
    # compute half-width of bars
    bw = nothing
    hw = if bw == nothing
        ignorenan_mean(diff(x))
    else
        Float64[_cycle(bw,i)*0.5 for i=1:length(x)]
    end

    # make fillto a vector... default fills to 0
    fillto = plotattributes[:fillrange]
    if fillto == nothing
        fillto = 0
    end
    # create the bar shapes by adding x/y segments
    positions, scales = Array{Point2f0}(undef, ny), Array{Vec2f0}(undef, ny)
    m = Reactive.value(kw_args[:model])
    sx, sy = m[1,1], m[2,2]
    for i=1:ny
        center = x[i]
        hwi = abs(_cycle(hw,i)); yi = y[i]; fi = _cycle(fillto,i)
        if Plots.isvertical(plotattributes)
            sz = (hwi*sx, yi*sy)
        else
            sz = (yi*sx, hwi*2*sy)
        end
        positions[i] = (center-hwi*0.5, fi)
        scales[i] = sz
    end

    kw_args[:scale] = scales
    kw_args[:offset] = Vec2f0(0)
    visualize((GLVisualize.RECTANGLE, positions), Style(:default), kw_args)
    #[]
end

const _box_halfwidth = 0.4

notch_width(q2, q4, N) = 1.58 * (q4-q2)/sqrt(N)

function gl_boxplot(plotattributes, kw_args)
    kwbox = copy(kw_args)
    range = 1.5; notch = false
    x, y = plotattributes[:x], plotattributes[:y]
    glabels = sort(collect(unique(x)))
    warning = false
    outliers_x, outliers_y = zeros(0), zeros(0)

    box_pos = Point2f0[]
    box_scale = Vec2f0[]
    outliers = Point2f0[]
    t_segments = Point2f0[]
    m = Reactive.value(kw_args[:model])
    sx, sy = m[1,1], m[2,2]
    for (i,glabel) in enumerate(glabels)
        # filter y
        values = y[filter(i -> _cycle(x,i) == glabel, 1:length(y))]
        # compute quantiles
        q1,q2,q3,q4,q5 = quantile(values, range(0, stop=1, length=5))
        # notch
        n = Plots.notch_width(q2, q4, length(values))
        # warn on inverted notches?
        if notch && !warning && ( (q2>(q3-n)) || (q4<(q3+n)) )
            @warn("Boxplot's notch went outside hinges. Set notch to false.")
            warning = true # Show the warning only one time
        end

        # make the shape
        center = Plots.discrete_value!(plotattributes[:subplot][:xaxis], glabel)[1]
        hw = plotattributes[:bar_width] == nothing ? Plots._box_halfwidth*2 : _cycle(plotattributes[:bar_width], i)
        l, m, r = center - hw/2, center, center + hw/2

        # internal nodes for notches
        L, R = center - 0.5 * hw, center + 0.5 * hw
        # outliers
        if Float64(range) != 0.0  # if the range is 0.0, the whiskers will extend to the data
            limit = range*(q4-q2)
            inside = Float64[]
            for value in values
                if (value < (q2 - limit)) || (value > (q4 + limit))
                    push!(outliers, (center, value))
                else
                    push!(inside, value)
                end
            end
            # change q1 and q5 to show outliers
            # using maximum and minimum values inside the limits
            q1, q5 = ignorenan_extrema(inside)
        end
        # Box
        if notch
            push!(t_segments, (m, q1), (l, q1), (r, q1), (m, q1), (m, q2))# lower T
            push!(box_pos, (l, q2));push!(box_scale, (hw*sx, n*sy)) # lower box
            push!(box_pos, (l, q4));push!(box_scale, (hw*sx, n*sy)) # upper box
            push!(t_segments, (m, q5), (l, q5), (r, q5), (m, q5), (m, q4))# upper T

        else
            push!(t_segments, (m, q2), (m, q1), (l, q1), (r, q1))# lower T
            push!(box_pos, (l, q2)); push!(box_scale, (hw*sx, (q3-q2)*sy)) # lower box
            push!(box_pos, (l, q4)); push!(box_scale, (hw*sx, (q3-q4)*sy)) # upper box
            push!(t_segments, (m, q4), (m, q5), (r, q5), (l, q5))# upper T
        end
    end
    kwbox = Dict{Symbol, Any}(
        :scale => box_scale,
        :model => kw_args[:model],
        :offset => Vec2f0(0),
    )
    extract_marker(plotattributes, kw_args)
    outlier_kw = Dict(
        :model => kw_args[:model],
        :color =>  scalar_color(plotattributes, :fill),
        :stroke_width => Float32(plotattributes[:markerstrokewidth]),
        :stroke_color => scalar_color(plotattributes, :markerstroke),
    )
    lines_kw = Dict(
        :model => kw_args[:model],
        :stroke_width =>  plotattributes[:linewidth],
        :stroke_color =>  scalar_color(plotattributes, :fill),
    )
    vis1 = GLVisualize.visualize((GLVisualize.RECTANGLE, box_pos), Style(:default), kwbox)
    vis2 = GLVisualize.visualize((GLVisualize.CIRCLE, outliers), Style(:default), outlier_kw)
    vis3 = GLVisualize.visualize(t_segments, Style(:linesegment), lines_kw)
    [vis1, vis2, vis3]
end

# ---------------------------------------------------------------------------
function gl_viewport(bb, rect)
    l, b, bw, bh = bb
    rw, rh = rect.w, rect.h
    GLVisualize.SimpleRectangle(
        round(Int, rw * l),
        round(Int, rh * b),
        round(Int, rw * bw),
        round(Int, rh * bh)
    )
end

function to_modelmatrix(rect, subrect, rel_plotarea, sp)
    xmin, xmax = Plots.axis_limits(sp[:xaxis])
    ymin, ymax = Plots.axis_limits(sp[:yaxis])
    mini, maxi = Vec3f0(xmin, ymin, 0), Vec3f0(xmax, ymax, 1)
    if Plots.is3d(sp)
        zmin, zmax = Plots.axis_limits(sp[:zaxis])
        mini, maxi = Vec3f0(xmin, ymin, zmin), Vec3f0(xmax, ymax, zmax)
        s = Vec3f0(1) ./ (maxi-mini)
        return GLAbstraction.scalematrix(s)*GLAbstraction.translationmatrix(-mini)
    end
    l, b, bw, bh = rel_plotarea
    w, h = rect.w*bw, rect.h*bh
    x, y = rect.w*l - subrect.x, rect.h*b - subrect.y
    t = -mini
    s = Vec3f0(w, h, 1) ./ (maxi-mini)
    GLAbstraction.translationmatrix(Vec3f0(x,y,0))*GLAbstraction.scalematrix(s)*GLAbstraction.translationmatrix(t)
end

# ----------------------------------------------------------------

function scale_for_annotations!(series::Series, scaletype::Symbol = :pixels)
    anns = series[:series_annotations]
    if anns != nothing && anns.baseshape != nothing
        # we use baseshape to overwrite the markershape attribute
        # with a list of custom shapes for each
        msw, msh = anns.scalefactor
        offsets = Array{Vec2f0}(undef, length(anns.strs))
        series[:markersize] = map(1:length(anns.strs)) do i
            str = _cycle(anns.strs, i)
            # get the width and height of the string (in mm)
            sw, sh = text_size(str, anns.font.pointsize)

            # how much to scale the base shape?
            # note: it's a rough assumption that the shape fills the unit box [-1,-1,1,1],
            # so we scale the length-2 shape by 1/2 the total length
            xscale = 0.5to_pixels(sw) * 1.8
            yscale = 0.5to_pixels(sh) * 1.8

            # we save the size of the larger direction to the markersize list,
            # and then re-scale a copy of baseshape to match the w/h ratio
            s = Vec2f0(xscale, yscale)
            offsets[i] = -s
            s
        end
        series[:offset] = offsets
    end
    return
end

function _display(plt::Plot{GLVisualizeBackend}, visible = true)
    screen = create_window(plt, visible)
    sw, sh = plt[:size]
    sw, sh = sw*px, sh*px

    for sp in plt.subplots
        _3d = Plots.is3d(sp)
        # camera = :perspective
        # initialize the sub-screen for this subplot
        rel_bbox = Plots.bbox_to_pcts(bbox(sp), sw, sh)
        sub_area = map(screen.area) do rect
            Plots.gl_viewport(rel_bbox, rect)
        end
        c = plt[:background_color_outside]
        sp_screen = GLVisualize.Screen(
            screen, color = c,
            area = sub_area
        )
        sp.o = sp_screen
        cam = get!(sp_screen.cameras, :perspective) do
            inside = sp_screen.inputs[:mouseinside]
            theta = _3d ? nothing : Signal(Vec3f0(0)) # surpress rotation for 2D (nothing will get usual rotation controle)
            GLAbstraction.PerspectiveCamera(
                sp_screen.inputs, Vec3f0(3), Vec3f0(0),
                keep = inside, theta = theta
            )
        end

        rel_plotarea = Plots.bbox_to_pcts(plotarea(sp), sw, sh)
        model_m = map(Plots.to_modelmatrix,
            screen.area, sub_area,
            Signal(rel_plotarea), Signal(sp)
        )

        # loop over the series and add them to the subplot
        if !_3d
            axis = gl_draw_axes_2d(sp, model_m, Reactive.value(sub_area))
            GLVisualize._view(axis, sp_screen, camera=:perspective)
            cam.projectiontype.value = GLVisualize.ORTHOGRAPHIC
            Reactive.run_till_now() # make sure Reactive.push! arrives
            GLAbstraction.center!(cam,
                GeometryTypes.AABB(
                    Vec3f0(-20), Vec3f0((GeometryTypes.widths(sp_screen)+40f0)..., 1)
                )
            )
        else
            axis = gl_draw_axes_3d(sp, model_m)
            GLVisualize._view(axis, sp_screen, camera=:perspective)
            push!(cam.projectiontype, GLVisualize.PERSPECTIVE)
        end
        for series in Plots.series_list(sp)

            plotattributes = series.plotattributes
            st = plotattributes[:seriestype]; kw_args = KW() # exctract kw

            kw_args[:model] = model_m # add transformation
            if !_3d # 3D is treated differently, since we need boundingboxes for camera
                kw_args[:boundingbox] = nothing # don't calculate bb, we dont need it
            end
            scale_for_annotations!(series)
            if st in (:surface, :wireframe)
                x, y, z = extract_surface(plotattributes)
                extract_gradient(plotattributes, kw_args, :fill)
                z = Plots.transpose_z(plotattributes, z, false)
                if isa(x, AbstractMatrix) && isa(y, AbstractMatrix)
                    x, y = Plots.transpose_z(plotattributes, x, false), Plots.transpose_z(plotattributes, y, false)
                end
                if st == :wireframe
                    kw_args[:wireframe] = true
                    kw_args[:stroke_color] = plotattributes[:linecolor]
                    kw_args[:stroke_width] = Float32(plotattributes[:linewidth]/100f0)
                end
                vis = GL.gl_surface(x, y, z, kw_args)
            elseif (st in (:path, :path3d, :straightline)) && plotattributes[:linewidth] > 0
                kw = copy(kw_args)
                points = Plots.extract_points(plotattributes)
                extract_linestyle(plotattributes, kw)
                vis = GL.gl_lines(points, kw)
                if plotattributes[:markershape] != :none
                    kw = copy(kw_args)
                    extract_stroke(plotattributes, kw)
                    extract_marker(plotattributes, kw)
                    vis2 = GL.gl_scatter(copy(points), kw)
                    vis = [vis; vis2]
                end
                if plotattributes[:fillrange] != nothing
                    kw = copy(kw_args)
                    fr = plotattributes[:fillrange]
                    ps = if all(x-> x >= 0, diff(plotattributes[:x])) # if is monotonic
                        vcat(points, Point2f0[(points[i][1], _cycle(fr, i)) for i=length(points):-1:1])
                    else
                        points
                    end
                    extract_c(plotattributes, kw, :fill)
                    vis = [GL.gl_poly(ps, kw), vis]
                end
            elseif st in (:scatter, :scatter3d) #|| plotattributes[:markershape] != :none
                extract_marker(plotattributes, kw_args)
                points = extract_points(plotattributes)
                vis = GL.gl_scatter(points, kw_args)
            elseif st == :shape
                extract_c(plotattributes, kw_args, :fill)
                vis = GL.gl_shape(plotattributes, kw_args)
            elseif st == :contour
                x,y,z = extract_surface(plotattributes)
                z = transpose_z(plotattributes, z, false)
                extract_extrema(plotattributes, kw_args)
                extract_gradient(plotattributes, kw_args, :fill)
                kw_args[:fillrange] = plotattributes[:fillrange]
                kw_args[:levels] = plotattributes[:levels]

                vis = GL.gl_contour(x,y,z, kw_args)
            elseif st == :heatmap
                x,y,z = extract_surface(plotattributes)
                extract_gradient(plotattributes, kw_args, :fill)
                extract_extrema(plotattributes, kw_args)
                extract_limits(sp, plotattributes, kw_args)
                vis = GL.gl_heatmap(x,y,z, kw_args)
            elseif st == :bar
                extract_c(plotattributes, kw_args, :fill)
                extract_stroke(plotattributes, kw_args, :marker)
                vis = gl_bar(plotattributes, kw_args)
            elseif st == :image
                extract_extrema(plotattributes, kw_args)
                vis = GL.gl_image(plotattributes[:z].surf, kw_args)
            elseif st == :boxplot
                 extract_c(plotattributes, kw_args, :fill)
                 vis = gl_boxplot(plotattributes, kw_args)
             elseif st == :volume
                  volume = plotattributes[:y]
                  _plotattributes = copy(plotattributes)
                  _plotattributes[:y] = 0:1
                  _plotattributes[:x] = 0:1
                  kw_args = KW()
                  extract_gradient(_plotattributes, kw_args, :fill)
                  vis = visualize(volume.v, Style(:default), kw_args)
             else
                error("failed to display plot type $st")
            end

            isa(vis, Array) && isempty(vis) && continue # nothing to see here

            GLVisualize._view(vis, sp_screen, camera=:perspective)
            if haskey(plotattributes, :hover) && !(plotattributes[:hover] in (false, :none, nothing))
                hover(vis, plotattributes[:hover], sp_screen)
            end
            if isdefined(:GLPlot) && isdefined(Main.GLPlot, :(register_plot!))
                del_signal = Main.GLPlot.register_plot!(vis, sp_screen, create_gizmo=false)
                append!(_glplot_deletes, del_signal)
            end
            anns = series[:series_annotations]
            for (x, y, str, font) in EachAnn(anns, plotattributes[:x], plotattributes[:y])
                txt_args = Dict{Symbol, Any}(:model => eye(GLAbstraction.Mat4f0))
                x, y = Reactive.value(model_m) * GeometryTypes.Vec{4, Float32}(x, y, 0, 1)
                extract_font(font, txt_args)
                t = glvisualize_text(Point2f0(x, y), PlotText(str, font), txt_args)
                GLVisualize._view(t, sp_screen, camera = :perspective)
            end

        end
        generate_legend(sp, sp_screen, model_m)
        if _3d
            GLAbstraction.center!(sp_screen)
        end
        GLAbstraction.post_empty()
        yield()
    end
end

function _show(io::IO, ::MIME"image/png", plt::Plot{GLVisualizeBackend})
    _display(plt, false)
    GLWindow.poll_glfw()
    if Base.n_avail(Reactive._messages) > 0
        Reactive.run_till_now()
    end
    yield()
    GLWindow.render_frame(GLWindow.rootscreen(plt.o))
    GLWindow.swapbuffers(plt.o)
    buff = GLWindow.screenbuffer(plt.o)
    png = map(RGB{U8}, buff)
    FileIO.save(FileIO.Stream(FileIO.DataFormat{:PNG}, io), png)
end

function gl_image(img, kw_args)
    rect = kw_args[:primitive]
    kw_args[:primitive] = GeometryTypes.SimpleRectangle{Float32}(rect.x, rect.y, rect.w, rect.h)
    visualize(img, Style(:default), kw_args)
end

function handle_segment(lines, line_segments, points::Vector{P}, segment) where P
    (isempty(segment) || length(segment) < 2) && return
    if length(segment) == 2
         append!(line_segments, view(points, segment))
    elseif length(segment) == 3
        p = view(points, segment)
        push!(line_segments, p[1], p[2], p[2], p[3])
    else
        append!(lines, view(points, segment))
        push!(lines, P(NaN))
    end
end

function gl_lines(points, kw_args)
    result = []
    isempty(points) && return result
    P = eltype(points)
    lines = P[]
    line_segments = P[]
    last = 1
    for (i,p) in enumerate(points)
        if isnan(p) || i==length(points)
            _i = isnan(p) ? i-1 : i
            handle_segment(lines, line_segments, points, last:_i)
            last = i+1
        end
    end
    if !isempty(lines)
        pop!(lines) # remove last NaN
        push!(result, visualize(lines, Style(:lines), kw_args))
    end
    if !isempty(line_segments)
        push!(result, visualize(line_segments, Style(:linesegment), kw_args))
    end
    return result
end

function gl_shape(plotattributes, kw_args)
    points = Plots.extract_points(plotattributes)
    result = []
    for rng in iter_segments(plotattributes[:x], plotattributes[:y])
        ps = points[rng]
        meshes = gl_poly(ps, kw_args)
        append!(result, meshes)
    end
    result
end

function gl_scatter(points, kw_args)
    prim = get(kw_args, :primitive, GeometryTypes.Circle)
    if isa(prim, GLNormalMesh)
        if haskey(kw_args, :model)
            p = get(kw_args, :perspective, eye(GeometryTypes.Mat4f0))
            kw_args[:scale] = GLAbstraction.const_lift(kw_args[:model], kw_args[:scale], p) do m, sc, p
                s  = Vec3f0(m[1,1], m[2,2], m[3,3])
                ps = Vec3f0(p[1,1], p[2,2], p[3,3])
                r  = sc ./ (s .* ps)
                r
            end
        end
    else # 2D prim
        kw_args[:scale] = to_vec(Vec2f0, kw_args[:scale])
    end

    if haskey(kw_args, :stroke_width)
        s = Reactive.value(kw_args[:scale])
        sw = kw_args[:stroke_width]
        if sw*5 > _cycle(Reactive.value(s), 1)[1] # restrict marker stroke to 1/10th of scale (and handle arrays of scales)
            kw_args[:stroke_width] = s[1] / 5f0
        end
    end
    kw_args[:scale_primitive] = false
    if isa(prim, String)
        kw_args[:position] = points
        if !isa(kw_args[:scale], Vector) # if not vector, we can assume it's relative scale
            kw_args[:relative_scale] = kw_args[:scale]
            delete!(kw_args, :scale)
        end
        return visualize(prim, Style(:default), kw_args)
    end

    visualize((prim, points), Style(:default), kw_args)
end

function gl_poly(points, kw_args)
    last(points) == first(points) && pop!(points)
    polys = GeometryTypes.split_intersections(points)
    result = []
    for poly in polys
        mesh = GLNormalMesh(poly) # make polygon
        if !isempty(GeometryTypes.faces(mesh)) # check if polygonation has any faces
            push!(result, GLVisualize.visualize(mesh, Style(:default), kw_args))
        else
            @warn("Couldn't draw the polygon: $points")
        end
    end
    result
end

function gl_surface(x,y,z, kw_args)
    if isa(x, AbstractRange) && isa(y, AbstractRange)
        main = z
        kw_args[:ranges] = (x, y)
    else
        if isa(x, AbstractMatrix) && isa(y, AbstractMatrix)
            main = map(s->map(Float32, s), (x, y, z))
        elseif isa(x, AbstractVector) || isa(y, AbstractVector)
            x = Float32[x[i] for i = 1:size(z,1), j = 1:size(z,2)]
            y = Float32[y[j] for i = 1:size(z,1), j = 1:size(z,2)]
            main = (x, y, map(Float32, z))
        else
            error("surface: combination of types not supported: $(typeof(x)) $(typeof(y)) $(typeof(z))")
        end
        if get(kw_args, :wireframe, false)
            points = map(Point3f0, zip(vec(x), vec(y), vec(z)))
            faces = Cuint[]
            idx = (i,j) -> CartesianIndices(size(z), i, j) - 1
            for i=1:size(z,1), j=1:size(z,2)

                i < size(z,1) && push!(faces, idx(i, j), idx(i+1, j))
                j < size(z,2) && push!(faces, idx(i, j), idx(i, j+1))

            end
            color = get(kw_args, :stroke_color, RGBA{Float32}(0,0,0,1))
            kw_args[:color] = color
            kw_args[:thickness] = Float32(get(kw_args, :stroke_width, 1f0))
            kw_args[:indices] = faces
            delete!(kw_args, :stroke_color)
            delete!(kw_args, :stroke_width)

            return visualize(points, Style(:linesegment), kw_args)
        end
    end
    return visualize(main, Style(:surface), kw_args)
end

function gl_contour(x, y, z, kw_args)
    if kw_args[:fillrange] != nothing

        delete!(kw_args, :intensity)
        I = GLVisualize.Intensity{Float32}
        main = [I(z[j,i]) for i=1:size(z, 2), j=1:size(z, 1)]
        return visualize(main, Style(:default), kw_args)

    else
        h = kw_args[:levels]
        T = eltype(z)
        levels = Contour.contours(map(T, x), map(T, y), z, h)
        result = Point2f0[]
        zmin, zmax = get(kw_args, :limits, Vec2f0(ignorenan_extrema(z)))
        cmap = get(kw_args, :color_map, get(kw_args, :color, RGBA{Float32}(0,0,0,1)))
        colors = RGBA{Float32}[]
        for c in levels.contours
            for elem in c.lines
                append!(result, elem.vertices)
                push!(result, Point2f0(NaN32))
                col = GLVisualize.color_lookup(cmap, c.level, zmin, zmax)
                append!(colors, fill(col, length(elem.vertices) + 1))
            end
        end
        kw_args[:color] = colors
        kw_args[:color_map] = nothing
        kw_args[:color_norm] = nothing
        kw_args[:intensity] = nothing
        return visualize(result, Style(:lines),kw_args)
    end
end

function gl_heatmap(x,y,z, kw_args)
    get!(kw_args, :color_norm, Vec2f0(ignorenan_extrema(z)))
    get!(kw_args, :color_map, Plots.make_gradient(cgrad()))
    delete!(kw_args, :intensity)
    I = GLVisualize.Intensity{Float32}
    heatmap = I[z[j,i] for i=1:size(z, 2), j=1:size(z, 1)]
    tex = GLAbstraction.Texture(heatmap, minfilter=:nearest)
    kw_args[:stroke_width] = 0f0
    kw_args[:levels] = 1f0
    visualize(tex, Style(:default), kw_args)
end

"""
Ugh, so much special casing (╯°□°）╯︵ ┻━┻
"""
function label_scatter(plotattributes, w, ho)
    kw = KW()
    extract_stroke(plotattributes, kw)
    extract_marker(plotattributes, kw)
    kw[:scale] = Vec2f0(w/2)
    kw[:offset] = Vec2f0(-w/4)
    if haskey(kw, :intensity)
        cmap = kw[:color_map]
        norm = kw[:color_norm]
        kw[:color] = GLVisualize.color_lookup(cmap, kw[:intensity][1], norm)
        delete!(kw, :intensity)
        delete!(kw, :color_map)
        delete!(kw, :color_norm)
    else
        color = get(kw, :color, nothing)
        kw[:color] = isa(color, Array) ? first(color) : color
    end
    strcolor = get(kw, :stroke_color, RGBA{Float32}(0,0,0,0))
    kw[:stroke_color] = isa(strcolor, Array) ? first(strcolor) : strcolor
    p = get(kw, :primitive, GeometryTypes.Circle)
    if isa(p, GLNormalMesh)
        bb = GeometryTypes.AABB{Float32}(GeometryTypes.vertices(p))
        bbw = GeometryTypes.widths(bb)
        if isapprox(bbw[3], 0)
            bbw = Vec3f0(bbw[1], bbw[2], 1)
        end
        mini = minimum(bb)
        m = GLAbstraction.translationmatrix(-mini)
        m *= GLAbstraction.scalematrix(1 ./ bbw)
        kw[:primitive] = m * p
        kw[:scale] = Vec3f0(w/2)
        delete!(kw, :offset)
    end
    if isa(p, Array)
        kw[:primitive] = GeometryTypes.Circle
    end
    GL.gl_scatter(Point2f0[(w/2, ho)], kw)
end

function make_label(sp, series, i)
    GL = Plots
    w, gap, ho = 20f0, 5f0, 5
    result = []
    plotattributes = series.plotattributes
    st = plotattributes[:seriestype]
    kw_args = KW()
    if (st in (:path, :path3d, :straightline)) && plotattributes[:linewidth] > 0
        points = Point2f0[(0, ho), (w, ho)]
        kw = KW()
        extract_linestyle(plotattributes, kw)
        append!(result, GL.gl_lines(points, kw))
        if plotattributes[:markershape] != :none
            push!(result, label_scatter(plotattributes, w, ho))
        end
    elseif st in (:scatter, :scatter3d) #|| plotattributes[:markershape] != :none
        push!(result, label_scatter(plotattributes, w, ho))
    else
        extract_c(plotattributes, kw_args, :fill)
        if isa(kw_args[:color], AbstractVector)
            kw_args[:color] = first(kw_args[:color])
        end
        push!(result, visualize(
            GeometryTypes.SimpleRectangle(-w/2, ho-w/4, w/2, w/2),
            Style(:default), kw_args
        ))
    end
    labeltext = if isa(series[:label], Array)
        i += 1
        series[:label][i]
    else
        series[:label]
    end
    ft = legendfont(sp)
    font = Plots.Font(ft.family, ft.pointsize, :left, :bottom, 0.0, ft.color)
    xy = Point2f0(w+gap, 0.0)
    kw = Dict(:model => text_model(font, xy), :scale_primitive=>false)
    extract_font(font, kw)
    t = PlotText(labeltext, font)
    push!(result, glvisualize_text(xy, t, kw))
    GLAbstraction.Context(result...), i
end

function generate_legend(sp, screen, model_m)
    legend = GLAbstraction.Context[]
    if sp[:legend] != :none
        i = 0
        for series in series_list(sp)
            should_add_to_legend(series) || continue
            result, i = make_label(sp, series, i)
            push!(legend, result)
        end
        if isempty(legend)
            return
        end
        list = visualize(legend, gap=Vec3f0(0,5,0))
        bb = GLAbstraction._boundingbox(list)
        wx,wy,_ = GeometryTypes.widths(bb)
        xmin, _ = Plots.axis_limits(sp[:xaxis])
        _, ymax = Plots.axis_limits(sp[:yaxis])
        area = map(model_m) do m
            p = m * GeometryTypes.Vec4f0(xmin, ymax, 0, 1)
            h = round(Int, wy)+20
            w = round(Int, wx)+20
            x,y = round(Int, p[1])+30, round(Int, p[2]-h)-30
            GeometryTypes.SimpleRectangle(x, y, w, h)
        end
        sscren = GLWindow.Screen(
            screen, area = area,
            color = sp[:background_color_legend],
            stroke = (2f0, RGBA(0.3, 0.3, 0.3, 0.9))
        )
        GLAbstraction.translate!(list, Vec3f0(10,10,0))
        GLVisualize._view(list, sscren, camera=:fixed_pixel)
    end
    return
end


text:
import CSV
using DataFrames: DataFrame, Not, All
include("utils.jl")

function load_results()
    if isfile("result.csv")
        df = CSV.read("result.csv")
        df.date
        Dict((df.expID .=> df.testID)
                        .=>
                        zip(df.prec, df.recall, df.shd,
                            df.time, df.date))
    else
        @warn "No result.csv"
        Dict()
    end
end

function parse_spec_string(s)
    re = r"d=(\d+)_k=(\d+)_gtype=(.*)_noise=(.*)_mat=([^_]*)(_mec=MLP)?(_raw)?"
    m = match(re, s)
    if isnothing(m) @show s end
    d = parse(Int, m.captures[1])
    k = parse(Int, m.captures[2])
    gtype, noise, mat, mec, raw = m.captures[3:end]
    if isnothing(mec)
        mec = "Linear"
    else
        mec = "MLP"
    end
    if isnothing(raw)
        raw = false
    else
        raw = true
    end
    # FIXME nothing seems to match any expressions, see
    # https://github.com/kmsquire/Match.jl/issues/60
    #
    # mec = @match mec begin
    #     nothing => "Linear"
    #     # CAUTION This is substring, thus this won't match
    #     # var::String => "MLP"
    #     _ => "MLP"
    # end
    return d, k, gtype, noise, mat, mec, raw
end

# FIXME these type definitions do not seem to be allowed inside function scope
MInt = Union{Missing, Int}
MString = Union{Missing, String}

function results2df(results)
    df = DataFrame(model=String[],
                   train_d=MInt[], train_k=MInt[],
                   train_gtype=MString[],
                   train_noise=MString[],
                   train_mat=MString[],
                   train_mec=MString[],

                   test_d=Int[], test_k=Int[],
                   test_gtype=String[],
                   test_noise=String[],
                   test_mat=String[],
                   test_mec=MString[],
                   raw=Bool[],

                   prec=Float32[], recall=Float32[], shd=Float32[],
                   time=Float32[])
    for item in results
        # 1. find the model name
        m = match(r"(.*)-d=.*", item[1][1])
        if isnothing(m)
            model = item[1][1]
            train_d = missing
            train_k = missing
            train_gtype = missing
            train_noise = missing
            train_mat = missing
            train_mec = missing
        else
            # model = item[1][1]
            model = m.captures[1]
            (train_d, train_k, train_gtype, train_noise, train_mat,
             train_mec, _) = parse_spec_string(item[1][1])
        end

        (test_d, test_k, test_gtype, test_noise, test_mat,
         test_mec, raw) = parse_spec_string(item[1][2])

        prec = round3(item[2][1]) * 100
        recall = round3(item[2][2]) * 100
        shd = round(item[2][3], digits=1)
        t = round(item[2][4], digits=1)
        # prec, recall, shd = round3.(item[2][1:end-1])

        push!(df, (model,
                   train_d, train_k, train_gtype, train_noise, train_mat, train_mec,
                   test_d, test_k, test_gtype, test_noise, test_mat, test_mec, raw,
                   prec, recall, shd, t))
    end
    return df
end

function table_raw(df)
    selector = (
        # (df.raw .== 1) .&
        (df.train_gtype .== "SF")
        .& in.(df.test_d, Ref([10, 20, 50, 100, 200, 300, 400]))
        # .& in.(df.train_d, Ref([10, 20, 50, 100, 200, 300, 400]))
    )
    sort(df[selector,
            All(:model, :train_d, :prec, :recall, :shd, :time, :test_d)],
         [:train_d])
end

function table_transfer_k(df)
    # - train on k=1, test on k=1,2,4
    selector = ((df.model .== "deep-EQ")
                .& (df.test_noise .== "Gaussian")
                .& (df.train_gtype .== "SF")
                .& (df.test_gtype .== "SF")
                .& (df.train_mec .== "Linear")
                .& (df.test_mec .== "Linear")
                .& (df.train_mat .== "COR")
                .& (df.test_mat .== "COR"))
    sort(df[selector,
            Not(All(:train_k, :train_noise,
                    :test_d,
                    :train_mec, :test_mec, :test_noise,
                    :train_mat, :test_mat))],
         [:train_gtype, :model, :test_gtype, :train_d, :test_k])
end

function table_transfer_noise(df)
    # - train on Gaussian, test on Exp, Gumbel
    selector = ((df.model .== "deep-EQ")
                .& in.(df.train_gtype, Ref(["SF", "ER"]))
                .& (df.train_gtype .== df.test_gtype)
                .& (df.test_k .== 1)
                .& (df.train_mec .== "Linear")
                .& (df.test_mec .== "Linear")
                .& (df.train_mat .== "COR")
                .& (df.test_mat .== "COR"))
    sort(df[selector,
            Not(All(:train_k,
                    # CAUTION I need to show this to tell that the training and
                    # testing noise are different :train_noise,
                    :test_d, :test_k,
                    :train_mec, :test_mec,
                    :test_gtype,
                    :train_mat, :test_mat))],
         [:train_gtype, :model, :train_d, :test_noise])
end

function table_transfer_gtype(df)
    # - train on ER, test on SF
    selector = ((df.model .== "deep-EQ")
                # TODO ER2 ER4 SF2 SF4
                .& in.(df.train_gtype, Ref(["ER", "SF"]))
                .& (df.train_gtype .!= df.test_gtype)
                .& (df.test_noise .== "Gaussian")
                .& (df.test_k .== 1)
                .& (df.train_mec .== "Linear")
                .& (df.test_mec .== "Linear")
                .& (df.train_mat .== "COR")
                .& (df.test_mat .== "COR"))
    sort(df[selector,
            Not(All(:train_k, :train_noise,
                    :test_d, :test_k,
                    :train_mec, :test_mec, :test_noise,
                    :train_mat, :test_mat))],
         [:train_gtype, :model, :test_gtype, :train_d])
end

function table_cmp_baseline(df)
    # comparison of EQ with FC and CNN (and random) baselines
    # only show d=10,15,20, SF and ER graphs, k=1
    sel1 = in.(df.model, Ref(["FC", "EQ", "CNN"]))
    # sel2 = in.(df.train_d, Ref([10, 15, 20,8,16,32]))
    sel2 = in.(df.train_d, Ref([10, 20, 50, 100]))
    selector = (sel1 .& sel2
                .& (df.train_gtype .== df.test_gtype)
                .& in.(df.train_gtype, Ref(["SF"]))
                .& (df.test_noise .== "Gaussian")
                .& (df.train_mec .== "Linear")
                .& (df.test_mec .== "Linear")
                .& (df.train_mat .== "COR")
                .& (df.test_mat .== "COR")
                .& (df.test_k .== 1))

    sort(df[selector,
            All(:model, :train_d, :prec, :recall, :shd)],
         [:model, :train_d])
end

function table_ER24(df)
    # - ER2 ER4 SF2 SF4
    # TODO CNN, FC data for these graphs
    selector = ((df.model .== "deep-EQ")
                .& in.(df.train_gtype, Ref(["SF", "SF2", "SF4", "ER", "ER2", "ER4"]))
                .& (df.train_gtype .== df.test_gtype)
                .& (df.test_k .== 1)
                .& (df.test_noise .== "Gaussian")
                .& (df.train_mec .== "Linear")
                .& (df.test_mec .== "Linear")
                .& (df.train_mat .== "COR")
                .& (df.test_mat .== "COR"))
    sort(df[selector,
            Not(All(:train_k, :train_noise,
                    :test_d, :test_k,
                    :train_mec, :test_mec, :test_noise,
                    :train_mat, :test_mat))],
         [:train_gtype, :model, :test_gtype, :train_d])
end

function table_ensemble(df)
    selector = (in.(df.model, Ref(["deep-EQ-ER-ensemble",
                                   "deep-EQ-SF-ensemble"]))
                .& (df.test_k .== 1)
                .& (df.test_mec .== "Linear")
                .& (df.test_mat .== "COR"))
    sort(df[selector,
            Not(All(:train_k, :train_noise,
                    :test_k,
                    r"train_",
                    :train_mec, :test_mec, :test_noise,
                    :train_mat, :test_mat))],
         [:model, :test_gtype, :test_d])
end

function test()
    # DONE close comparison between flat-CNN and deep-EQ
    # header=map((x)->replace(x,"_"=>"/"),
    #            names(tbl_baseline))
    # quotestrings=true,
    table_cmp_baseline(df)
    CSV.write("results/cmp-baseline.csv", table_cmp_baseline(df))
    CSV.write("results/er4.csv", table_ER24(df))

    # DONE add MLP testing results

    # 1.5 raw
    CSV.write("results/raw.csv", table_raw(df))

    # 2. transfering
    #
    # - train on k=1, test on k=1,2,4
    CSV.write("results/transfer_k.csv", table_transfer_k(df))
    # - train on Gaussian, test on Exp, Gumbel
    CSV.write("results/transfer_noise.csv", table_transfer_noise(df))
    # - train on ER, test on SF
    CSV.write("results/transfer_gtype.csv", table_transfer_gtype(df))

    # 3. the ensemble models. The training setting is in model name
    select_ensemble_df(df, "deep-EQ-ER-ensemble")
    select_ensemble_df(df, "deep-EQ-SF-ensemble")

    select_ensemble_df(df, "deep-EQ-ER-COV-ensemble")
    select_ensemble_df(df, "deep-EQ-SF-COV-ensemble")

    CSV.write("results/ensemble.csv", table_ensemble(df))
end

max_stars_repo_name: lihebi/DAG-EQ

max_issues_repo_name: lihebi/DAG-EQ

max_forks_repo_name: lihebi/DAG-EQ


text:
include("Tests/Tests.jl")

using Test, JuMP, SDPA
factory = with_optimizer(SDPA.Optimizer)
config = MOI.Test.TestConfig(atol=1e-5, rtol=1e-5, query=false)
@testset "Linear" begin
    Tests.linear_test(factory, config)
end
@testset "SDP" begin
    Tests.sd_test(factory, config)
end


text:
module ALPA
using StatsBase

export addnode!, delnode!, addedge!, deledge!, ne,
    lpa_addnode!, lpa_deletenode!, lpa_addedge!, lpa_deleteedge!, loadgml, savegml, modularity, graphdiff, Graph, NeighComm, NodeStatus, getctnf
# package code goes here
include("wgraph.jl")
include("wevolve.jl")
include("graphio.jl")
include("utils.jl")
# package code goes here
N = 10000
c = NeighComm(collect(1:N), fill(-1,N), 1)
ns = NodeStatus(collect(1:N), fill(false,N), 1)
random_order = Array(Int,N)
#m = zeros(Int,N)
a = IntSet()

lpa_deletenode!(g, m, u) = lpa_deletenode!(g, m, c, u, a, random_order, ns)

lpa_addedge!(g,m,u,v) = lpa_addedge!(g,m,c,u,v,a,random_order,ns)

lpa_deleteedge!(g,m,u,v) = lpa_deleteedge!(g,m,c,u,v,a,random_order,ns)
end # module


text:
"""Abstract type for nodes in AST tree
"""
abstract type Node end

max_stars_repo_name: Sohva/TensorDSL

max_issues_repo_name: Sohva/TensorDSL

max_forks_repo_name: Sohva/TensorDSL


text:
using HypothesisTests, Test
using HypothesisTests: default_tail

@testset "T-test" begin

@testset "One sample" begin
	@test pvalue(OneSampleTTest(-5:5)) == 1

	tst = OneSampleTTest(-5:10)
	@test abs(pvalue(tst) - 0.0530) <= 1e-4
	@test abs(pvalue(tst; tail=:left) - 0.9735) <= 1e-4
	@test abs(pvalue(tst; tail=:right) - 0.0265) <= 1e-4
	@test default_tail(tst) == :both
	show(IOBuffer(), tst)

	tst = OneSampleTTest(mean(-5:10), std(-5:10), 16)
	@test abs(pvalue(tst) - 0.0530) <= 1e-4

	@test all(abs.([confint(tst)...;] - [-0.0369, 5.0369]) .<= 1e-4)
	@test all(abs.([confint(tst, level=0.9)...;] - [0.4135, 4.5865]) .<= 1e-4)
	c = confint(tst; tail=:left)
	@test c[1] == -Inf
	@test abs(c[2] - 4.5865) .<= 1e-4
	c = confint(tst; tail=:right)
	@test abs(c[1] - 0.4135) .<= 1e-4
	@test c[2] == Inf
	show(IOBuffer(), tst)

	tst = OneSampleTTest(-10:5)
	@test abs(pvalue(tst) - 0.0530) <= 1e-4
	@test abs(pvalue(tst; tail=:left) - 0.0265) <= 1e-4
	@test abs(pvalue(tst; tail=:right) - 0.9735) <= 1e-4
	@test all(abs.([confint(tst)...] - [-5.0369, 0.0369]) .<= 1e-4)
	@test abs.(confint(tst; tail=:left)[2] - (-0.4135)) .<= 1e-4
	@test abs.(confint(tst; tail=:right)[1] - (-4.5865)) .<= 1e-4
	show(IOBuffer(), tst)
end

@testset "Paired" begin
	@test abs(pvalue(OneSampleTTest([1, 1, 2, 1, 0], [0, 1, 1, 1, 0])) - 0.1778) <= 1e-4
end

@testset "Two sample" begin
	# From http://en.wikipedia.org/w/index.php?title=Student%27s_t-test&oldid=526762741

	a1 = [30.02, 29.99, 30.11, 29.97, 30.01, 29.99]
	a2 = [29.89, 29.93, 29.72, 29.98, 30.02, 29.98]

	tst = EqualVarianceTTest(a1, a2)
	@test tst.df == 10
	@test abs(tst.t - 1.959) <= 1e-3
	@test abs(pvalue(tst) - 0.078) <= 1e-3
	@test all(abs.([confint(tst)...] - [-0.0131, 0.2031]) .<= 1e-4)
	@test default_tail(tst) == :both
	show(IOBuffer(), tst)

	tst = UnequalVarianceTTest(a1, a2)
	@test abs(tst.df - 7.03) <= 0.01
	@test abs(tst.t - 1.959) <= 1e-3
	@test abs(pvalue(tst) - 0.091) <= 1e-3
	@test all(abs.([confint(tst)...] - [-0.0196, 0.2096]) .<= 1e-4)
	@test default_tail(tst) == :both
	show(IOBuffer(), tst)
end
end


text:
module ImageSegmentation

import Base: show

using LinearAlgebra, Statistics
using DataStructures, StaticArrays, ImageCore, ImageFiltering, ImageMorphology, LightGraphs, SimpleWeightedGraphs, RegionTrees, Distances, StaticArrays, Clustering, MetaGraphs
import Clustering: kmeans, fuzzy_cmeans

include("compat.jl")
include("core.jl")
include("region_growing.jl")
include("felzenszwalb.jl")
include("fast_scanning.jl")
include("watershed.jl")
include("region_merging.jl")
include("meanshift.jl")
include("clustering.jl")
include("merge_segments.jl")

export
    #accessor methods
    labels_map,
    segment_labels,
    segment_pixel_count,
    segment_mean,

    # methods
    seeded_region_growing,
    unseeded_region_growing,
    felzenszwalb,
    fast_scanning,
    watershed,
    hmin_transform,
    region_adjacency_graph,
    remove_segment,
    remove_segment!,
    prune_segments,
    region_tree,
    region_splitting,
    meanshift,
    kmeans,
    fuzzy_cmeans,
    merge_segments,

    # types
    SegmentedImage,
    ImageEdge

@deprecate rem_segment  remove_segment
@deprecate rem_segment! remove_segment!

end # module


text:
# ---
# title: 670. Maximum Swap
# id: problem670
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array, Math
# link: <https://leetcode.com/problems/maximum-swap/description/>
# hidden: true
# ---
#
# Given a non-negative integer, you could swap two digits **at most** once to
# get the maximum valued number. Return the maximum valued number you could get.
#
# **Example 1:**
#
#
#
#     Input: 2736
#     Output: 7236
#     Explanation: Swap the number 2 and the number 7.
#
#
# **Example 2:**
#
#
#
#     Input: 9973
#     Output: 9973
#     Explanation: No swap.
#
#
# **Note:**
#
#   1. The given number is in the range [0, 108]
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
# module for storing q value at each state

using Random
using StatsBase

mutable struct StateInfo
  q
  epsilon

  # init
  function StateInfo(action_num=0; epsilon=0.3)
    self = new()

    self.q = zeros(action_num)
    self.epsilon = epsilon

    return self
  end
end

# get index of action which has maximum q value
# greedy policy
function greedy(self::StateInfo)
  return argmax(self.q)
end

# select another action following probability epsilon
function epsilon_greedy(self::StateInfo)
  if rand(1)[1] < self.epsilon
    return rand(1:length(self.q), 1)[1]
  else
    return greedy(self)
  end
end

# call greedy()
function _pi(self::StateInfo)
  return epsilon_greedy(self)
end

# get max q value
function max_q(self::StateInfo)
  return maximum(self.q)
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
using DReal
using Base.Test
reset_ctx!()

A = Var(Bool,"aa")
B = Var(Bool,"bb")
C = Var(Bool,"cc")

formula = add!((A & B) | C)
@show is_satisfiable()
@show a, b, c = model(A, B, C)
# @test (a & b) | c


text:
#=
   Test of credible interval in Turing.jl

   This is a port of the WebPPL model in
   https://mhtess.github.io/bdappl/chapters/03-simpleModels.html
   """
   A simple generative process
   Imagine we are investigating the origins of prosocial behavior: Some scientists believe humans are
   born with a self-interested tendency and only acquire prosociality through instruction (perhaps mediated
   by language). Others believe altruism in innate; as soon as children are physically and cognitively
   able to help, they will. To investigate this, we are going to put pre-verbal infants in a context where
   a person is struggling to complete their goal (e.g., as in this classic study
   [https://www.youtube.com/watch?v=kfGAen6QiUE Toddler Altruism]).

   We plan to conduct this experiment on 20 infants. Each child may or may not help; what we are interested
   in is the population’s propensity to help. A simple generative model of this task would be that infants'
   propensity to help is like the bias of a coin, and whether or not they help in the moment is the
   result of flipping that coin.
   """

   credible_interval is defined in jl_utils.jl

   Output of the model:
"""
Summary Statistics
                parameters      mean       std   naive_se      mcse         ess      rhat   ess_per_sec
                    Symbol   Float64   Float64    Float64   Float64     Float64   Float64       Float64

         priorDistribution    0.4966    0.2861     0.0029    0.0054   3065.5383    1.0000     1283.7262
        propensity_to_help    0.7260    0.0937     0.0009    0.0013   5215.2559    0.9999     2183.9430
       posteriorPredictive   14.5179    2.6933     0.0269    0.0323   6745.5678    1.0000     2824.7771
  prior_propensity_to_help    0.4966    0.2861     0.0029    0.0054   3065.5383    1.0000     1283.7262
           priorPredictive    9.9341    6.0131     0.0601    0.1088   3370.6685    0.9999     1411.5027

Quantiles
                parameters      2.5%     25.0%     50.0%     75.0%     97.5%
                    Symbol   Float64   Float64   Float64   Float64   Float64

         priorDistribution    0.0203    0.2559    0.4953    0.7390    0.9736
        propensity_to_help    0.5260    0.6663    0.7309    0.7937    0.8867
       posteriorPredictive    9.0000   13.0000   15.0000   17.0000   19.0000
  prior_propensity_to_help    0.0203    0.2559    0.4953    0.7390    0.9736
           priorPredictive    0.0000    5.0000   10.0000   15.0000   20.0000

Credible interval for posteriorPredictive with mass 0.9: (10.000000..18.000000)
Credible interval for posteriorPredictive with mass 0.95: (9.000000..19.000000)
Credible interval for posteriorPredictive with mass 0.99: (8.000000..20.000000)
Credible interval for posteriorPredictive with mass 0.1: (13.000000..13.000000)
Credible interval for posteriorPredictive with mass 0.01: (9.000000..9.000000)
"""

=#

using Turing
using Random
include("jl_utils.jl")

@model function credible_interval_test(k=15,n=20)
    priorDistribution ~ Uniform(0,1)

    # propensity to help, or true population proportion who would help
    propensity_to_help ~ Uniform(0,1)

    # Observed k children who help
    # Assuming each child's response is independent of each other
    k ~ Binomial(n, propensity_to_help)

    # predict what the next n will do
    posteriorPredictive ~ Binomial(n, propensity_to_help)

    # duplicate model structure and parameters but omit observe
    prior_propensity_to_help ~ Dirac(priorDistribution)
    priorPredictive ~ Binomial(n, prior_propensity_to_help)

end

# observed data
k = 15 # number of children who help
n = 20 # total number of children
model = credible_interval_test(k,n)

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 10_000)
# chns = sample(model, MH(), 10_000)
# chns = sample(model, PG(15), 10_000)
# chns = sample(model, IS(), 10_000)
chns = sample(model, SMC(), 10_000)
# chns = sample(model, SMC(), MCMCThreads(), 10_000, 4)

display(chns)

credible_interval(chns, "posteriorPredictive",0.90)
credible_interval(chns, "posteriorPredictive",0.95)
credible_interval(chns, "posteriorPredictive",0.99)
credible_interval(chns, "posteriorPredictive",0.1)
credible_interval(chns, "posteriorPredictive",0.01)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
using QMTK
using QMTK.Consts.Pauli
using Compat.Test

chain = Chain(Fixed, 4)
t = kronsum(bonds(chain, 1)) do i, j
    :(σ₁[$i] * σ₁[$j])
end
h = kronprod(σ₁, σ₁, σ₀, σ₀) + kronprod(σ₀, σ₁, σ₁, σ₀) + kronprod(σ₀, σ₀, σ₁, σ₁)

@test t == h

chain = Chain(Periodic, 4)
t = kronsum(sites(chain)) do i
    :(σ₁[$i])
end

h = kronprod(σ₁, σ₀, σ₀, σ₀) + kronprod(σ₀, σ₁, σ₀, σ₀) +
    kronprod(σ₀, σ₀, σ₁, σ₀) + kronprod(σ₀, σ₀, σ₀, σ₁)
@test t == h


text:
#=
bfgs.jl

    Quasi-Newton type update step routines using the
    Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm, that uses approximations
    to the Hessian using past approximations as well as the gradient.

@author: Quint Wiersma <q.wiersma@vu.nl>

@date: 2022/01/24
=#

# BFGS state
mutable struct BFGSState{Tv, Tf, Tm}
    x::Tv       # current state
    x_prev::Tv  # previous state
    ∇f_prev::Tv # previous gradient of f(x)
    f_prev::Tf  # previous f(x)
    s::Tv       # change in state
    y::Tv       # change in gradient
    p::Tv       # search direction
    Hi::Tm      # inverse Hessian approx
    u::Tv       # buffer
end

"""
    update_state!(state, ls, f, ∇f!)

Update state using the BFGS inverse Hessian approximation for the search
direction with an inexact backtracking linesearch, storing the result in `x`.

#### Arguments
  - `state::BFGSState`  : state variables
  - `ls::BackTrack`     : line search parameters
  - `f::Function`       : ``f(x)``
  - `∇f!::Function`     : gradient of `f`
"""
function update_state!(state::BFGSState, ls::BackTrack, f::Function, ∇f!::Function)
    # Infer type
    T= eltype(state.p)

    # Current gradient
    ∇f!(state.∇f_prev, state.x)

    # Search direction
    mul!(state.p, state.Hi, state.∇f_prev, T(-1), zero(T))

    # Backtracking line Search
    f_x= backtrack!(ls, state.x, state.x_prev, state.p, state.f_prev,
                    state.∇f_prev, f)

    # Store current objective function value
    state.f_prev= f_x

    return nothing
end

"""
    update_h!(state)

Update the inverse Hessian approximation using the BFGS algorithm, storing the
result in `Hi`.

#### Arguments
  - `state::BFGSState`  : state variables
"""
function update_h!(state::BFGSState)
    # Curvature condition
    curv= dot(state.y, state.s)

    if curv > 0
        # Cache
        mul!(state.u, state.Hi, state.y)

        # Constants
        α= (curv + dot(state.u, state.y)) * inv(curv * curv)
        β= inv(curv)

        # BFGS update, using Sherman-Morrison-Woodbury
        n= length(state.x)
        # Small dimensions: loop unfolding
        if n < 50
            for i in 1:n
                s_i= state.s[i]
                u_i= state.u[i]
                @simd for j in 1:n
                    @inbounds state.Hi[i,j]+= α * s_i * state.s[j] -
                                                β * (u_i * state.s[j] +
                                                        state.u[j] * s_i)
                end
            end
        # Large dimensions: BLAS routines
        else
            BLAS.ger!(α, state.s, state.s, state.Hi)
            BLAS.ger!(-β, state.u, state.s, state.Hi)
            BLAS.ger!(-β, state.s, state.u, state.Hi)
        end
    end

    return nothing
end

"""
    bfgs!(state, ls, f, ∇f!)

Update state and inverse Hessian approximation using the BFGS algorithm.

#### Arguments
  - `state::BFGSState`  : state variables
  - `ls::BackTrack`     : line search parameters
  - `f::Function`       : ``f(x)``
  - `∇f!::Function`     : gradient of `f`
"""
function bfgs!(state::BFGSState, ls::BackTrack, f::Function, ∇f!::Function)
    # Update state
    update_state!(state, ls, f, ∇f!)

    # Change in state
    @. state.s= ls.α * state.p

    # Change in gradient
    ∇f!(state.y, state.x)
    @. state.y= state.y - state.∇f_prev

    # Update inverse Hessian approx
    update_h!(state)

    return nothing
end


text:
using JuMP

### Place optimization problem definition and solve in a function
function solveOptProb1(a,b)
	m = Model()
	@variable(m, 0 <= x <= 2 )
	@variable(m, 0 <= y <= 30 )

	@objective(m, Max, a*x + 3*y )
	@constraint(m, 1x + 5y <= b )

	status = solve(m)

	return getobjectivevalue(m), [getvalue(x); getvalue(y)], status

end

### Test function

(obj, z, status) = solveOptProb1(1,1)

### Iterate over possible values of a and b
aOpt = 0:5
bOpt = 1:3

aN = length(aOpt)
bN = length(bOpt)

objValues = zeros(aN,bN)
zValues = zeros(2,aN,bN)

for i in 1:length(aOpt)
	for j in 1:length(bOpt)
		(objValues[i,j], zValues[:,i,j], status) = solveOptProb1(aOpt[i],bOpt[j])
	end
end

println("objValues = ")
println(objValues)
println(" ")

println("zValues = ")
println(zValues)
println(" ")


text:
module Experiments

using Reexport
using RecipesBase

@reexport using AccuracyAtTopPrimal
@reexport using BSON
@reexport using CSV
@reexport using DatasetProvider
@reexport using DataFrames
@reexport using Flux
@reexport using Plots

using BenchmarkTools
using DrWatson
using EvalMetrics
using HypothesisTests
using Measurements
using PaperUtils
using Query

using ProgressMeter
using Random
using Statistics

using AccuracyAtTopPrimal: AbstractThreshold
using Base.Iterators: partition
using BSON: load
using DrWatson: datadir
using Flux: gpu, cpu, params
using StatsBase: sample

export Data, Train, Model
export trainmodel, davemodel, loadmodel, loaddata, run_simulations, collect_metrics, selectmetric, selectbest, betterzero, crit_diag, run_benchmarks, summary, restore_path, loadrow, plotroc, wilcoxontest, corelationmatrix, mergecsv

modeldir(args...) = datadir("models", args...)

include("corrmatrix.jl")
include("data_train_model.jl")
include("train.jl")
include("evaluation.jl")

end # module