text:
# https://github.com/jheinen/GR.jl

# significant contributions by @jheinen

import GR
export GR

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

gr_linetype(k) = (auto = 1, solid = 1, dash = 2, dot = 3, dashdot = 4, dashdotdot = -1)[k]

gr_markertype(k) = (
    auto = 1,
    none = -1,
    circle = -1,
    rect = -7,
    diamond = -13,
    utriangle = -3,
    dtriangle = -5,
    ltriangle = -18,
    rtriangle = -17,
    pentagon = -21,
    hexagon = -22,
    heptagon = -23,
    octagon = -24,
    cross = 2,
    xcross = 5,
    (+) = 2,
    x = 5,
    star4 = -25,
    star5 = -26,
    star6 = -27,
    star7 = -28,
    star8 = -29,
    vline = -30,
    hline = -31,
)[k]

gr_halign(k) = (left = 1, hcenter = 2, right = 3)[k]
gr_valign(k) = (top = 1, vcenter = 3, bottom = 5)[k]

const gr_font_family = Dict(
    # compat:
    "times" => 101,
    "helvetica" => 105,
    "courier" => 109,
    "bookman" => 114,
    "newcenturyschlbk" => 118,
    "avantgarde" => 122,
    "palatino" => 126,
    "serif-roman" => 232,
    "sans-serif" => 233,
    # https://gr-framework.org/fonts.html:
    "times roman" => 101,
    "times italic" => 102,
    "times bold" => 103,
    "times bold italic" => 104,
    "helvetica" => 105,
    "helvetica oblique" => 106,
    "helvetica bold" => 107,
    "helvetica bold oblique" => 108,
    "courier" => 109,
    "courier oblique" => 110,
    "courier bold" => 111,
    "courier bold oblique" => 112,
    "symbol" => 113,
    "bookman light" => 114,
    "bookman light italic" => 115,
    "bookman demi" => 116,
    "bookman demi italic" => 117,
    "new century schoolbook roman" => 118,
    "new century schoolbook italic" => 119,
    "new century schoolbook bold" => 120,
    "new century schoolbook bold italic" => 121,
    "avantgarde book" => 122,
    "avantgarde book oblique" => 123,
    "avantgarde demi" => 124,
    "avantgarde demi oblique" => 125,
    "palatino roman" => 126,
    "palatino italic" => 127,
    "palatino bold" => 128,
    "palatino bold italic" => 129,
    "zapf chancery medium italic" => 130,
    "zapf dingbats" => 131,
    "computer modern" => 232,
    "dejavu sans" => 233,
)

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

gr_color(c) = gr_color(c, color_type(c))

gr_color(c, ::Type{<:AbstractRGB}) = UInt32(
    round(UInt, clamp(alpha(c) * 255, 0, 255)) << 24 +
    round(UInt, clamp(blue(c) * 255, 0, 255)) << 16 +
    round(UInt, clamp(green(c) * 255, 0, 255)) << 8 +
    round(UInt, clamp(red(c) * 255, 0, 255)),
)
function gr_color(c, ::Type{<:AbstractGray})
    g = round(UInt, clamp(gray(c) * 255, 0, 255))
    α = round(UInt, clamp(alpha(c) * 255, 0, 255))
    rgba = UInt32(α << 24 + g << 16 + g << 8 + g)
end
gr_color(c, ::Type) = gr_color(RGBA(c), RGB)

function gr_getcolorind(c)
    gr_set_transparency(float(alpha(c)))
    convert(Int, GR.inqcolorfromrgb(red(c), green(c), blue(c)))
end

gr_set_linecolor(c) = GR.setlinecolorind(gr_getcolorind(_cycle(c, 1)))
gr_set_fillcolor(c) = GR.setfillcolorind(gr_getcolorind(_cycle(c, 1)))
gr_set_markercolor(c) = GR.setmarkercolorind(gr_getcolorind(_cycle(c, 1)))
gr_set_bordercolor(c) = GR.setbordercolorind(gr_getcolorind(_cycle(c, 1)))
gr_set_textcolor(c) = GR.settextcolorind(gr_getcolorind(_cycle(c, 1)))
gr_set_transparency(α::Real) = GR.settransparency(clamp(α, 0, 1))
gr_set_transparency(::Nothing) = GR.settransparency(1)
gr_set_transparency(c, α) = gr_set_transparency(α)
gr_set_transparency(c::Colorant, ::Nothing) = gr_set_transparency(c)
gr_set_transparency(c::Colorant) = GR.settransparency(alpha(c))

gr_set_arrowstyle(s::Symbol) = GR.setarrowstyle(
    get(
        (
            simple = 1,
            hollow = 3,
            filled = 4,
            triangle = 5,
            filledtriangle = 6,
            closed = 6,
            open = 5,
        ),
        s,
        1,
    ),
)

gr_set_fillstyle(::Nothing) = GR.setfillintstyle(GR.INTSTYLE_SOLID)
function gr_set_fillstyle(s::Symbol)
    GR.setfillintstyle(GR.INTSTYLE_HATCH)
    GR.setfillstyle(get(((/) = 9, (\) = 10, (|) = 7, (-) = 8, (+) = 11, (x) = 6), s, 9))
end

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

# draw line segments, splitting x/y into contiguous/finite segments
# note: this can be used for shapes by passing func `GR.fillarea`
function gr_polyline(x, y, func = GR.polyline; arrowside = :none, arrowstyle = :simple)
    iend = 0
    n = length(x)
    while iend < n - 1
        # set istart to the first index that is finite
        istart = -1
        for j in (iend + 1):n
            if isfinite(x[j]) && isfinite(y[j])
                istart = j
                break
            end
        end

        if istart > 0
            # iend is the last finite index
            iend = -1
            for j in (istart + 1):n
                if isfinite(x[j]) && isfinite(y[j])
                    iend = j
                else
                    break
                end
            end
        end

        # if we found a start and end, draw the line segment, otherwise we're done
        if istart > 0 && iend > 0
            func(x[istart:iend], y[istart:iend])
            if arrowside in (:head, :both)
                gr_set_arrowstyle(arrowstyle)
                GR.drawarrow(x[iend - 1], y[iend - 1], x[iend], y[iend])
            end
            if arrowside in (:tail, :both)
                gr_set_arrowstyle(arrowstyle)
                GR.drawarrow(x[istart + 1], y[istart + 1], x[istart], y[istart])
            end
        else
            break
        end
    end
end

function gr_polyline3d(x, y, z, func = GR.polyline3d)
    iend = 0
    n = length(x)
    while iend < n - 1
        # set istart to the first index that is finite
        istart = -1
        for j in (iend + 1):n
            if isfinite(x[j]) && isfinite(y[j]) && isfinite(z[j])
                istart = j
                break
            end
        end

        if istart > 0
            # iend is the last finite index
            iend = -1
            for j in (istart + 1):n
                if isfinite(x[j]) && isfinite(y[j]) && isfinite(z[j])
                    iend = j
                else
                    break
                end
            end
        end

        # if we found a start and end, draw the line segment, otherwise we're done
        if istart > 0 && iend > 0
            func(x[istart:iend], y[istart:iend], z[istart:iend])
        else
            break
        end
    end
end

gr_inqtext(x, y, s) = gr_inqtext(x, y, string(s))

function gr_inqtext(x, y, s::AbstractString)
    if length(s) >= 2 && s[1] == '$' && s[end] == '$'
        GR.inqmathtex(x, y, s[2:(end - 1)])
    elseif occursin('\\', s) || occursin("10^{", s)
        GR.inqtextext(x, y, s)
    else
        GR.inqtext(x, y, s)
    end
end

gr_text(x, y, s) = gr_text(x, y, string(s))

function gr_text(x, y, s::AbstractString)
    if length(s) >= 2 && s[1] == '$' && s[end] == '$'
        GR.mathtex(x, y, s[2:(end - 1)])
    elseif occursin('\\', s) || occursin("10^{", s)
        GR.textext(x, y, s)
    else
        GR.text(x, y, s)
    end
end

function gr_polaraxes(rmin::Real, rmax::Real, sp::Subplot)
    GR.savestate()
    xaxis = sp[:xaxis]
    yaxis = sp[:yaxis]

    α = 0:45:315
    a = α .+ 90
    sinf = sind.(a)
    cosf = cosd.(a)
    rtick_values, rtick_labels = get_ticks(sp, yaxis, update = false)

    #draw angular grid
    if xaxis[:grid]
        gr_set_line(
            xaxis[:gridlinewidth],
            xaxis[:gridstyle],
            xaxis[:foreground_color_grid],
            sp,
        )
        gr_set_transparency(xaxis[:foreground_color_grid], xaxis[:gridalpha])
        for i in eachindex(α)
            GR.polyline([sinf[i], 0], [cosf[i], 0])
        end
    end

    #draw radial grid
    if yaxis[:grid]
        gr_set_line(
            yaxis[:gridlinewidth],
            yaxis[:gridstyle],
            yaxis[:foreground_color_grid],
            sp,
        )
        gr_set_transparency(yaxis[:foreground_color_grid], yaxis[:gridalpha])
        for i in eachindex(rtick_values)
            r = (rtick_values[i] - rmin) / (rmax - rmin)
            if r <= 1.0 && r >= 0.0
                GR.drawarc(-r, r, -r, r, 0, 359)
            end
        end
        GR.drawarc(-1, 1, -1, 1, 0, 359)
    end

    #prepare to draw ticks
    gr_set_transparency(1)
    GR.setlinecolorind(90)
    GR.settextalign(GR.TEXT_HALIGN_CENTER, GR.TEXT_VALIGN_HALF)

    #draw angular ticks
    if xaxis[:showaxis]
        GR.drawarc(-1, 1, -1, 1, 0, 359)
        for i in eachindex(α)
            x, y = GR.wctondc(1.1 * sinf[i], 1.1 * cosf[i])
            GR.textext(x, y, string((360 - α[i]) % 360, "^o"))
        end
    end

    #draw radial ticks
    if yaxis[:showaxis]
        for i in eachindex(rtick_values)
            r = (rtick_values[i] - rmin) / (rmax - rmin)
            if r <= 1.0 && r >= 0.0
                x, y = GR.wctondc(0.05, r)
                gr_text(x, y, _cycle(rtick_labels, i))
            end
        end
    end
    GR.restorestate()
end

# using the axis extrema and limit overrides, return the min/max value for this axis
gr_x_axislims(sp::Subplot) = axis_limits(sp, :x)
gr_y_axislims(sp::Subplot) = axis_limits(sp, :y)
gr_z_axislims(sp::Subplot) = axis_limits(sp, :z)
gr_xy_axislims(sp::Subplot) = gr_x_axislims(sp)..., gr_y_axislims(sp)...

function gr_fill_viewport(vp::AVec{Float64}, c)
    GR.savestate()
    GR.selntran(0)
    GR.setscale(0)
    GR.setfillintstyle(GR.INTSTYLE_SOLID)
    gr_set_fillcolor(c)
    GR.fillrect(vp...)
    GR.selntran(1)
    GR.restorestate()
end

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

# draw ONE Shape
function gr_draw_marker(series, xi, yi, clims, i, msize, strokewidth, shape::Shape)
    sx, sy = coords(shape)
    # convert to ndc coords (percentages of window) ...
    w, h = get_size(series)
    f = msize / (w + h)

    xi, yi = GR.wctondc(xi, yi)

    # ... convert back to world coordinates
    xs_ys = GR.ndctowc.(xi .+ sx .* f, yi .+ sy .* f)
    xs, ys = getindex.(xs_ys, 1), getindex.(xs_ys, 2)

    # draw the interior
    mc = get_markercolor(series, clims, i)
    gr_set_fill(mc)
    gr_set_transparency(mc, get_markeralpha(series, i))
    GR.fillarea(xs, ys)

    # draw the shapes
    msc = get_markerstrokecolor(series, i)
    gr_set_line(strokewidth, :solid, msc, series)
    gr_set_transparency(msc, get_markerstrokealpha(series, i))
    GR.polyline(xs, ys)
end

function gr_nominal_size(s)
    w, h = get_size(s)
    min(w, h) / 500
end

# draw ONE symbol marker
function gr_draw_marker(series, xi, yi, clims, i, msize, strokewidth, shape::Symbol)
    GR.setborderwidth(strokewidth)
    gr_set_bordercolor(get_markerstrokecolor(series, i))
    gr_set_markercolor(get_markercolor(series, clims, i))
    gr_set_transparency(get_markeralpha(series, i))
    GR.setmarkertype(gr_markertype(shape))
    GR.setmarkersize(0.3msize / gr_nominal_size(series))
    GR.polymarker([xi], [yi])
end

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

function gr_set_line(lw, style, c, s) # s can be Subplot or Series
    GR.setlinetype(gr_linetype(style))
    GR.setlinewidth(get_thickness_scaling(s) * max(0, lw / gr_nominal_size(s)))
    gr_set_linecolor(c)
end

function gr_set_fill(c) #, a)
    gr_set_fillcolor(c) #, a)
    GR.setfillintstyle(GR.INTSTYLE_SOLID)
end

# this stores the conversion from a font pointsize to "percentage of window height"
# (which is what GR uses). `s` can be a Series, Subplot or Plot
gr_point_mult(s) = 1.5 * get_thickness_scaling(s) * px / pt / maximum(get_size(s))

# set the font attributes.
function gr_set_font(
    f::Font,
    s;
    halign = f.halign,
    valign = f.valign,
    color = f.color,
    rotation = f.rotation,
)
    family = lowercase(f.family)
    GR.setcharheight(gr_point_mult(s) * f.pointsize)
    GR.setcharup(sind(-rotation), cosd(-rotation))
    if haskey(gr_font_family, family)
        GR.settextfontprec(
            gr_font_family[family],
            gr_font_family[family] >= 200 ? 3 : GR.TEXT_PRECISION_STRING,
        )
    end
    gr_set_textcolor(color)
    GR.settextalign(gr_halign(halign), gr_valign(valign))
end

function gr_w3tondc(x, y, z)
    xw, yw, zw = GR.wc3towc(x, y, z)
    x, y = GR.wctondc(xw, yw)
    return x, y
end

# --------------------------------------------------------------------------------------
# viewport plot area
function gr_viewport_from_bbox(
    sp::Subplot{GRBackend},
    bb::BoundingBox,
    w,
    h,
    viewport_canvas,
)
    viewport = zeros(4)
    viewport[1] = viewport_canvas[2] * (left(bb) / w)
    viewport[2] = viewport_canvas[2] * (right(bb) / w)
    viewport[3] = viewport_canvas[4] * (1.0 - bottom(bb) / h)
    viewport[4] = viewport_canvas[4] * (1.0 - top(bb) / h)
    if hascolorbar(sp)
        viewport[2] -= 0.1 * (1 + RecipesPipeline.is3d(sp) / 2)
    end
    viewport
end

# change so we're focused on the viewport area
function gr_set_viewport_cmap(sp::Subplot, viewport_plotarea)
    GR.setviewport(
        viewport_plotarea[2] + (RecipesPipeline.is3d(sp) ? 0.07 : 0.02),
        viewport_plotarea[2] + (RecipesPipeline.is3d(sp) ? 0.10 : 0.05),
        viewport_plotarea[3],
        viewport_plotarea[4],
    )
end

function gr_set_viewport_polar(viewport_plotarea)
    xmin, xmax, ymin, ymax = viewport_plotarea
    ymax -= 0.05 * (xmax - xmin)
    xcenter = 0.5 * (xmin + xmax)
    ycenter = 0.5 * (ymin + ymax)
    r = 0.5 * NaNMath.min(xmax - xmin, ymax - ymin)
    GR.setviewport(xcenter - r, xcenter + r, ycenter - r, ycenter + r)
    GR.setwindow(-1, 1, -1, 1)
    r
end

struct GRColorbar
    gradients
    fills
    lines
    GRColorbar() = new([], [], [])
end

function gr_update_colorbar!(cbar::GRColorbar, series::Series)
    style = colorbar_style(series)
    style === nothing && return
    list =
        style == cbar_gradient ? cbar.gradients :
        style == cbar_fill ? cbar.fills :
        style == cbar_lines ? cbar.lines : error("Unknown colorbar style: $style.")
    push!(list, series)
end

function gr_contour_levels(series::Series, clims)
    levels = contour_levels(series, clims)
    if isfilledcontour(series)
        # GR implicitly uses the maximal z value as the highest level
        levels = levels[1:(end - 1)]
    end
    levels
end

function gr_colorbar_colors(series::Series, clims)
    if iscontour(series)
        levels = gr_contour_levels(series, clims)
        if isfilledcontour(series)
            # GR.contourf uses a color range according to supplied levels
            zrange = ignorenan_extrema(levels)
        else
            # GR.contour uses a color range according to data range
            zrange = clims
        end
        colors = 1000 .+ 255 .* (levels .- zrange[1]) ./ (zrange[2] - zrange[1])
    else
        colors = 1000:1255
    end
    round.(Int, colors)
end

function _cbar_unique(values, propname)
    out = last(values)
    if any(x != out for x in values)
        @warn "Multiple series with different $propname share a colorbar. " *
              "Colorbar may not reflect all series correctly."
    end
    out
end

# add the colorbar
function gr_draw_colorbar(cbar::GRColorbar, sp::Subplot, clims, viewport_plotarea)
    GR.savestate()
    xmin, xmax = gr_xy_axislims(sp)[1:2]
    zmin, zmax = clims[1:2]
    gr_set_viewport_cmap(sp, viewport_plotarea)
    GR.setscale(0)
    GR.setwindow(xmin, xmax, zmin, zmax)
    if !isempty(cbar.gradients)
        series = cbar.gradients
        gr_set_gradient(_cbar_unique(get_colorgradient.(series), "color"))
        gr_set_transparency(_cbar_unique(get_fillalpha.(series), "fill alpha"))
        GR.cellarray(xmin, xmax, zmax, zmin, 1, 256, 1000:1255)
    end

    if !isempty(cbar.fills)
        series = cbar.fills
        GR.setfillintstyle(GR.INTSTYLE_SOLID)
        gr_set_gradient(_cbar_unique(get_colorgradient.(series), "color"))
        gr_set_transparency(_cbar_unique(get_fillalpha.(series), "fill alpha"))
        levels = _cbar_unique(contour_levels.(series, Ref(clims)), "levels")
        # GR implicitly uses the maximal z value as the highest level
        if levels[end] < clims[2]
            @warn "GR: highest contour level less than maximal z value is not supported."
            # replace levels, rather than assign to levels[end], to ensure type
            # promotion in case levels is an integer array
            levels = [levels[1:(end - 1)]; clims[2]]
        end
        colors = gr_colorbar_colors(last(series), clims)
        for (from, to, color) in zip(levels[1:(end - 1)], levels[2:end], colors)
            GR.setfillcolorind(color)
            GR.fillrect(xmin, xmax, from, to)
        end
    end

    if !isempty(cbar.lines)
        series = cbar.lines
        gr_set_gradient(_cbar_unique(get_colorgradient.(series), "color"))
        gr_set_line(
            _cbar_unique(get_linewidth.(series), "line width"),
            _cbar_unique(get_linestyle.(series), "line style"),
            _cbar_unique(get_linecolor.(series, Ref(clims)), "line color"),
            sp,
        )
        gr_set_transparency(_cbar_unique(get_linealpha.(series), "line alpha"))
        levels = _cbar_unique(contour_levels.(series, Ref(clims)), "levels")
        colors = gr_colorbar_colors(last(series), clims)
        for (line, color) in zip(levels, colors)
            GR.setlinecolorind(color)
            GR.polyline([xmin, xmax], [line, line])
        end
    end

    ztick = 0.5 * GR.tick(zmin, zmax)
    gr_set_line(1, :solid, plot_color(:black), sp)
    GR.axes(0, ztick, xmax, zmin, 0, 1, 0.005)

    title = if isa(sp[:colorbar_title], PlotText)
        sp[:colorbar_title]
    else
        text(sp[:colorbar_title], colorbartitlefont(sp))
    end
    gr_set_font(title.font, sp)
    GR.settextalign(GR.TEXT_HALIGN_CENTER, GR.TEXT_VALIGN_TOP)
    GR.setcharup(-1, 0)
    gr_text(viewport_plotarea[2] + 0.1, gr_view_ycenter(viewport_plotarea), title.str)

    GR.restorestate()
end

gr_view_xcenter(viewport_plotarea) = 0.5 * (viewport_plotarea[1] + viewport_plotarea[2])
gr_view_ycenter(viewport_plotarea) = 0.5 * (viewport_plotarea[3] + viewport_plotarea[4])

gr_view_xposition(viewport_plotarea, position) =
    viewport_plotarea[1] + position * (viewport_plotarea[2] - viewport_plotarea[1])
gr_view_yposition(viewport_plotarea, position) =
    viewport_plotarea[3] + position * (viewport_plotarea[4] - viewport_plotarea[3])

function position(symb)
    if symb == :top || symb == :right
        return 0.95
    elseif symb == :left || symb == :bottom
        return 0.05
    end
    return 0.5
end

function alignment(symb)
    if symb == :top || symb == :right
        return GR.TEXT_HALIGN_RIGHT
    elseif symb == :left || symb == :bottom
        return GR.TEXT_HALIGN_LEFT
    end
    return GR.TEXT_HALIGN_CENTER
end

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

function gr_set_gradient(c)
    grad = _as_gradient(c)
    for (i, z) in enumerate(range(0, stop = 1, length = 256))
        c = grad[z]
        GR.setcolorrep(999 + i, red(c), green(c), blue(c))
    end
    grad
end

function gr_set_gradient(series::Series)
    color = get_colorgradient(series)
    color !== nothing && gr_set_gradient(color)
end

# this is our new display func... set up the viewport_canvas, compute bounding boxes, and display each subplot
function gr_display(plt::Plot, fmt = "")
    GR.clearws()

    dpi_factor = fmt == "png" ? plt[:dpi] / Plots.DPI : 1

    # collect some monitor/display sizes in meters and pixels
    display_width_meters, display_height_meters, display_width_px, display_height_px =
        GR.inqdspsize()
    display_width_ratio = display_width_meters / display_width_px
    display_height_ratio = display_height_meters / display_height_px

    # compute the viewport_canvas, normalized to the larger dimension
    viewport_canvas = Float64[0, 1, 0, 1]
    w, h = get_size(plt)
    if w > h
        ratio = float(h) / w
        msize = display_width_ratio * w * dpi_factor
        GR.setwsviewport(0, msize, 0, msize * ratio)
        GR.setwswindow(0, 1, 0, ratio)
        viewport_canvas[3] *= ratio
        viewport_canvas[4] *= ratio
    else
        ratio = float(w) / h
        msize = display_height_ratio * h * dpi_factor
        GR.setwsviewport(0, msize * ratio, 0, msize)
        GR.setwswindow(0, ratio, 0, 1)
        viewport_canvas[1] *= ratio
        viewport_canvas[2] *= ratio
    end

    # fill in the viewport_canvas background
    gr_fill_viewport(viewport_canvas, plt[:background_color_outside])

    # subplots:
    for sp in plt.subplots
        gr_display(sp, w * px, h * px, viewport_canvas)
    end

    GR.updatews()
end

function gr_set_tickfont(sp, letter)
    axis = sp[get_attr_symbol(letter, :axis)]

    # invalidate alignment changes for small rotations (|θ| < 45°)
    trigger(rot) = abs(sind(rot)) < abs(cosd(rot)) ? 0 : sign(rot)

    rot = axis[:rotation]
    if letter === :x || (RecipesPipeline.is3d(sp) && letter === :y)
        halign = (:left, :hcenter, :right)[trigger(rot) + 2]
        valign = (axis[:mirror] ? :bottom : :top, :vcenter)[trigger(abs(rot)) + 1]
    else
        halign = (axis[:mirror] ? :left : :right, :hcenter)[trigger(abs(rot)) + 1]
        valign = (:top, :vcenter, :bottom)[trigger(rot) + 2]
    end
    gr_set_font(
        tickfont(axis),
        sp,
        halign = halign,
        valign = valign,
        rotation = axis[:rotation],
        color = axis[:tickfontcolor],
    )
end

function gr_text_size(str)
    GR.savestate()
    GR.selntran(0)
    xs, ys = gr_inqtext(0, 0, string(str))
    l, r = extrema(xs)
    b, t = extrema(ys)
    w = r - l
    h = t - b
    GR.restorestate()
    return w, h
end

function gr_text_size(str, rot)
    GR.savestate()
    GR.selntran(0)
    xs, ys = gr_inqtext(0, 0, string(str))
    l, r = extrema(xs)
    b, t = extrema(ys)
    w = text_box_width(r - l, t - b, rot)
    h = text_box_height(r - l, t - b, rot)
    GR.restorestate()
    return w, h
end

text_box_width(w, h, rot) = abs(cosd(rot)) * w + abs(cosd(rot + 90)) * h
text_box_height(w, h, rot) = abs(sind(rot)) * w + abs(sind(rot + 90)) * h

function gr_get_ticks_size(ticks, rot)
    w, h = 0.0, 0.0
    for (cv, dv) in zip(ticks...)
        wi, hi = gr_text_size(dv, rot)
        w = NaNMath.max(w, wi)
        h = NaNMath.max(h, hi)
    end
    return w, h
end

function labelfunc(scale::Symbol, backend::GRBackend)
    texfunc = labelfunc_tex(scale)
    # replace dash with \minus (U+2212)
    label -> replace(texfunc(label), "-" => "−")
end

function gr_axis_height(sp, axis)
    GR.savestate()
    ticks = get_ticks(sp, axis, update = false)
    gr_set_font(tickfont(axis), sp)
    h = (
        ticks in (nothing, false, :none) ? 0 :
        last(gr_get_ticks_size(ticks, axis[:rotation]))
    )
    if axis[:guide] != ""
        gr_set_font(guidefont(axis), sp)
        h += last(gr_text_size(axis[:guide]))
    end
    GR.restorestate()
    return h
end

function gr_axis_width(sp, axis)
    GR.savestate()
    ticks = get_ticks(sp, axis, update = false)
    gr_set_font(tickfont(axis), sp)
    w = (
        ticks in (nothing, false, :none) ? 0 :
        first(gr_get_ticks_size(ticks, axis[:rotation]))
    )
    if axis[:guide] != ""
        gr_set_font(guidefont(axis), sp)
        w += last(gr_text_size(axis[:guide]))
    end
    GR.restorestate()
    return w
end

function _update_min_padding!(sp::Subplot{GRBackend})
    dpi = sp.plt[:thickness_scaling]
    ENV["GKS_ENCODING"] = "utf8"
    if !haskey(ENV, "GKSwstype")
        if isijulia()
            ENV["GKSwstype"] = "svg"
        end
    end
    # Add margin given by the user
    leftpad   = 2mm + sp[:left_margin]
    toppad    = 2mm + sp[:top_margin]
    rightpad  = 2mm + sp[:right_margin]
    bottompad = 2mm + sp[:bottom_margin]
    # Add margin for title
    if sp[:title] != ""
        gr_set_font(titlefont(sp), sp)
        l = last(last(gr_text_size(sp[:title])))
        h = 1mm + get_size(sp)[2] * l * px
        toppad += h
    end

    if RecipesPipeline.is3d(sp)
        xaxis, yaxis, zaxis = sp[:xaxis], sp[:yaxis], sp[:zaxis]
        xticks, yticks, zticks =
            get_ticks(sp, xaxis), get_ticks(sp, yaxis), get_ticks(sp, zaxis)
        # Add margin for x and y ticks
        h = 0mm
        if !isempty(first(xticks))
            gr_set_font(
                tickfont(xaxis),
                sp,
                halign = (:left, :hcenter, :right)[sign(xaxis[:rotation]) + 2],
                valign = (xaxis[:mirror] ? :bottom : :top),
                rotation = xaxis[:rotation],
            )
            l = 0.01 + last(gr_get_ticks_size(xticks, xaxis[:rotation]))
            h = max(h, 1mm + get_size(sp)[2] * l * px)
        end
        if !isempty(first(yticks))
            gr_set_font(
                tickfont(yaxis),
                sp,
                halign = (:left, :hcenter, :right)[sign(yaxis[:rotation]) + 2],
                valign = (yaxis[:mirror] ? :bottom : :top),
                rotation = yaxis[:rotation],
            )
            l = 0.01 + last(gr_get_ticks_size(yticks, yaxis[:rotation]))
            h = max(h, 1mm + get_size(sp)[2] * l * px)
        end
        if h > 0mm
            if xaxis[:mirror] || yaxis[:mirror]
                toppad += h
            end
            if !xaxis[:mirror] || !yaxis[:mirror]
                bottompad += h
            end
        end

        if !isempty(first(zticks))
            gr_set_font(
                tickfont(zaxis),
                sp,
                halign = (zaxis[:mirror] ? :left : :right),
                valign = (:top, :vcenter, :bottom)[sign(zaxis[:rotation]) + 2],
                rotation = zaxis[:rotation],
                color = zaxis[:tickfontcolor],
            )
            l = 0.01 + first(gr_get_ticks_size(zticks, zaxis[:rotation]))
            w = 1mm + get_size(sp)[1] * l * px
            if zaxis[:mirror]
                rightpad += w
            else
                leftpad += w
            end
        end

        # Add margin for x or y label
        h = 0mm
        if xaxis[:guide] != ""
            gr_set_font(guidefont(sp[:xaxis]), sp)
            l = last(gr_text_size(sp[:xaxis][:guide]))
            h = max(h, 1mm + get_size(sp)[2] * l * px)
        end
        if yaxis[:guide] != ""
            gr_set_font(guidefont(sp[:yaxis]), sp)
            l = last(gr_text_size(sp[:yaxis][:guide]))
            h = max(h, 1mm + get_size(sp)[2] * l * px)
        end
        if h > 0mm
            if (
                xaxis[:guide_position] == :top ||
                (xaxis[:guide_position] == :auto && xaxis[:mirror] == true)
            )
                toppad += h
            else
                bottompad += h
            end
        end
        # Add margin for z label
        if zaxis[:guide] != ""
            gr_set_font(guidefont(sp[:zaxis]), sp)
            l = last(gr_text_size(sp[:zaxis][:guide]))
            w = 1mm + get_size(sp)[2] * l * px
            if (
                zaxis[:guide_position] == :right ||
                (zaxis[:guide_position] == :auto && zaxis[:mirror] == true)
            )
                rightpad += w
            else
                leftpad += w
            end
        end
    else
        # Add margin for x and y ticks
        xticks, yticks = get_ticks(sp, sp[:xaxis]), get_ticks(sp, sp[:yaxis])
        if !isempty(first(xticks))
            gr_set_tickfont(sp, :x)
            l = 0.01 + last(gr_get_ticks_size(xticks, sp[:xaxis][:rotation]))
            h = 1mm + get_size(sp)[2] * l * px
            if sp[:xaxis][:mirror]
                toppad += h
            else
                bottompad += h
            end
        end
        if !isempty(first(yticks))
            gr_set_tickfont(sp, :y)
            l = 0.01 + first(gr_get_ticks_size(yticks, sp[:yaxis][:rotation]))
            w = 1mm + get_size(sp)[1] * l * px
            if sp[:yaxis][:mirror]
                rightpad += w
            else
                leftpad += w
            end
        end

        # Add margin for x label
        if sp[:xaxis][:guide] != ""
            gr_set_font(guidefont(sp[:xaxis]), sp)
            l = last(gr_text_size(sp[:xaxis][:guide]))
            h = 1mm + get_size(sp)[2] * l * px
            if (
                sp[:xaxis][:guide_position] == :top ||
                (sp[:xaxis][:guide_position] == :auto && sp[:xaxis][:mirror] == true)
            )
                toppad += h
            else
                bottompad += h
            end
        end
        # Add margin for y label
        if sp[:yaxis][:guide] != ""
            gr_set_font(guidefont(sp[:yaxis]), sp)
            l = last(gr_text_size(sp[:yaxis][:guide]))
            w = 1mm + get_size(sp)[2] * l * px
            if (
                sp[:yaxis][:guide_position] == :right ||
                (sp[:yaxis][:guide_position] == :auto && sp[:yaxis][:mirror] == true)
            )
                rightpad += w
            else
                leftpad += w
            end
        end
    end
    if sp[:colorbar_title] != ""
        rightpad += 4mm
    end
    sp.minpad = Tuple(dpi * [leftpad, toppad, rightpad, bottompad])
end

function is_equally_spaced(v)
    d = collect(v[2:end] .- v[1:(end - 1)])
    all(d .≈ d[1])
end

remap(x, lo, hi) = (x - lo) / (hi - lo)
function get_z_normalized(z, clims...)
    isnan(z) && return 256 / 255
    return remap(clamp(z, clims...), clims...)
end

function gr_clims(args...)
    lo, hi = get_clims(args...)
    if lo == hi
        if lo == 0
            hi = one(hi)
        elseif lo < 0
            hi = zero(hi)
        else
            lo = zero(lo)
        end
    end
    return lo, hi
end

function gr_display(sp::Subplot{GRBackend}, w, h, viewport_canvas)
    _update_min_padding!(sp)

    # the viewports for this subplot
    viewport_subplot = gr_viewport_from_bbox(sp, bbox(sp), w, h, viewport_canvas)
    viewport_plotarea = gr_viewport_from_bbox(sp, plotarea(sp), w, h, viewport_canvas)

    # update viewport_plotarea
    leg = gr_get_legend_geometry(viewport_plotarea, sp)
    gr_update_viewport_legend!(viewport_plotarea, sp, leg)
    gr_update_viewport_ratio!(viewport_plotarea, sp)

    # fill in the plot area background
    gr_fill_plotarea(sp, viewport_plotarea)

    # set our plot area view
    GR.setviewport(viewport_plotarea...)

    # set the scale flags and window
    gr_set_window(sp, viewport_plotarea)

    # draw the axes
    gr_draw_axes(sp, viewport_plotarea)
    gr_add_title(sp, viewport_plotarea, viewport_subplot)

    # this needs to be here to point the colormap to the right indices
    GR.setcolormap(1000 + GR.COLORMAP_COOLWARM)

    # init the colorbar
    cbar = GRColorbar()

    for series in series_list(sp)
        gr_add_series(sp, series)
        gr_update_colorbar!(cbar, series)
    end

    # draw the colorbar
    hascolorbar(sp) && gr_draw_colorbar(cbar, sp, gr_clims(sp), viewport_plotarea)

    # add the legend
    gr_add_legend(sp, leg, viewport_plotarea)

    # add annotations
    for ann in sp[:annotations]
        x, y, val = locate_annotation(sp, ann...)
        x, y = if RecipesPipeline.is3d(sp)
            gr_w3tondc(x, y, z)
        else
            GR.wctondc(x, y)
        end
        gr_set_font(val.font, sp)
        gr_text(x, y, val.str)
    end
end

## Legend

function gr_add_legend(sp, leg, viewport_plotarea)
    if sp[:legend] ∉ (:none, :inline)
        GR.savestate()
        GR.selntran(0)
        GR.setscale(0)
        gr_set_font(legendfont(sp), sp)
        if leg.w > 0
            xpos, ypos = gr_legend_pos(sp, leg, viewport_plotarea)
            GR.setfillintstyle(GR.INTSTYLE_SOLID)
            gr_set_fillcolor(sp[:background_color_legend])
            GR.fillrect(
                xpos - leg.leftw,
                xpos + leg.textw + leg.rightw,
                ypos + leg.dy,
                ypos - leg.h,
            ) # Allocating white space for actual legend width here
            gr_set_line(1, :solid, sp[:foreground_color_legend], sp)
            GR.drawrect(
                xpos - leg.leftw,
                xpos + leg.textw + leg.rightw,
                ypos + leg.dy,
                ypos - leg.h,
            ) # Drawing actual legend width here
            i = 0
            if sp[:legendtitle] !== nothing
                GR.settextalign(GR.TEXT_HALIGN_CENTER, GR.TEXT_VALIGN_HALF)
                gr_set_font(legendtitlefont(sp), sp)
                gr_text(xpos - 0.03 + 0.5 * leg.w, ypos, string(sp[:legendtitle]))
                ypos -= leg.dy
                gr_set_font(legendfont(sp), sp)
            end
            for series in series_list(sp)
                clims = gr_clims(sp, series)
                should_add_to_legend(series) || continue
                st = series[:seriestype]
                lc = get_linecolor(series, clims)
                gr_set_line(sp[:legendfontsize] / 8, get_linestyle(series), lc, sp)

                if (
                    (st == :shape || series[:fillrange] !== nothing) &&
                    series[:ribbon] === nothing
                )
                    fc = get_fillcolor(series, clims)
                    gr_set_fill(fc)
                    fs = get_fillstyle(series, i)
                    gr_set_fillstyle(fs)
                    l, r = xpos - leg.width_factor * 3.5, xpos - leg.width_factor / 2
                    b, t = ypos - 0.4 * leg.dy, ypos + 0.4 * leg.dy
                    x = [l, r, r, l, l]
                    y = [b, b, t, t, b]
                    gr_set_transparency(fc, get_fillalpha(series))
                    gr_polyline(x, y, GR.fillarea)
                    lc = get_linecolor(series, clims)
                    gr_set_transparency(lc, get_linealpha(series))
                    gr_set_line(get_linewidth(series), get_linestyle(series), lc, sp)
                    st == :shape && gr_polyline(x, y)
                end

                if st in (:path, :straightline, :path3d)
                    gr_set_transparency(lc, get_linealpha(series))
                    if series[:fillrange] === nothing || series[:ribbon] !== nothing
                        GR.polyline(
                            [xpos - leg.width_factor * 3.5, xpos - leg.width_factor / 2],
                            [ypos, ypos],
                        )
                    else
                        GR.polyline(
                            [xpos - leg.width_factor * 3.5, xpos - leg.width_factor / 2],
                            [ypos + 0.4 * leg.dy, ypos + 0.4 * leg.dy],
                        )
                    end
                end

                if series[:markershape] != :none
                    ms = first(series[:markersize])
                    msw = first(series[:markerstrokewidth])
                    s, sw = if ms > 0
                        0.8 * sp[:legendfontsize], 0.8 * sp[:legendfontsize] * msw / ms
                    else
                        0, 0.8 * sp[:legendfontsize] * msw / 8
                    end
                    gr_draw_markers(series, xpos - leg.width_factor * 2, ypos, clims, s, sw)
                end

                lab = series[:label]
                GR.settextalign(GR.TEXT_HALIGN_LEFT, GR.TEXT_VALIGN_HALF)
                gr_set_textcolor(plot_color(sp[:legendfontcolor]))
                gr_text(xpos, ypos, string(lab))
                ypos -= leg.dy
            end
        end
        GR.selntran(1)
        GR.restorestate()
    end
end

function gr_legend_pos(sp::Subplot, leg, viewport_plotarea)
    s = sp[:legend]
    s isa Real && return gr_legend_pos(s, leg, viewport_plotarea)
    if s isa Tuple{<:Real,Symbol}
        if s[2] !== :outer
            return gr_legend_pos(s[1], leg, viewport_plotarea)
        end

        xaxis, yaxis = sp[:xaxis], sp[:yaxis]
        xmirror =
            xaxis[:guide_position] == :top ||
            (xaxis[:guide_position] == :auto && xaxis[:mirror] == true)
        ymirror =
            yaxis[:guide_position] == :right ||
            (yaxis[:guide_position] == :auto && yaxis[:mirror] == true)
        axisclearance = [
            !ymirror * gr_axis_width(sp, sp[:yaxis]),
            ymirror * gr_axis_width(sp, sp[:yaxis]),
            !xmirror * gr_axis_height(sp, sp[:xaxis]),
            xmirror * gr_axis_height(sp, sp[:xaxis]),
        ]
        return gr_legend_pos(s[1], leg, viewport_plotarea; axisclearance)
    end
    s isa Symbol || return gr_legend_pos(s, viewport_plotarea)
    str = string(s)
    if str == "best"
        str = "topright"
    end
    if occursin("outer", str)
        xaxis, yaxis = sp[:xaxis], sp[:yaxis]
        xmirror =
            xaxis[:guide_position] == :top ||
            (xaxis[:guide_position] == :auto && xaxis[:mirror] == true)
        ymirror =
            yaxis[:guide_position] == :right ||
            (yaxis[:guide_position] == :auto && yaxis[:mirror] == true)
    end
    if occursin("right", str)
        if occursin("outer", str)
            # As per https://github.com/jheinen/GR.jl/blob/master/src/jlgr.jl#L525
            xpos =
                viewport_plotarea[2] +
                leg.xoffset +
                leg.leftw +
                ymirror * gr_axis_width(sp, sp[:yaxis])
        else
            xpos = viewport_plotarea[2] - leg.rightw - leg.textw - leg.xoffset
        end
    elseif occursin("left", str)
        if occursin("outer", str)
            xpos =
                viewport_plotarea[1] - !ymirror * gr_axis_width(sp, sp[:yaxis]) -
                leg.xoffset * 2 - leg.rightw - leg.textw
        else
            xpos = viewport_plotarea[1] + leg.leftw + leg.xoffset
        end
    else
        xpos =
            (viewport_plotarea[2] - viewport_plotarea[1]) / 2 +
            viewport_plotarea[1] +
            leg.leftw - leg.rightw - leg.textw - leg.xoffset * 2
    end
    if occursin("top", str)
        if s == :outertop
            ypos =
                viewport_plotarea[4] +
                leg.yoffset +
                leg.h +
                xmirror * gr_axis_height(sp, sp[:xaxis])
        else
            ypos = viewport_plotarea[4] - leg.yoffset - leg.dy
        end
    elseif occursin("bottom", str)
        if s == :outerbottom
            ypos =
                viewport_plotarea[3] - leg.yoffset - leg.dy -
                !xmirror * gr_axis_height(sp, sp[:xaxis])
        else
            ypos = viewport_plotarea[3] + leg.yoffset + leg.h
        end
    else
        # Adding min y to shift legend pos to correct graph (#2377)
        ypos =
            (viewport_plotarea[4] - viewport_plotarea[3] + leg.h) / 2 + viewport_plotarea[3]
    end
    return xpos, ypos
end

function gr_legend_pos(v::Tuple{S,T}, viewport_plotarea) where {S<:Real,T<:Real}
    xpos = v[1] * (viewport_plotarea[2] - viewport_plotarea[1]) + viewport_plotarea[1]
    ypos = v[2] * (viewport_plotarea[4] - viewport_plotarea[3]) + viewport_plotarea[3]
    (xpos, ypos)
end

function gr_legend_pos(theta::Real, leg, viewport_plotarea; axisclearance = nothing)
    xcenter = +(viewport_plotarea[1:2]...) / 2
    ycenter = +(viewport_plotarea[3:4]...) / 2

    if isnothing(axisclearance)
        # Inner
        # rectangle where the anchor can legally be
        xmin = viewport_plotarea[1] + leg.xoffset + leg.leftw
        xmax = viewport_plotarea[2] - leg.xoffset - leg.rightw - leg.textw
        ymin = viewport_plotarea[3] + leg.yoffset + leg.h
        ymax = viewport_plotarea[4] - leg.yoffset - leg.dy
    else
        # Outer
        xmin =
            viewport_plotarea[1] - leg.xoffset - leg.rightw - leg.textw - axisclearance[1]
        xmax = viewport_plotarea[2] + leg.xoffset + leg.leftw + axisclearance[2]
        ymin = viewport_plotarea[3] - leg.yoffset - leg.dy - axisclearance[3]
        ymax = viewport_plotarea[4] + leg.yoffset + leg.h + axisclearance[4]
    end
    return legend_pos_from_angle(theta, xmin, xcenter, xmax, ymin, ycenter, ymax)
end

function gr_get_legend_geometry(viewport_plotarea, sp)
    legendn = legendw = dy = 0
    if sp[:legend] != :none
        GR.savestate()
        GR.selntran(0)
        GR.setscale(0)
        if sp[:legendtitle] !== nothing
            gr_set_font(legendtitlefont(sp), sp)
            legendn += 1
            tbx, tby = gr_inqtext(0, 0, string(sp[:legendtitle]))
            legendw = tbx[3] - tbx[1]
            dy = tby[3] - tby[1]
        end
        gr_set_font(legendfont(sp), sp)
        for series in series_list(sp)
            should_add_to_legend(series) || continue
            legendn += 1
            tbx, tby = gr_inqtext(0, 0, string(series[:label]))
            legendw = max(legendw, tbx[3] - tbx[1]) # Holds text width right now
            dy = max(dy, tby[3] - tby[1])
        end

        GR.setscale(1)
        GR.selntran(1)
        GR.restorestate()
    end

    legend_width_factor = (viewport_plotarea[2] - viewport_plotarea[1]) / 45 # Determines the width of legend box
    legend_textw = legendw
    legend_rightw = legend_width_factor
    legend_leftw = legend_width_factor * 4
    total_legendw = legend_textw + legend_leftw + legend_rightw

    x_legend_offset = (viewport_plotarea[2] - viewport_plotarea[1]) / 30
    y_legend_offset = (viewport_plotarea[4] - viewport_plotarea[3]) / 30

    dy *= get(sp[:extra_kwargs], :legend_hfactor, 1)

    legendh = dy * legendn

    return (
        w = legendw,
        h = legendh,
        dy = dy,
        leftw = legend_leftw,
        textw = legend_textw,
        rightw = legend_rightw,
        xoffset = x_legend_offset,
        yoffset = y_legend_offset,
        width_factor = legend_width_factor,
    )
end

## Viewport, window and scale

function gr_update_viewport_legend!(viewport_plotarea, sp, leg)
    s = sp[:legend]

    xaxis, yaxis = sp[:xaxis], sp[:yaxis]
    xmirror =
        xaxis[:guide_position] == :top ||
        (xaxis[:guide_position] == :auto && xaxis[:mirror] == true)
    ymirror =
        yaxis[:guide_position] == :right ||
        (yaxis[:guide_position] == :auto && yaxis[:mirror] == true)

    if s isa Tuple{<:Real,Symbol}
        if s[2] === :outer
            (x, y) = gr_legend_pos(sp, leg, viewport_plotarea) # Dry run, to figure out
            if x < viewport_plotarea[1]
                viewport_plotarea[1] +=
                    leg.leftw +
                    leg.textw +
                    leg.rightw +
                    leg.xoffset +
                    !ymirror * gr_axis_width(sp, sp[:yaxis])
            elseif x > viewport_plotarea[2]
                viewport_plotarea[2] -= leg.leftw + leg.textw + leg.rightw + leg.xoffset
            end
            if y < viewport_plotarea[3]
                viewport_plotarea[3] +=
                    leg.h + leg.dy + leg.yoffset + !xmirror * gr_axis_height(sp, sp[:xaxis])
            elseif y > viewport_plotarea[4]
                viewport_plotarea[4] -= leg.h + leg.dy + leg.yoffset
            end
        end
    end
    leg_str = string(s)
    if occursin("outer", leg_str)
        if occursin("right", leg_str)
            viewport_plotarea[2] -= leg.leftw + leg.textw + leg.rightw + leg.xoffset
        elseif occursin("left", leg_str)
            viewport_plotarea[1] +=
                leg.leftw +
                leg.textw +
                leg.rightw +
                leg.xoffset +
                !ymirror * gr_axis_width(sp, sp[:yaxis])
        elseif occursin("top", leg_str)
            viewport_plotarea[4] -= leg.h + leg.dy + leg.yoffset
        elseif occursin("bottom", leg_str)
            viewport_plotarea[3] +=
                leg.h + leg.dy + leg.yoffset + !xmirror * gr_axis_height(sp, sp[:xaxis])
        end
    end
    if s === :inline
        if sp[:yaxis][:mirror]
            viewport_plotarea[1] += leg.w
        else
            viewport_plotarea[2] -= leg.w
        end
    end
end

function gr_update_viewport_ratio!(viewport_plotarea, sp)
    ratio = get_aspect_ratio(sp)
    if ratio != :none
        xmin, xmax, ymin, ymax = gr_xy_axislims(sp)
        if ratio == :equal
            ratio = 1
        end
        viewport_ratio =
            (viewport_plotarea[2] - viewport_plotarea[1]) /
            (viewport_plotarea[4] - viewport_plotarea[3])
        window_ratio = (xmax - xmin) / (ymax - ymin) / ratio
        if window_ratio < viewport_ratio
            viewport_center = 0.5 * (viewport_plotarea[1] + viewport_plotarea[2])
            viewport_size =
                (viewport_plotarea[2] - viewport_plotarea[1]) * window_ratio /
                viewport_ratio
            viewport_plotarea[1] = viewport_center - 0.5 * viewport_size
            viewport_plotarea[2] = viewport_center + 0.5 * viewport_size
        elseif window_ratio > viewport_ratio
            viewport_center = 0.5 * (viewport_plotarea[3] + viewport_plotarea[4])
            viewport_size =
                (viewport_plotarea[4] - viewport_plotarea[3]) * viewport_ratio /
                window_ratio
            viewport_plotarea[3] = viewport_center - 0.5 * viewport_size
            viewport_plotarea[4] = viewport_center + 0.5 * viewport_size
        end
    end
end

function gr_set_window(sp, viewport_plotarea)
    if ispolar(sp)
        gr_set_viewport_polar(viewport_plotarea)
    else
        xmin, xmax, ymin, ymax = gr_xy_axislims(sp)
        needs_3d = needs_any_3d_axes(sp)
        if needs_3d
            zmin, zmax = gr_z_axislims(sp)
            zok = zmax > zmin
        else
            zok = true
        end

        scaleop = 0
        if xmax > xmin && ymax > ymin && zok
            sp[:xaxis][:scale] == :log10 && (scaleop |= GR.OPTION_X_LOG)
            sp[:yaxis][:scale] == :log10 && (scaleop |= GR.OPTION_Y_LOG)
            needs_3d && sp[:zaxis][:scale] == :log10 && (scaleop |= GR.OPTION_Z_LOG)
            sp[:xaxis][:flip] && (scaleop |= GR.OPTION_FLIP_X)
            sp[:yaxis][:flip] && (scaleop |= GR.OPTION_FLIP_Y)
            needs_3d && sp[:zaxis][:flip] && (scaleop |= GR.OPTION_FLIP_Z)
            # NOTE: setwindow sets the "data coordinate" limits of the current "viewport"
            GR.setwindow(xmin, xmax, ymin, ymax)
            GR.setscale(scaleop)
        end
    end
end

function gr_fill_plotarea(sp, viewport_plotarea)
    if !RecipesPipeline.is3d(sp)
        gr_fill_viewport(viewport_plotarea, plot_color(sp[:background_color_inside]))
    end
end

## Axes

function gr_draw_axes(sp, viewport_plotarea)
    GR.setlinewidth(sp.plt[:thickness_scaling])

    if RecipesPipeline.is3d(sp)
        # set space
        xmin, xmax, ymin, ymax = gr_xy_axislims(sp)
        zmin, zmax = gr_z_axislims(sp)

        camera = round.(Int, sp[:camera])

        warn_invalid(val) =
            if val < 0 || val > 90
                @warn "camera: $(val)° ∉ [0°, 90°]"
            end
        warn_invalid.(camera)

        GR.setspace(zmin, zmax, camera...)

        # fill the plot area
        gr_set_fill(plot_color(sp[:background_color_inside]))
        plot_area_x = [xmin, xmin, xmin, xmax, xmax, xmax, xmin]
        plot_area_y = [ymin, ymin, ymax, ymax, ymax, ymin, ymin]
        plot_area_z = [zmin, zmax, zmax, zmax, zmin, zmin, zmin]
        x_bg, y_bg =
            RecipesPipeline.unzip(GR.wc3towc.(plot_area_x, plot_area_y, plot_area_z))
        GR.fillarea(x_bg, y_bg)

        for letter in (:x, :y, :z)
            gr_draw_axis_3d(sp, letter, viewport_plotarea)
        end
    elseif ispolar(sp)
        r = gr_set_viewport_polar(viewport_plotarea)
        #rmin, rmax = GR.adjustrange(ignorenan_minimum(r), ignorenan_maximum(r))
        rmin, rmax = axis_limits(sp, :y)
        gr_polaraxes(rmin, rmax, sp)
    elseif sp[:framestyle] != :none
        for letter in (:x, :y)
            gr_draw_axis(sp, letter, viewport_plotarea)
        end
    end
end

function gr_draw_axis(sp, letter, viewport_plotarea)
    ax = axis_drawing_info(sp, letter)
    axis = sp[get_attr_symbol(letter, :axis)]

    # draw segments
    gr_draw_grid(sp, axis, ax.grid_segments)
    gr_draw_minorgrid(sp, axis, ax.minorgrid_segments)
    gr_draw_spine(sp, axis, ax.segments)
    gr_draw_border(sp, axis, ax.border_segments)
    gr_draw_ticks(sp, axis, ax.tick_segments)

    # labels
    gr_label_ticks(sp, letter, ax.ticks)
    gr_label_axis(sp, letter, viewport_plotarea)
end

function gr_draw_axis_3d(sp, letter, viewport_plotarea)
    ax = axis_drawing_info_3d(sp, letter)
    axis = sp[get_attr_symbol(letter, :axis)]

    # draw segments
    gr_draw_grid(sp, axis, ax.grid_segments, gr_polyline3d)
    gr_draw_minorgrid(sp, axis, ax.minorgrid_segments, gr_polyline3d)
    gr_draw_spine(sp, axis, ax.segments, gr_polyline3d)
    gr_draw_border(sp, axis, ax.border_segments, gr_polyline3d)
    gr_draw_ticks(sp, axis, ax.tick_segments, gr_polyline3d)

    # labels
    GR.setscale(0)
    gr_label_ticks_3d(sp, letter, ax.ticks)
    gr_label_axis_3d(sp, letter)
    gr_set_window(sp, viewport_plotarea)
end

function gr_draw_grid(sp, axis, segments, func = gr_polyline)
    if axis[:grid]
        gr_set_line(
            axis[:gridlinewidth],
            axis[:gridstyle],
            axis[:foreground_color_grid],
            sp,
        )
        gr_set_transparency(axis[:foreground_color_grid], axis[:gridalpha])
        func(coords(segments)...)
    end
end

function gr_draw_minorgrid(sp, axis, segments, func = gr_polyline)
    if axis[:minorgrid]
        gr_set_line(
            axis[:minorgridlinewidth],
            axis[:minorgridstyle],
            axis[:foreground_color_minor_grid],
            sp,
        )
        gr_set_transparency(axis[:foreground_color_minor_grid], axis[:minorgridalpha])
        func(coords(segments)...)
    end
end

function gr_draw_spine(sp, axis, segments, func = gr_polyline)
    if axis[:showaxis]
        gr_set_line(1, :solid, axis[:foreground_color_border], sp)
        gr_set_transparency(1.0)
        GR.setclip(0)
        func(coords(segments)...)
        GR.setclip(1)
    end
end

function gr_draw_border(sp, axis, segments, func = gr_polyline)
    intensity = sp[:framestyle] == :semi ? 0.5 : 1
    if sp[:framestyle] in (:box, :semi)
        GR.setclip(0)
        gr_set_line(intensity, :solid, axis[:foreground_color_border], sp)
        gr_set_transparency(axis[:foreground_color_border], intensity)
        func(coords(segments)...)
        GR.setclip(1)
    end
end

function gr_draw_ticks(sp, axis, segments, func = gr_polyline)
    if axis[:showaxis]
        if sp[:framestyle] in (:zerolines, :grid)
            gr_set_line(1, :solid, axis[:foreground_color_grid], sp)
            gr_set_transparency(
                axis[:foreground_color_grid],
                axis[:tick_direction] == :out ? axis[:gridalpha] : 0,
            )
        else
            gr_set_line(1, :solid, axis[:foreground_color_axis], sp)
        end
        GR.setclip(0)
        func(coords(segments)...)
        GR.setclip(1)
    end
end

function gr_label_ticks(sp, letter, ticks)
    axis = sp[get_attr_symbol(letter, :axis)]
    isy = letter === :y
    oletter = isy ? :x : :y
    oaxis = sp[get_attr_symbol(oletter, :axis)]
    oamin, oamax = axis_limits(sp, oletter)
    gr_set_tickfont(sp, letter)
    out_factor = ifelse(axis[:tick_direction] === :out, 1.5, 1)
    x_offset = isy ? -1.5e-2 * out_factor : 0
    y_offset = isy ? 0 : -8e-3 * out_factor

    ov = sp[:framestyle] == :origin ? 0 : xor(oaxis[:flip], axis[:mirror]) ? oamax : oamin
    sgn = axis[:mirror] ? -1 : 1
    for (cv, dv) in zip(ticks...)
        x, y = GR.wctondc(reverse_if((cv, ov), isy)...)
        gr_text(x + sgn * x_offset, y + sgn * y_offset, dv)
    end
end

function gr_label_ticks(sp, letter, ticks::Nothing) end

function gr_label_ticks_3d(sp, letter, ticks)
    near_letter = letter in (:x, :z) ? :y : :x
    far_letter = letter in (:x, :y) ? :z : :x

    ax = sp[get_attr_symbol(letter, :axis)]
    nax = sp[get_attr_symbol(near_letter, :axis)]
    fax = sp[get_attr_symbol(far_letter, :axis)]

    amin, amax = axis_limits(sp, letter)
    namin, namax = axis_limits(sp, near_letter)
    famin, famax = axis_limits(sp, far_letter)
    n0, n1 = letter === :y ? (namax, namin) : (namin, namax)

    # find out which axes we are dealing with
    i = findfirst(==(letter), (:x, :y, :z))
    letters = axes_shift((:x, :y, :z), 1 - i)
    asyms = get_attr_symbol.(letters, :axis)

    # get axis objects, ticks and minor ticks
    # regardless of the `letter` we now use the convention that `x` in variable names refer to
    # the first axesm `y` to the second, etc ...
    ylims, zlims = axis_limits.(Ref(sp), letters[2:3])
    xax, yax, zax = getindex.(Ref(sp), asyms)

    gr_set_tickfont(sp, letter)
    nt = sp[:framestyle] == :origin ? 0 : ax[:mirror] ? n1 : n0
    ft = sp[:framestyle] == :origin ? 0 : ax[:mirror] ? famax : famin

    xoffset = if letter === :x
        (sp[:yaxis][:mirror] ? 1 : -1) * 1e-2 * (sp[:xaxis][:tick_direction] == :out ? 1.5 : 1)
    elseif letter === :y
        (sp[:yaxis][:mirror] ? -1 : 1) * 1e-2 * (sp[:yaxis][:tick_direction] == :out ? 1.5 : 1)
    else
        (sp[:zaxis][:mirror] ? 1 : -1) * 1e-2 * (sp[:zaxis][:tick_direction] == :out ? 1.5 : 1)
    end
    yoffset = if letter === :x
        (sp[:xaxis][:mirror] ? 1 : -1) * 1e-2 * (sp[:xaxis][:tick_direction] == :out ? 1.5 : 1)
    elseif letter === :y
        (sp[:yaxis][:mirror] ? 1 : -1) * 1e-2 * (sp[:yaxis][:tick_direction] == :out ? 1.5 : 1)
    else
        0
    end

    cvs, dvs = ticks
    ax[:flip] && reverse!(cvs)

    for (cv, dv) in zip((cvs, dvs)...)
        xi, yi = gr_w3tondc(sort_3d_axes(cv, nt, ft, letter)...)
        gr_text(xi + xoffset, yi + yoffset, dv)
    end
end

function gr_label_axis(sp, letter, viewport_plotarea)
    axis = sp[get_attr_symbol(letter, :axis)]
    mirror = axis[:mirror]
    # guide
    if axis[:guide] != ""
        GR.savestate()
        gr_set_font(guidefont(axis), sp)
        guide_position = axis[:guide_position]
        angle = float(axis[:guidefontrotation])  # github.com/JuliaPlots/Plots.jl/issues/3089
        if letter === :y
            angle += 180.0  # default angle = 0. should yield GR.setcharup(-1, 0) i.e. 180°
            GR.setcharup(cosd(angle), sind(angle))
            ypos = gr_view_yposition(viewport_plotarea, position(axis[:guidefontvalign]))
            yalign = alignment(axis[:guidefontvalign])
            if guide_position === :right || (guide_position == :auto && mirror)
                GR.settextalign(yalign, GR.TEXT_VALIGN_BOTTOM)
                xpos = viewport_plotarea[2] + 0.03 + mirror * gr_axis_width(sp, axis)
            else
                GR.settextalign(yalign, GR.TEXT_VALIGN_TOP)
                xpos = viewport_plotarea[1] - 0.03 - !mirror * gr_axis_width(sp, axis)
            end
        else
            angle += 90.0  # default angle = 0. should yield GR.setcharup(0, 1) i.e. 90°
            GR.setcharup(cosd(angle), sind(angle))
            xpos = gr_view_xposition(viewport_plotarea, position(axis[:guidefonthalign]))
            xalign = alignment(axis[:guidefonthalign])
            if guide_position === :top || (guide_position == :auto && mirror)
                GR.settextalign(xalign, GR.TEXT_VALIGN_TOP)
                ypos =
                    viewport_plotarea[4] +
                    0.015 +
                    (mirror ? gr_axis_height(sp, axis) : 0.015)
            else
                GR.settextalign(xalign, GR.TEXT_VALIGN_BOTTOM)
                ypos =
                    viewport_plotarea[3] - 0.015 -
                    (mirror ? 0.015 : gr_axis_height(sp, axis))
            end
        end
        gr_text(xpos, ypos, axis[:guide])
        GR.restorestate()
    end
end

function gr_label_axis_3d(sp, letter)
    ax = sp[get_attr_symbol(letter, :axis)]
    if ax[:guide] != ""
        near_letter = letter in (:x, :z) ? :y : :x
        far_letter = letter in (:x, :y) ? :z : :x

        nax = sp[get_attr_symbol(near_letter, :axis)]
        fax = sp[get_attr_symbol(far_letter, :axis)]

        amin, amax = axis_limits(sp, letter)
        namin, namax = axis_limits(sp, near_letter)
        famin, famax = axis_limits(sp, far_letter)
        n0, n1 = letter === :y ? (namax, namin) : (namin, namax)

        GR.savestate()
        gr_set_font(
            guidefont(ax),
            sp,
            halign = (:left, :hcenter, :right)[sign(ax[:rotation]) + 2],
            valign = ax[:mirror] ? :bottom : :top,
            rotation = ax[:rotation],
            # color = ax[:guidefontcolor],
        )
        ag = (amin + amax) / 2
        ng = ax[:mirror] ? n1 : n0
        fg = ax[:mirror] ? famax : famin
        x, y = gr_w3tondc(sort_3d_axes(ag, ng, fg, letter)...)
        if letter in (:x, :y)
            h = gr_axis_height(sp, ax)
            x_offset = letter === :x ? -h : h
            y_offset = -h
        else
            x_offset = -0.03 - gr_axis_width(sp, ax)
            y_offset = 0
        end
        letter === :z && GR.setcharup(-1, 0)
        sgn = ax[:mirror] ? -1 : 1
        gr_text(x + sgn * x_offset, y + sgn * y_offset, ax[:guide])
        GR.restorestate()
    end
end

function gr_add_title(sp, viewport_plotarea, viewport_subplot)
    if sp[:title] != ""
        GR.savestate()
        gr_set_font(titlefont(sp), sp)
        loc = sp[:titlelocation]
        if loc == :left
            xpos = viewport_plotarea[1]
            halign = GR.TEXT_HALIGN_LEFT
        elseif loc == :right
            xpos = viewport_plotarea[2]
            halign = GR.TEXT_HALIGN_RIGHT
        else
            xpos = gr_view_xcenter(viewport_plotarea)
            halign = GR.TEXT_HALIGN_CENTER
        end
        GR.settextalign(halign, GR.TEXT_VALIGN_TOP)
        gr_text(xpos, viewport_subplot[4], sp[:title])
        GR.restorestate()
    end
end

## Series

function gr_add_series(sp, series)
    st = series[:seriestype]

    # update the current stored gradient
    gr_set_gradient(series)

    GR.savestate()

    x, y, z = (handle_surface(series[letter]) for letter in (:x, :y, :z))
    xscale, yscale = sp[:xaxis][:scale], sp[:yaxis][:scale]
    frng = series[:fillrange]

    # recompute data
    if ispolar(sp) && z === nothing
        rmin, rmax = axis_limits(sp, :y)
        if frng !== nothing
            _, frng = convert_to_polar(x, frng, (rmin, rmax))
        end
        x, y = convert_to_polar(x, y, (rmin, rmax))
    end

    clims = gr_clims(sp, series)

    # add custom frame shapes to markershape?
    series_annotations_shapes!(series)
    # -------------------------------------------------------

    # draw the series
    if st in (:path, :scatter, :straightline)
        if st === :straightline
            x, y = straightline_data(series)
        end
        gr_draw_segments(series, x, y, frng, clims)
        if series[:markershape] !== :none
            gr_draw_markers(series, x, y, clims)
        end
    elseif st === :shape
        gr_draw_shapes(series, clims)
    elseif st in (:path3d, :scatter3d)
        gr_draw_segments_3d(series, x, y, z, clims)
        if st === :scatter3d || series[:markershape] !== :none
            # TODO: Do we need to transform to 2d coordinates here?
            x2, y2 = RecipesPipeline.unzip(map(GR.wc3towc, x, y, z))
            gr_draw_markers(series, x2, y2, clims)
        end
    elseif st === :contour
        gr_draw_contour(series, x, y, z, clims)
    elseif st in (:surface, :wireframe, :mesh3d)
        gr_draw_surface(series, x, y, z, clims)
    elseif st === :volume
        sp[:legend] = :none
        GR.gr3.clear()
        dmin, dmax = GR.gr3.volume(y.v, 0)
    elseif st === :heatmap
        # `z` is already transposed, so we need to reverse before passing its size.
        x, y = heatmap_edges(x, xscale, y, yscale, reverse(size(z)), ispolar(series))
        gr_draw_heatmap(series, x, y, z, clims)
    elseif st === :image
        gr_draw_image(series, x, y, z, clims)
    end

    # this is all we need to add the series_annotations text
    anns = series[:series_annotations]
    for (xi, yi, str, fnt) in EachAnn(anns, x, y)
        gr_set_font(fnt, sp)
        gr_text(GR.wctondc(xi, yi)..., str)
    end

    if sp[:legend] == :inline && should_add_to_legend(series)
        gr_set_font(legendfont(sp), sp)
        gr_set_textcolor(plot_color(sp[:legendfontcolor]))
        if sp[:yaxis][:mirror]
            (_, i) = sp[:xaxis][:flip] ? findmax(x) : findmin(x)
            GR.settextalign(GR.TEXT_HALIGN_RIGHT, GR.TEXT_VALIGN_HALF)
            offset = -0.01
        else
            (_, i) = sp[:xaxis][:flip] ? findmin(x) : findmax(x)
            GR.settextalign(GR.TEXT_HALIGN_LEFT, GR.TEXT_VALIGN_HALF)
            offset = 0.01
        end
        (x_l, y_l) = GR.wctondc(x[i], y[i])
        gr_text(x_l + offset, y_l, series[:label])
    end
    GR.restorestate()
end

function gr_draw_segments(series, x, y, fillrange, clims)
    st = series[:seriestype]
    if x !== nothing && length(x) > 1
        segments = series_segments(series, st; check = true)
        # do area fill
        if fillrange !== nothing
            GR.setfillintstyle(GR.INTSTYLE_SOLID)
            fr_from, fr_to = (is_2tuple(fillrange) ? fillrange : (y, fillrange))
            for segment in segments
                i, rng = segment.attr_index, segment.range
                fc = get_fillcolor(series, clims, i)
                gr_set_fillcolor(fc)
                fs = get_fillstyle(series, i)
                gr_set_fillstyle(fs)
                fx = _cycle(x, vcat(rng, reverse(rng)))
                fy = vcat(_cycle(fr_from, rng), _cycle(fr_to, reverse(rng)))
                gr_set_transparency(fc, get_fillalpha(series, i))
                GR.fillarea(fx, fy)
            end
        end

        # draw the line(s)
        if st in (:path, :straightline)
            for segment in segments
                i, rng = segment.attr_index, segment.range
                lc = get_linecolor(series, clims, i)
                gr_set_line(get_linewidth(series, i), get_linestyle(series, i), lc, series)
                arrowside = isa(series[:arrow], Arrow) ? series[:arrow].side : :none
                arrowstyle = isa(series[:arrow], Arrow) ? series[:arrow].style : :simple
                gr_set_fillcolor(lc)
                gr_set_transparency(lc, get_linealpha(series, i))
                gr_polyline(x[rng], y[rng]; arrowside = arrowside, arrowstyle = arrowstyle)
            end
        end
    end
end

function gr_draw_segments_3d(series, x, y, z, clims)
    if series[:seriestype] === :path3d && length(x) > 1
        lz = series[:line_z]
        segments = series_segments(series, :path3d; check = true)
        for segment in segments
            i, rng = segment.attr_index, segment.range
            lc = get_linecolor(series, clims, i)
            gr_set_line(get_linewidth(series, i), get_linestyle(series, i), lc, series)
            gr_set_transparency(lc, get_linealpha(series, i))
            GR.polyline3d(x[rng], y[rng], z[rng])
        end
    end
end

function gr_draw_markers(
    series::Series,
    x,
    y,
    clims,
    msize = series[:markersize],
    strokewidth = series[:markerstrokewidth],
)
    isempty(x) && return
    GR.setfillintstyle(GR.INTSTYLE_SOLID)

    shapes = series[:markershape]
    if shapes != :none
        for segment in series_segments(series, :scatter)
            i = segment.attr_index
            rng = intersect(eachindex(x), segment.range)
            if !isempty(rng)
                ms = get_thickness_scaling(series) * _cycle(msize, i)
                msw = get_thickness_scaling(series) * _cycle(strokewidth, i)
                shape = _cycle(shapes, i)
                for j in rng
                    gr_draw_marker(
                        series,
                        _cycle(x, j),
                        _cycle(y, j),
                        clims,
                        i,
                        ms,
                        msw,
                        shape,
                    )
                end
            end
        end
    end
end

function gr_draw_shapes(series, clims)
    x, y = shape_data(series)
    for segment in series_segments(series, :shape)
        i, rng = segment.attr_index, segment.range
        if length(rng) > 1
            # connect to the beginning
            rng = vcat(rng, rng[1])

            # get the segments
            xseg, yseg = x[rng], y[rng]

            # draw the interior
            fc = get_fillcolor(series, clims, i)
            gr_set_fill(fc)
            fs = get_fillstyle(series, i)
            gr_set_fillstyle(fs)
            gr_set_transparency(fc, get_fillalpha(series, i))
            GR.fillarea(xseg, yseg)

            # draw the shapes
            lc = get_linecolor(series, clims, i)
            gr_set_line(get_linewidth(series, i), get_linestyle(series, i), lc, series)
            gr_set_transparency(lc, get_linealpha(series, i))
            GR.polyline(xseg, yseg)
        end
    end
end

function gr_draw_contour(series, x, y, z, clims)
    GR.setspace(clims[1], clims[2], 0, 90)
    gr_set_line(get_linewidth(series), get_linestyle(series), get_linecolor(series), series)
    gr_set_transparency(get_fillalpha(series))
    is_lc_black = let black = plot_color(:black)
        plot_color(series[:linecolor]) in (black, [black])
    end
    h = gr_contour_levels(series, clims)
    if series[:fillrange] !== nothing
        if series[:fillcolor] != series[:linecolor] && !is_lc_black
            @warn "GR: filled contour only supported with black contour lines"
        end
        GR.contourf(x, y, h, z, series[:contour_labels] == true ? 1 : 0)
    else
        coff = is_lc_black ? 0 : 1000
        GR.contour(x, y, h, z, coff + (series[:contour_labels] == true ? 1 : 0))
    end
end

function gr_draw_surface(series, x, y, z, clims)
    e_kwargs = series[:extra_kwargs]
    st = series[:seriestype]
    if st === :surface
        if ndims(x) == ndims(y) == ndims(z) == 2
            GR.gr3.surface(x', y', z, GR.OPTION_3D_MESH)
        else
            fillalpha = get_fillalpha(series)
            fillcolor = get_fillcolor(series)
            # NOTE: setting nx = 0 or ny = 0 disables GR.gridit interpolation
            nx, ny = get(e_kwargs, :nx, 200), get(e_kwargs, :ny, 200)
            if length(x) == length(y) == length(z) && nx > 0 && ny > 0
                x, y, z = GR.gridit(x, y, z, nx, ny)
            end
            d_opt = get(e_kwargs, :display_option, GR.OPTION_COLORED_MESH)
            if (!isnothing(fillalpha) && fillalpha < 1) || alpha(first(fillcolor)) < 1
                gr_set_transparency(fillcolor, fillalpha)
                GR.surface(x, y, z, d_opt)
            else
                GR.gr3.surface(x, y, z, d_opt)
            end
        end
    elseif st === :wireframe
        GR.setfillcolorind(0)
        GR.surface(x, y, z, get(e_kwargs, :display_option, GR.OPTION_FILLED_MESH))
    elseif st === :mesh3d
        @warn "GR: mesh3d is experimental (no face colors)"
        gr_set_line(
            get_linewidth(series),
            get_linestyle(series),
            get_linecolor(series),
            series,
        )
        GR.polyline3d(mesh3d_triangles(x, y, z, series[:connections])...)
    else
        throw(ArgumentError("Not handled !"))
    end
end

function gr_draw_heatmap(series, x, y, z, clims)
    fillgrad = _as_gradient(series[:fillcolor])
    GR.setspace(clims..., 0, 90)
    w, h = length(x) - 1, length(y) - 1
    if !ispolar(series) && is_uniformly_spaced(x) && is_uniformly_spaced(y)
        # For uniformly spaced data use GR.drawimage, which can be
        # much faster than GR.nonuniformcellarray, especially for
        # pdf output, and also supports alpha values.
        # Note that drawimage draws uniformly spaced data correctly
        # even on log scales, where it is visually non-uniform.
        colors = plot_color.(get(fillgrad, z, clims), series[:fillalpha])
        rgba = gr_color.(colors)
        GR.drawimage(first(x), last(x), last(y), first(y), w, h, rgba)
    else
        if something(series[:fillalpha], 1) < 1
            @warn "GR: transparency not supported in non-uniform heatmaps. Alpha values ignored."
        end
        z_normalized = get_z_normalized.(z, clims...)
        rgba = Int32[round(Int32, 1000 + _i * 255) for _i in z_normalized]
        if !ispolar(series)
            GR.nonuniformcellarray(x, y, w, h, rgba)
        else
            if y[1] < 0
                @warn "'y[1] < 0' (rmin) is not yet supported."
            end
            xmin, xmax, ymin, ymax = gr_xy_axislims(series[:subplot])
            GR.setwindow(-ymax, ymax, -ymax, ymax)
            GR.nonuniformpolarcellarray(rad2deg.(x), y, w, h, rgba)
        end
    end
end

function gr_draw_image(series, x, y, z, clims)
    w, h = size(z)
    xmin, xmax = ignorenan_extrema(x)
    ymin, ymax = ignorenan_extrema(y)
    rgba = gr_color.(z)
    GR.drawimage(xmin, xmax, ymax, ymin, w, h, rgba)
end

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

for (mime, fmt) in (
    "application/pdf" => "pdf",
    "image/png" => "png",
    "application/postscript" => "ps",
    "image/svg+xml" => "svg",
)
    @eval function _show(io::IO, ::MIME{Symbol($mime)}, plt::Plot{GRBackend})
        ENV["GKS_ENCODING"] = "utf8"
        GR.emergencyclosegks()
        filepath = tempname() * "." * $fmt
        env = get(ENV, "GKSwstype", "0")
        ENV["GKSwstype"] = $fmt
        ENV["GKS_FILEPATH"] = filepath
        gr_display(plt, $fmt)
        GR.emergencyclosegks()
        write(io, read(filepath, String))
        rm(filepath)
        if env != "0"
            ENV["GKSwstype"] = env
        else
            pop!(ENV, "GKSwstype")
        end
    end
end

function _display(plt::Plot{GRBackend})
    ENV["GKS_ENCODING"] = "utf8"
    if plt[:display_type] == :inline
        GR.emergencyclosegks()
        filepath = tempname() * ".pdf"
        ENV["GKSwstype"] = "pdf"
        ENV["GKS_FILEPATH"] = filepath
        gr_display(plt)
        GR.emergencyclosegks()
        content = string(
            "\033]1337;File=inline=1;preserveAspectRatio=0:",
            base64encode(open(read, filepath)),
            "\a",
        )
        println(content)
        rm(filepath)
    else
        ENV["GKS_DOUBLE_BUF"] = true
        gr_display(plt)
    end
end

closeall(::GRBackend) = GR.emergencyclosegks()


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

using Markdown
using InteractiveUtils

# ╔═╡ 090f6f0f-fd67-419b-a35e-5392a59f39b1
using PlutoTest

# ╔═╡ e89c666e-b0e7-11eb-22cf-7f6fd1e299bb
md"""
# [Day 11: Corporate Policy](https://adventofcode.com/2015/day/11)

Santa's previous password expired, and he needs help choosing a new one.

To help him remember his new password after the old one expires, Santa has devised a method of coming up with a password based on the previous one. Corporate policy dictates that passwords must be exactly eight lowercase letters (for security reasons), so he finds his new password by incrementing his old password string repeatedly until it is valid.

Incrementing is just like counting with numbers: `xx`, `xy`, `xz`, `ya`, `yb`, and so on. Increase the rightmost letter one step; if it was `z`, it wraps around to `a`, and repeat with the next letter to the left until one doesn't wrap around.

Unfortunately for Santa, a new Security-Elf recently started, and he has imposed some additional password requirements:

-   Passwords must include one increasing straight of at least three letters, like `abc`, `bcd`, `cde`, and so on, up to `xyz`. They cannot skip letters; `abd` doesn't count.
-   Passwords may not contain the letters `i`, `o`, or `l`, as these letters can be mistaken for other characters and are therefore confusing.
-   Passwords must contain at least two different, non-overlapping pairs of letters, like `aa`, `bb`, or `zz`.

For example:

-   `hijklmmn` meets the first requirement (because it contains the straight `hij`) but fails the second requirement requirement (because it contains `i` and `l`).
-   `abbceffg` meets the third requirement (because it repeats `bb` and `ff`) but fails the first requirement.
-   `abbcegjk` fails the third requirement, because it only has one double letter (`bb`).
-   The next password after `abcdefgh` is `abcdffaa`.
-   The next password after `ghijklmn` is `ghjaabcc`, because you eventually skip all the passwords that start with `ghi...`, since `i` is not allowed.

**Given Santa's current password (your puzzle input), what should his next password be?**
"""

# ╔═╡ ec540096-d9d1-4b5f-b6a2-16c470322424
puzzle_input = "hxbxwxba"

# ╔═╡ 4e20e380-e928-4764-a982-3c0ca6bb942c
function next(str)
	chars = collect(str)
	vmin = Int('a')
	vmax = Int('z')

	for idx in reverse(1:length(str))
		v = Int(chars[idx]) + 1
		if v > vmax
			v -= vmax-vmin + 1
			chars[idx] = v
		else
			chars[idx] = v
			break
		end
	end

	join(chars)
end

# ╔═╡ 7a68015d-0282-4cb0-bebc-c265d5ca72d3
@assert next("xx") == "xy"

# ╔═╡ c4d6ad5b-f500-4435-a8eb-76c1f33f9775
@assert next("xy") == "xz"

# ╔═╡ 8f2717d2-ed5e-4e4e-abcd-6ff2982d8cd6
@assert next("xz") == "ya"

# ╔═╡ 58f76b47-542f-438b-ac71-694a76d9d77b
@assert next("ya") == "yb"

# ╔═╡ 9a03d3c5-3bc6-48c7-87cc-49ef6cbd31d4
function has_straight(s)
	for idx in 1:length(s)-2
		a = Int(s[idx])
		b = Int(s[idx+1])
		c = Int(s[idx+2])

		if b == a + 1 && c == b + 1
			return true
		end
	end
	false
end

# ╔═╡ b10b13e9-de70-4c09-8b7a-9cd530d1e3ac
@assert has_straight("hijklmmn")

# ╔═╡ 2d9e27ff-5557-4c7e-b6a8-b3accd15f7f0
@assert !has_straight("abbceffg")

# ╔═╡ 26342478-769e-4ef1-a3f3-ddf5aad2e281
has_no_confusing(s) = !occursin("i", s) && !occursin("o", s) && !occursin("l", s)

# ╔═╡ a470ee13-e98c-423b-8398-e6c8774983c5
@assert !has_no_confusing("hijklmmn")

# ╔═╡ dd4c2b0b-d8af-411f-9290-d36398772267
@assert has_no_confusing("abbceffg")

# ╔═╡ a0681275-076f-400f-8862-74c0cba5ffc1
# Passwords must contain at least two different, non-overlapping pairs of letters, like aa, bb, or zz.

# ╔═╡ 2144ecf3-79c2-4de1-8aa4-52b4e6b57208
function has_two_nonoverlapping_pair(s)
	for idx1 in 1:length(s)-1
		if s[idx1] == s[idx1+1]
			for idx2 in idx1+2:length(s)-1
				if s[idx2] == s[idx2+1]
					return true
				end
			end
		end
	end
	false
end

# ╔═╡ d103a1b1-7ade-41da-bc1e-b01337333675
@assert has_two_nonoverlapping_pair("abbceffg")

# ╔═╡ 5a13d4d9-7451-4fb0-8aaf-8d47947c10e4
check1(s) = has_straight(s) && has_no_confusing(s) && has_two_nonoverlapping_pair(s)

# ╔═╡ de4134e2-cc3f-4872-9244-4596d8525292
function next_valid(s)
	while true
		s = next(s)
		if check1(s)
			return s
		end
	end
end

# ╔═╡ 6b671ffc-ccc1-4b84-9ddb-7f7137995348
@assert next_valid("abcdefgh") == "abcdffaa"

# ╔═╡ a33d8ddb-5ef5-4c0c-891c-eef1368fbcc9
@assert next_valid("ghijklmn") == "ghjaabcc"

# ╔═╡ cbcd58ec-5f81-4ca4-8a4f-2a570db79f2a
part1 = next_valid(puzzle_input)

# ╔═╡ 2d3ab22d-11fe-4ee3-9ee1-9410ae03498a
@test part1 == "hxbxxyzz"

# ╔═╡ bd64f7b1-4074-4eeb-925d-eb8cc7a524bf
md"Your puzzle answer was `hxbxxyzz`."

# ╔═╡ aa05c55c-e79d-4846-95e1-823f0a5ea2e3
md"""
# Part Two
Santa's password expired again. **What's the next one?**
"""

# ╔═╡ 1ef77364-44a3-4819-b0c3-11cca842e53c
part2 = next_valid("hxbxxyzz")

# ╔═╡ fb47e75f-a3a0-4ee9-97c8-8a7e3d3f7831
@test part2 == "hxcaabcc"

# ╔═╡ bd67cd4e-af24-403c-a3ba-9753f11dc21c
md"Your puzzle answer was `hxcaabcc`."

# ╔═╡ Cell order:
# ╠═090f6f0f-fd67-419b-a35e-5392a59f39b1
# ╟─e89c666e-b0e7-11eb-22cf-7f6fd1e299bb
# ╠═ec540096-d9d1-4b5f-b6a2-16c470322424
# ╠═4e20e380-e928-4764-a982-3c0ca6bb942c
# ╠═7a68015d-0282-4cb0-bebc-c265d5ca72d3
# ╠═c4d6ad5b-f500-4435-a8eb-76c1f33f9775
# ╠═8f2717d2-ed5e-4e4e-abcd-6ff2982d8cd6
# ╠═58f76b47-542f-438b-ac71-694a76d9d77b
# ╠═9a03d3c5-3bc6-48c7-87cc-49ef6cbd31d4
# ╠═b10b13e9-de70-4c09-8b7a-9cd530d1e3ac
# ╠═2d9e27ff-5557-4c7e-b6a8-b3accd15f7f0
# ╠═26342478-769e-4ef1-a3f3-ddf5aad2e281
# ╠═a470ee13-e98c-423b-8398-e6c8774983c5
# ╠═dd4c2b0b-d8af-411f-9290-d36398772267
# ╠═a0681275-076f-400f-8862-74c0cba5ffc1
# ╠═2144ecf3-79c2-4de1-8aa4-52b4e6b57208
# ╠═d103a1b1-7ade-41da-bc1e-b01337333675
# ╠═5a13d4d9-7451-4fb0-8aaf-8d47947c10e4
# ╠═de4134e2-cc3f-4872-9244-4596d8525292
# ╠═6b671ffc-ccc1-4b84-9ddb-7f7137995348
# ╠═a33d8ddb-5ef5-4c0c-891c-eef1368fbcc9
# ╠═cbcd58ec-5f81-4ca4-8a4f-2a570db79f2a
# ╠═2d3ab22d-11fe-4ee3-9ee1-9410ae03498a
# ╟─bd64f7b1-4074-4eeb-925d-eb8cc7a524bf
# ╟─aa05c55c-e79d-4846-95e1-823f0a5ea2e3
# ╠═1ef77364-44a3-4819-b0c3-11cca842e53c
# ╠═fb47e75f-a3a0-4ee9-97c8-8a7e3d3f7831
# ╟─bd67cd4e-af24-403c-a3ba-9753f11dc21c

ext: jl


text:
@recipe function f(mach::MLJBase.Machine{<:EitherTunedModel})
    rep = report(mach)
    measurement = string(typeof(rep.best_history_entry.measure[1]))
    r = rep.plotting
    z = r.measurements
    X = r.parameter_values
    guides = r.parameter_names
    scales = r.parameter_scales
    n = size(X, 2)
    indices = LinearIndices((n, n))'

    layout := (n, n)
    size --> (275*n, 250*n)
    label := ""

    framestyle := :none
    for i in 1:n-1, j in i+1:n
        @series begin
            subplot := indices[i, j]
            ()
        end
    end

    framestyle := :box
    seriestype := :scatter
    for i in 1:n
        x = X[:, i]
        xsc = scales[i]
        xsc == :none && (x = string.(x))
        xscale := (xsc in (:custom, :linear, :none) ? :identity : xsc)

        @series begin
            subplot := indices[i, i]
            xguide := (i == n) ? guides[i] : ""
            yguide := measurement
            x, z
        end

        for j in 1:i-1
            y = X[:, j]
            ysc = scales[j]
            ysc == :none && (y = string.(y))
            ms = get(plotattributes, :markersize, 4)

            @series begin
                subplot := indices[i, j]
                yscale := (ysc in (:custom, :linear, :none) ? :identity : ysc)
                xguide := (i == n) ? guides[j] : ""
                yguide := (j == 1) ? guides[i] : ""
                markersize := _getmarkersize(ms, z)
                marker_z := z
                colorbar := false
                y, x
            end
        end
    end
end

function _getmarkersize(ms, z)
    ret = sqrt.(abs.(z))
    minz, maxz = extrema(x for x in ret if isfinite(x))
    ret .-= minz
    ret ./= maxz - minz
    4ms .* ret .+ 1
end

max_stars_repo_name: alan-turing-institute/MLJTuning

max_issues_repo_name: alan-turing-institute/MLJTuning

max_forks_repo_name: alan-turing-institute/MLJTuning


text:
using tele
using NCDatasets, Revise, Statistics, PyPlot
using CFTime, DSP, FFTW
using DataFrames, NaNMath
#analyzing data from https://psl.noaa.gov/data/gridded/data.noaa.ersst.v5.html

#FILE PATHS
data_dir = "/home/brynn/Documents/JP/860/tsa_data/"
ersst_p = data_dir * "sst.mnmean.nc"
# cesm_p2 = data_dir * "tos_Omon_CCSM4_past1000_r1i1p1_135001-185012.nc"
cesm_p2 = data_dir * "tos_Omon_CCSM4_historical_r1i1p1_185001-200512.nc"

#READ NETCDFS TO BOUNDING BOX
bbox = [71, 55, 290, 350]
(late, lone, timee, datae) = read_netcdf(ersst_p, "sst", bbox)
(latc, lonc, timec, datac) = read_netcdf(cesm_p2, "tos", bbox)

#FILTER CHARACTERISTICS
responsetype = Lowpass(0.1*10^(-7),fs = 3.80265176 * 10^(-7))
designmethod = Butterworth(4)

#Compute raw average time series over bounding box
datac = datac .- 273.15 #convert to celsius
tavge_raw = spatial_mean(datae)
tavgc_raw = spatial_mean(datac)

#ersst
dataep = zeros(length(timee), length(lone), length(late))
permutedims!(dataep, datae, [3,1,2])
filte = filt(digitalfilter(responsetype, designmethod), dataep)
permutedims!(datae, filte, [2,3,1])
anome = anomaly(datae, timee)
#filte holds filtered anomalized data

#cesm
timec = reinterpret(typeof(timee[begin]), timec) #convert from DatetimeNoLeap to DateTime

datacp = zeros(length(timec), length(lonc), length(latc))
permutedims!(datacp, datac, [3,1,2])
filtc = filt(digitalfilter(responsetype, designmethod), datacp)
permutedims!(datac, filtc, [2,3,1])
anomc = anomaly(datac, timec)
#filtc holds filtered, anomalized data

#filtered, anomalized time series
tavge_anom = spatial_mean(anome)
tavgc_anom = spatial_mean(anomc)

#Make contour plots
savetype = ".pdf"
contplot(latc,lonc,temp_mean(datac), "images/datac"*savetype)
contplot(late,lone,temp_mean(datae), "images/datae"*savetype)
contplot(latc,lonc,temp_mean(anomc), "images/anomc"*savetype)
contplot(late,lone,temp_mean(anome), "images/anome"*savetype)

#Make time series plots
filtseriesplot(tavgc_raw, tavgc_anom, timec, tavge_raw, tavge_anom, timee, "", "images/timeseries"*savetype)
# filtseriesplot(tavgc_anom, tavgc_anom_filt, timec, tavge_anom, tavge_anom_filt, timee, "", "images/anom"*savetype)

#COMPUTE FFTs
ffte_raw = fft(tavge_raw)
ffte_anom = fft(tavge_anom)
ffte_freq = FFTW.fftfreq(length(timee), 3.80265176 * 10^(-7))

fftc_raw = fft(tavgc_raw)
fftc_anom = fft(tavgc_anom)
fftc_freq = FFTW.fftfreq(length(timec), 3.80265176 * 10^(-7))

#PLot FFT
fftplot(ffte_freq[:], [ffte_raw, ffte_anom], "images/ffte"*savetype, "ERSST")
fftplot(fftc_freq[:], [fftc_raw, fftc_anom], "images/fftc"*savetype, "CESM")

#Plot FFT
figure()
ax = subplot()
halfway = 1006
plot(ffte_freq[begin+1:halfway], abs.(ffte_anom)[begin+1:halfway], alpha = 0.5)
halfway = 936
plot(fftc_freq[begin+1:halfway], abs.(fftc_anom)[begin+1:halfway], alpha = 0.5)
text(7 * 10 ^-10, 400, "f = 5.5 × 10^(-10) Hz")
text(7 * 10 ^-10, 380, "T = 47 years")
ax.set_xlim(ffte_freq[begin+1], 2 * 0.1*10^(-7))
ax.set_ylim(0, 500)
xlabel("Frequency [Hz]")
ylabel("Intensity")
legend(["ERSST", "CESM"])
savefig("images/fftec"*savetype)
close()

#histogram
#fig2b - boxplotGPCC Total # of Precips Obs
figure(figsize = (6, 8))
boxdict = boxplot([tavge_raw, tavgc_raw], notch = true, labels =  ["ERSST", "CESM (1850 - 2005)"])
grid(true)
ylabel("SST [°C]")
savefig("images/boxes_ersst"*savetype)
close()

#Extreme values analysis
df = DataFrame(min = [], max = [], mean = [], median = [], std = [])
push!(df, character(tavge_raw))
push!(df, character(tavgc_raw))

quant = quantile(tavgc_raw, [0.25, 0.5, 0.75])
thresh = (quant[3] - quant[1])* 1 + quant[3]
ind = findall(>(thresh), tavgc_raw)

extreme1 = temp_mean(datac[:,:, ind])
extreme2 = temp_mean(anomc[:,:, ind])

#Pot extreme values
contplot(latc,lonc,extreme1, "images/ext_raw"*savetype)
contplot(latc,lonc,extreme2, "images/ext_anom"*savetype)


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

@doc raw"""Represents a projected volume source

    IoK8sApiCoreV1ProjectedVolumeSource(;
        defaultMode=nothing,
        sources=nothing,
    )

    - defaultMode::Int32 : Mode bits to use on created files by default. Must be a value between 0 and 0777. Directories within the path are not affected by this setting. This might be in conflict with other options that affect the file mode, like fsGroup, and the result can be other mode bits set.
    - sources::Vector{IoK8sApiCoreV1VolumeProjection} : list of volume projections
"""
mutable struct IoK8sApiCoreV1ProjectedVolumeSource <: SwaggerModel
    defaultMode::Any # spec type: Union{ Nothing, Int32 } # spec name: defaultMode
    sources::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1VolumeProjection} } # spec name: sources

    function IoK8sApiCoreV1ProjectedVolumeSource(;defaultMode=nothing, sources=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1ProjectedVolumeSource, Symbol("defaultMode"), defaultMode)
        setfield!(o, Symbol("defaultMode"), defaultMode)
        validate_property(IoK8sApiCoreV1ProjectedVolumeSource, Symbol("sources"), sources)
        setfield!(o, Symbol("sources"), sources)
        o
    end
end # type IoK8sApiCoreV1ProjectedVolumeSource

const _property_map_IoK8sApiCoreV1ProjectedVolumeSource = Dict{Symbol,Symbol}(Symbol("defaultMode")=>Symbol("defaultMode"), Symbol("sources")=>Symbol("sources"))
const _property_types_IoK8sApiCoreV1ProjectedVolumeSource = Dict{Symbol,String}(Symbol("defaultMode")=>"Int32", Symbol("sources")=>"Vector{IoK8sApiCoreV1VolumeProjection}")
Base.propertynames(::Type{ IoK8sApiCoreV1ProjectedVolumeSource }) = collect(keys(_property_map_IoK8sApiCoreV1ProjectedVolumeSource))
Swagger.property_type(::Type{ IoK8sApiCoreV1ProjectedVolumeSource }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1ProjectedVolumeSource[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1ProjectedVolumeSource }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1ProjectedVolumeSource[property_name]

function check_required(o::IoK8sApiCoreV1ProjectedVolumeSource)
    (getproperty(o, Symbol("sources")) === nothing) && (return false)
    true
end

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


text:
function build_tree(p::DESPOTAlphaPlanner{S,A}, b0) where {S, A}
    D = DESPOTAlphaTree(p, b0)
    b = 1
    trial = 1
    start = CPUtime_us()

    Depth = sizehint!(Int[], 2000)
    sol = solver(p)
    while D.u[1]-D.l[1] > sol.epsilon_0 &&
          CPUtime_us()-start < sol.T_max*1e6 &&
          trial <= sol.max_trials
        b, as = explore!(D, 1, A[], p)
        backup!(D, b, as, p)
        push!(Depth, D.Delta[b])
        trial += 1
    end
    if (CPUtime_us()-start)*1e-6 > sol.T_max*sol.timeout_warning_threshold
        @warn(@sprintf("Surpass the time limit. The actual runtime is %3.1fs.
        Hyperparameters: K=%3d, C=%3d, bounds=%s",
        (CPUtime_us()-start)*1e-6, sol.K, sol.C, typeof(sol.bounds)))
        info_analysis(Dict(:tree=>D, :depth=>Depth))
    end
    return D, Depth
end

function explore!(D::DESPOTAlphaTree{S,A}, b::Int, as::Vector{A}, p::DESPOTAlphaPlanner{S,A}) where {S,A}
    sol = solver(p)
    while D.Delta[b] < sol.max_depth &&
        excess_uncertainty(D, b, sol.xi, p) > 0.0

        if isempty(D.children[b]) # a leaf
            expand!(D, b, as, p)
        end
        ba, b = next_best(D, b, p)
        push!(as, D.ba_action[ba])
    end
    if D.Delta[b] == sol.max_depth
        D.u[b] = D.l[b]
    end
    return b, as
end

function backup!(D::DESPOTAlphaTree, b::Int, as::Vector, p::DESPOTAlphaPlanner)
    disc = discount(p.pomdp)
    sol = solver(p)

    @assert b != 1
    ba = D.parent[b]
    as_ba = D.as_map[as]
    for sib in D.ba_children[ba]
        UpdateSiblings(D, b, as_ba, sib)
    end
    b = D.ba_parent[ba]
    while true
        r = D.weights[b]' * D.r[as_ba]
        D.ba_l[ba] = r + disc * sum(D.l[bp] * D.obs_prob[bp] for bp in D.ba_children[ba])
        D.ba_u[ba] = r + disc * sum(D.u[bp] * D.obs_prob[bp] for bp in D.ba_children[ba])

        if D.ba_l[ba] > D.l[b]
            D.l[b] = D.ba_l[ba]
            new_alpha = zeros(sol.K)
            for spi in 1:sol.K
                for (oi, bp) in enumerate(D.ba_children[ba])
                    new_alpha[spi] += D.L[as_ba][oi, spi] * D.alpha[bp][spi]
                end
                new_alpha[spi] *= disc
            end
            new_alpha .+= D.r[as_ba]
            D.alpha[b] = new_alpha
        end

        if D.ba_u[ba] < D.u[b]
            D.u[b] = maximum(D.ba_u[ba] for ba in D.children[b])
        end

        if sol.pp_update
            len_as = length(as)
            as_b = D.as_map[view(as, 1:(len_as-1))]
            new_pp = fill(-Inf, sol.K)
            for ba in D.children[b]
                as[len_as] = D.ba_action[ba]
                new_pp .= max.(new_pp, D.r[as_b] .+ disc .* D.pp_alpha_upper[D.as_map[as]])
            end
            D.pp_alpha_upper[as_b] .= min.(D.pp_alpha_upper[as_b], new_pp)
        end

        if b == 1
            break
        end
        as = view(as, 1:length(as)-1)
        ba = D.parent[b]
        as_ba = D.as_map[as]
        for sib in D.ba_children[ba]
            UpdateSiblings(D, b, as_ba, sib)
        end
        b = D.ba_parent[ba]
    end
    return nothing
end

function UpdateSiblings(D::DESPOTAlphaTree, trial::Int, as::Int, sib::Int)
    new_l = D.alpha[trial]' * D.weights[sib]
    if new_l > D.l[sib]
        D.l[sib] = new_l
        D.alpha[sib] = D.alpha[trial]
    end
    # @assert D.u[trial] >= D.weights[trial]' * D.pp_alpha_upper[as]
    theta = Inf
    for si in eachindex(D.weights[trial])
        new_theta = D.weights[sib][si] / (D.weights[trial][si] + 1e-10)
        if new_theta < theta
            theta = new_theta
        end
    end
    new_u = theta * D.u[trial]
    for si in eachindex(D.weights[trial])
        new_u += (D.weights[sib][si] - theta * D.weights[trial][si]) * D.pp_alpha_upper[as][si]
    end
    if new_u < D.u[sib]
        D.u[sib] = new_u
    end
end

function next_best(D::DESPOTAlphaTree, b::Int, p::DESPOTAlphaPlanner)
    max_u = -Inf
    best_ba = first(D.children[b])
    @inbounds for ba in D.children[b]
        if D.ba_u[ba] > max_u
            max_u = D.ba_u[ba]
            best_ba = ba
        end
    end

    max_eu = -Inf
    best_bp = first(D.ba_children[best_ba])
    tolerated_gap = solver(p).xi * max(D.u[1]-D.l[1], 0.0) / p.discounts[D.Delta[best_bp]+1]
    @inbounds for bp in D.ba_children[best_ba]
        eu = D.obs_prob[bp] * (D.u[bp]-D.l[bp] - tolerated_gap)
        if eu > max_eu
            max_eu = eu
            best_bp = bp
        end
    end

    return best_ba, best_bp
end

function excess_uncertainty(D::DESPOTAlphaTree, b::Int, ξ::Float64, p::DESPOTAlphaPlanner)
    return D.obs_prob[b] * (D.u[b]-D.l[b] - ξ * max(D.u[1]-D.l[1], 0.0) / p.discounts[D.Delta[b]+1])
end


text:
function precompile_constraint(model, F, S)
    Base.precompile(get, (model, ListOfConstraintIndices{F,S}))
    Base.precompile(get, (model, ListOfConstraintAttributesSet{F,S}))
    Base.precompile(get, (model, ConstraintSet, ConstraintIndex{F,S}))
    Base.precompile(get, (model, ConstraintSet, Vector{ConstraintIndex{F,S}}))
    Base.precompile(set, (model, ConstraintSet, ConstraintIndex{F,S}, S))

    Base.precompile(get, (model, ConstraintFunction, ConstraintIndex{F,S}))
    Base.precompile(set, (model, ConstraintFunction, ConstraintIndex{F,S}, F))
    Base.precompile(
        get,
        (model, ConstraintFunction, Vector{ConstraintIndex{F,S}}),
    )

    Base.precompile(get, (model, ConstraintDual, ConstraintIndex{F,S}))
    Base.precompile(get, (model, ConstraintPrimal, ConstraintIndex{F,S}))

    Base.precompile(
        get,
        (model, AbstractConstraintAttribute, ConstraintIndex{F,S}),
    )
    Base.precompile(
        get,
        (model, AbstractConstraintAttribute, Vector{ConstraintIndex{F,S}}),
    )

    Base.precompile(add_constraint, (model, F, S))
    Base.precompile(add_constraints, (model, Vector{F}, Vector{S}))
    Base.precompile(delete, (model, ConstraintIndex{F,S}))
    Base.precompile(is_valid, (model, ConstraintIndex{F,S}))
    Base.precompile(get, (model, ConstraintName, ConstraintIndex{F,S}))
    return Base.precompile(
        set,
        (model, ConstraintName, ConstraintIndex{F,S}, String),
    )
end

function precompile_variables(model)
    Base.precompile(delete, (model, VariableIndex))
    Base.precompile(delete, (model, Vector{VariableIndex}))
    Base.precompile(get, (model, AbstractVariableAttribute, VariableIndex))
    Base.precompile(
        get,
        (model, AbstractVariableAttribute, Vector{VariableIndex}),
    )
    Base.precompile(get, (model, VariableName, VariableIndex))
    Base.precompile(set, (model, VariableName, VariableIndex, String))
    Base.precompile(get, (model, VariablePrimalStart, VariableIndex))
    Base.precompile(set, (model, VariablePrimalStart, VariableIndex, Float64))
    Base.precompile(get, (model, VariablePrimalStart, Vector{VariableIndex}))
    Base.precompile(
        set,
        (model, VariablePrimalStart, Vector{VariableIndex}, Vector{Float64}),
    )
    Base.precompile(get, (model, VariablePrimal, VariableIndex))
    Base.precompile(get, (model, VariablePrimal, Vector{VariableIndex}))
    return Base.precompile(add_constrained_variables, (model, Reals))
end

function precompile_model(model, constraints)
    Base.precompile(empty!, (model,))
    Base.precompile(is_empty, (model,))
    Base.precompile(get, (model, ListOfConstraintTypesPresent))
    Base.precompile(optimize!, (model,))
    Base.precompile(add_variable, (model,))
    Base.precompile(add_variables, (model, Int))
    for attr in (
        ListOfVariableIndices,
        ListOfVariableAttributesSet,
        TerminationStatus,
        DualStatus,
        PrimalStatus,
        ObjectiveValue,
        Silent,
        TimeLimitSec,
        NumberOfVariables,
    )
        Base.precompile(get, (model, attr))
    end

    precompile_variables(model)
    for (F, S) in constraints
        precompile_constraint(model, F, S)
    end
    return Base.precompile(Tuple{typeof(add_constrained_variables),model,Reals})
end


text:
#=
# 2D heat equation, Dirichlet bc, CG
=#

### If the Femshop package has already been added, use this line #########
using Femshop # Note: to add the package, first do: ]add "https://github.com/paralab/femshop.git"

### If not, use these four lines (working from the examples directory) ###
# if !@isdefined(Femshop)
#     include("../Femshop.jl");
#     using .Femshop
# end
##########################################################################

init_femshop("heat2d");

# Optionally generate a log
useLog("heat2dlog")

# Set up the configuration (order doesn't matter)
domain(2)                   # dimension
solverType(CG)              # Use CG solver (default)
functionSpace(order=2)      # basis polynomial order
nodeType(LOBATTO)           # GLL elemental node arrangement (default)
timeStepper(RK4)            # time stepper (optional second arg is CFL#)

# Specify the problem
mesh(QUADMESH, elsperdim=10) # 10x10 elements

u = variable("u")            # make a scalar variable with symbol u
testSymbol("v")              # sets the symbol for a test function

T = 1;
timeInterval(T)              # The time interval is 0 to T
initial(u, "abs(x-0.5)+abs(y-0.5) < 0.2 ? 1 : 0") # initial condition

boundary(u, 1, DIRICHLET, 0)  # boundary condition for BID 1 is Dirichlet with value 0

# Write the weak form
coefficient("f", "0.5*sin(6*pi*x)*sin(6*pi*y)")
weakForm(u, "Dt(u*v) + 0.01 * dot(grad(u),grad(v)) - f*v")

solve(u);

# solution is stored in the variable's "values"
# using Plots
# pyplot();
# display(plot(Femshop.grid_data.allnodes[1,:], Femshop.grid_data.allnodes[2,:], u.values[:], st = :surface))

# Dump things to the log if desired
log_dump_config();
log_dump_prob();

finalize_femshop() # Finish writing and close any files


text:
module ClusterOrderTools

using Clustering, Random, Statistics
using ProgressMeter

export kmeansorder, evalclustering, average_heatmap, average_heatmap_vert

function kmeansorder(X, k, sr=1618)
    if k == 1
       return (KSI=[1], A=ones(size(X, 2)), SI=1:size(X, 2), k=1)
    end
    Random.seed!(sr)
    KM = kmeans(X, k)
    KSI = sortperm([findmax(@view KM.centers[:, i])[2] for i = 1:k], rev=false)
    KSSI = sortperm(KSI)
    A = KSSI[KM.assignments]
    SI = sortperm(A);

    (KSI=KSI, A=A, SI=SI, k=k, KM=KM, X=X)
end

function evalclustering(X, ks=2:2:40)
    KMK = @showprogress [kmeansorder(X, k) for k in ks]

    A = getfield.(KMK, :A)

    ri = [randindex(A[i], A[i+1]) for i = 1:(length(ks)-1)]
    vi = [varinfo(ks[i], A[i], ks[i+1], A[i+1]) for i = 1:(length(ks)-1)]
    ci = sum.(getfield.(getfield.(KMK, :KM), :costs))
    rtot = sum((X .- mean(X, dims=2)).^2)
    r² = 1 .- ci./rtot
    n, m = size(X)

    bic = n*log.(ci./(n*m)) + ks.*log(n)

    aic = ci + 2*size(X, 1).*ks

    (ks=ks, randindex=ri, varinfo=vi, costs=ci, rtot =rtot, rsquared=r², bic=bic, aic=aic)
end

function average_heatmap(H::Vector{T}, δ) where {T}
    n = length(H)

    w = cld(n, δ)
    AH = zeros(Float64, w)
    te = zeros(Int, w)
    for i = 1:n
        wi = cld(i, δ)
        te[wi] += 1
        AH[wi] += H[i]
    end
    if te[end] < 0.75δ
        AH[end] += AH[end-1]
        te[end] += te[end-1]
    end
    AH./te
end

average_heatmap(H, δh, δv) = average_heatmap_vert(average_heatmap(H, δh), δv)

function average_heatmap_vert(H, δ=2)

    n, m = size(H)

    w = cld(n, δ)
    AH = zeros(Float64, w, m)
    te = zeros(Int, w)
    for i = 1:m
        for j = 1:n
            wj = cld(j, δ)
            (i == 1) && (te[wj] += 1)
            AH[wj, i] += H[j, i]
        end
    end

    if te[end] < 0.75δ
        for i = 1:m
            AH[end, i] += AH[end-1, i]
        end
        te[end] += te[end-1]
    end

    AH./te
end

function average_heatmap(H, δ=2)

    n, m = size(H)
    w = cld(m, δ)

    AH = zeros(Float64, n, w)
    te = zeros(Int, w)
    for i = 1:m
        wi = cld(i, δ)
        te[wi] += 1
        for j = 1:n
            AH[j, wi] += H[j, i]
        end
    end
    if te[end] < 0.75δ
        for j = 1:n
            AH[j, end] += AH[j, end-1]
        end
        te[end] += te[end-1]
    end
    AH./te'
end

end # module

max_stars_repo_name: owensnick/ClusterOrderTools

max_issues_repo_name: owensnick/ClusterOrderTools

max_forks_repo_name: owensnick/ClusterOrderTools