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

# significant contributions by @jheinen

import GR
export GR

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

const gr_linetype = KW(
    :auto => 1,
    :solid => 1,
    :dash => 2,
    :dot => 3,
    :dashdot => 4,
    :dashdotdot => -1
)

const gr_markertype = KW(
    :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
)

const gr_halign = KW(
    :left => 1,
    :hcenter => 2,
    :right => 3
)

const gr_valign = KW(
    :top => 1,
    :vcenter => 3,
    :bottom => 5
)

const gr_font_family = Dict(
    "times" => 1,
    "helvetica" => 5,
    "courier" => 9,
    "bookman" => 14,
    "newcenturyschlbk" => 18,
    "avantgarde" => 22,
    "palatino" => 26
)

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

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_textcolor(c)   = GR.settextcolorind(gr_getcolorind(_cycle(c,1)))
gr_set_transparency(α::Real) = GR.settransparency(clamp(α, 0, 1))
function gr_set_transparency(::Nothing) 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)
    iend = 0
    n = length(x)
    while iend < n-1
        # set istart to the first index that is finite
        istart = -1
        for j = 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 = 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.drawarrow(x[iend-1], y[iend-1], x[iend], y[iend])
            end
            if arrowside in (:tail,:both)
                GR.drawarrow(x[istart+1], y[istart+1], x[istart], y[istart])
            end
        else
            break
        end
    end
end

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

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

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

function gr_text(x, y, s)
    if length(s) >= 2 && s[1] == '$' && s[end] == '$'
        GR.mathtex(x, y, s[2:end-1])
    elseif findfirst(isequal('\\'), s) != nothing || 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(yaxis)
    if yaxis[:formatter] in (:scientific, :auto) && yaxis[:ticks] in (:auto, :native)
        rtick_labels = convert_sci_unicode.(rtick_labels)
    end

    #draw angular grid
    if xaxis[:grid]
        gr_set_line(xaxis[:gridlinewidth], xaxis[:gridstyle], xaxis[:foreground_color_grid])
        gr_set_transparency(xaxis[:gridalpha])
        for i in 1:length(α)
            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])
        gr_set_transparency(yaxis[:gridalpha])
        for i in 1:length(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 1:length(α)
            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 1:length(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[:xaxis])
gr_y_axislims(sp::Subplot) = axis_limits(sp[:yaxis])
gr_z_axislims(sp::Subplot) = axis_limits(sp[:zaxis])
gr_xy_axislims(sp::Subplot) = gr_x_axislims(sp)..., gr_y_axislims(sp)...

function gr_lims(axis::Axis, adjust::Bool, expand = nothing)
    if expand != nothing
        expand_extrema!(axis, expand)
    end
    lims = axis_limits(axis)
    if adjust
        GR.adjustrange(lims...)
    else
        lims
    end
end

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

normalize_zvals(args...) = nothing
function normalize_zvals(zv::AVec, clims::NTuple{2, <:Real})
    vmin, vmax = ignorenan_extrema(zv)
    isfinite(clims[1]) && (vmin = clims[1])
    isfinite(clims[2]) && (vmax = clims[2])
    if vmin == vmax
        zeros(length(zv))
    else
        clamp.((zv .- vmin) ./ (vmax .- vmin), 0, 1)
    end
end

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

# draw ONE Shape
function gr_draw_marker(xi, yi, msize, shape::Shape)
    sx, sy = coords(shape)
    # convert to ndc coords (percentages of window)
    GR.selntran(0)
    w, h = gr_plot_size
    f = msize / (w + h)
    xi, yi = GR.wctondc(xi, yi)
    GR.fillarea(xi .+ sx .* f,
                yi .+ sy .* f)
    GR.selntran(1)
end

# draw ONE symbol marker
function gr_draw_marker(xi, yi, msize::Number, shape::Symbol)
    GR.setmarkertype(gr_markertype[shape])
    w, h = gr_plot_size
    GR.setmarkersize(0.3msize / ((w + h) * 0.001))
    GR.polymarker([xi], [yi])
end

# draw the markers, one at a time
function gr_draw_markers(series::Series, x, y, clims, msize = series[:markersize])

    isempty(x) && return
    GR.setfillintstyle(GR.INTSTYLE_SOLID)

    shapes = series[:markershape]
    if shapes != :none
        for i=1:length(x)
            msi = _cycle(msize, i)
            shape = _cycle(shapes, i)
            cfunc = isa(shape, Shape) ? gr_set_fillcolor : gr_set_markercolor

            # draw a filled in shape, slightly bigger, to estimate a stroke
            if series[:markerstrokewidth] > 0
                cfunc(get_markerstrokecolor(series, i))
                gr_set_transparency(get_markerstrokealpha(series, i))
                gr_draw_marker(x[i], y[i], msi + series[:markerstrokewidth], shape)
            end

            # draw the shape - don't draw filled area if marker shape is 1D
            if !(shape in (:hline, :vline, :+, :x))
                cfunc(get_markercolor(series, clims, i))
                gr_set_transparency(get_markeralpha(series, i))
                gr_draw_marker(x[i], y[i], msi, shape)
            end
        end
    end
end

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

function gr_set_line(lw, style, c) #, a)
    GR.setlinetype(gr_linetype[style])
    w, h = gr_plot_size
    GR.setlinewidth(_gr_thickness_scaling[1] * max(0, lw / ((w + h) * 0.001)))
    gr_set_linecolor(c) #, a)
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)
const _gr_point_mult = 0.0018 * ones(1)
const _gr_thickness_scaling = ones(1)

# set the font attributes... assumes _gr_point_mult has been populated already
function gr_set_font(f::Font; halign = f.halign, valign = f.valign,
                              color = f.color, rotation = f.rotation)
    family = lowercase(f.family)
    GR.setcharheight(_gr_point_mult[1] * f.pointsize)
    GR.setcharup(sind(-rotation), cosd(-rotation))
    if haskey(gr_font_family, family)
        GR.settextfontprec(100 + gr_font_family[family], GR.TEXT_PRECISION_STRING)
    end
    gr_set_textcolor(color)
    GR.settextalign(gr_halign[halign], gr_valign[valign])
end

function gr_nans_to_infs!(z)
    for (i,zi) in enumerate(z)
        if zi == NaN
            z[i] = Inf
        end
    end
end

# --------------------------------------------------------------------------------------
# viewport plot area

# this stays constant for a given subplot while displaying that subplot.
# values are [xmin, xmax, ymin, ymax].  they range [0,1].
const viewport_plotarea = zeros(4)

# the size of the current plot in pixels
const gr_plot_size = [600.0, 400.0]

const gr_colorbar_ratio = 0.1

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 is3d(sp)
        vp = viewport[:]
        extent = min(vp[2] - vp[1], vp[4] - vp[3])
        viewport[1] = 0.5 * (vp[1] + vp[2] - extent)
        viewport[2] = 0.5 * (vp[1] + vp[2] + extent)
        viewport[3] = 0.5 * (vp[3] + vp[4] - extent)
        viewport[4] = 0.5 * (vp[3] + vp[4] + extent)
    end
    if hascolorbar(sp)
        viewport[2] -= gr_colorbar_ratio
    end
    viewport
end

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

# reset the viewport to the plot area
function gr_set_viewport_plotarea()
    GR.setviewport(
        viewport_plotarea[1],
        viewport_plotarea[2],
        viewport_plotarea[3],
        viewport_plotarea[4]
    )
end

function gr_set_viewport_polar()
    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

# add the colorbar
function gr_colorbar(sp::Subplot, clims)
    xmin, xmax = gr_xy_axislims(sp)[1:2]
    gr_set_viewport_cmap(sp)
    l = zeros(Int32, 1, 256)
    l[1,:] = Int[round(Int, _i) for _i in range(1000, stop=1255, length=256)]
    GR.setscale(0)
    GR.setwindow(xmin, xmax, clims[1], clims[2])
    GR.cellarray(xmin, xmax, clims[2], clims[1], 1, length(l), l)
    ztick = 0.5 * GR.tick(clims[1], clims[2])
    GR.axes(0, ztick, xmax, clims[1], 0, 1, 0.005)

    gr_set_font(guidefont(sp[:yaxis]))
    GR.settextalign(GR.TEXT_HALIGN_CENTER, GR.TEXT_VALIGN_TOP)
    GR.setcharup(-1, 0)
    gr_text(viewport_plotarea[2] + gr_colorbar_ratio,
            gr_view_ycenter(), sp[:colorbar_title])

    gr_set_viewport_plotarea()
end

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

function gr_legend_pos(s::Symbol,w,h)
    str = string(s)
    if str == "best"
        str = "topright"
    end
    if occursin("right", str)
        xpos = viewport_plotarea[2] - 0.05 - w
    elseif occursin("left", str)
        xpos = viewport_plotarea[1] + 0.11
    else
        xpos = (viewport_plotarea[2]-viewport_plotarea[1])/2 - w/2 +.04
    end
    if occursin("top", str)
        ypos = viewport_plotarea[4] - 0.06
    elseif occursin("bottom", str)
        ypos = viewport_plotarea[3] + h + 0.06
    else
        ypos = (viewport_plotarea[4]-viewport_plotarea[3])/2 + h/2
    end
    (xpos,ypos)
end

function gr_legend_pos(v::Tuple{S,T},w,h) 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

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

const _gr_gradient_alpha = ones(256)

function gr_set_gradient(c)
    grad = isa(c, ColorGradient) ? c : cgrad()
    for (i,z) in enumerate(range(0, stop=1, length=256))
        c = grad[z]
        GR.setcolorrep(999+i, red(c), green(c), blue(c))
        _gr_gradient_alpha[i] = alpha(c)
    end
    grad
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()

    _gr_thickness_scaling[1] = plt[:thickness_scaling]
    dpi_factor = plt[:dpi] / Plots.DPI
    if fmt == "svg"
        dpi_factor *= 4
    end

    # 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 = plt[:size]
    gr_plot_size[:] = [w, h]
    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])

    # update point mult
    px_per_pt = px / pt
    _gr_point_mult[1] = 1.5 * _gr_thickness_scaling[1] * px_per_pt / max(h,w)

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

    GR.updatews()
end

function gr_set_xticks_font(sp)
    flip = sp[:yaxis][:flip]
    mirror = sp[:xaxis][:mirror]
    gr_set_font(tickfont(sp[:xaxis]),
                halign = (:left, :hcenter, :right)[sign(sp[:xaxis][:rotation]) + 2],
                valign = (mirror ? :bottom : :top),
                rotation = sp[:xaxis][:rotation])
    return flip, mirror
end

function gr_set_yticks_font(sp)
    flip = sp[:xaxis][:flip]
    mirror = sp[:yaxis][:mirror]
    gr_set_font(tickfont(sp[:yaxis]),
                halign = (mirror ? :left : :right),
                valign = (:top, :vcenter, :bottom)[sign(sp[:yaxis][:rotation]) + 2],
                rotation = sp[:yaxis][:rotation])
    return flip, mirror
end

function gr_get_ticks_size(ticks, i)
    GR.savestate()
    GR.selntran(0)
    l = 0.0
    for (cv, dv) in zip(ticks...)
        tb = gr_inqtext(0, 0, string(dv))[i]
        tb_min, tb_max = extrema(tb)
        l = max(l, tb_max - tb_min)
    end
    GR.restorestate()
    return l
end

function _update_min_padding!(sp::Subplot{GRBackend})
    dpi = sp.plt[:thickness_scaling]
    if !haskey(ENV, "GKSwstype")
        if isijulia()
            ENV["GKSwstype"] = "svg"
        end
    end
    # Add margin given by the user
    leftpad   = 4mm  + sp[:left_margin]
    toppad    = 2mm  + sp[:top_margin]
    rightpad  = 4mm  + sp[:right_margin]
    bottompad = 2mm  + sp[:bottom_margin]
    # Add margin for title
    if sp[:title] != ""
        toppad += 5mm
    end
    # Add margin for x and y ticks
    xticks, yticks = axis_drawing_info(sp)[1:2]
    if !(xticks in (nothing, false, :none))
        flip, mirror = gr_set_xticks_font(sp)
        l = gr_get_ticks_size(xticks, 2)
        if mirror
            toppad += 1mm + gr_plot_size[2] * l * px
        else
            bottompad += 1mm + gr_plot_size[2] * l * px
        end
    end
    if !(yticks in (nothing, false, :none))
        flip, mirror = gr_set_yticks_font(sp)
        l = gr_get_ticks_size(yticks, 1)
        if mirror
            rightpad += 1mm + gr_plot_size[1] * l * px
        else
            leftpad += 1mm + gr_plot_size[1] * l * px
        end
    end
    # Add margin for x label
    if sp[:xaxis][:guide] != ""
        if sp[:xaxis][:guide_position] == :top || (sp[:xaxis][:guide_position] == :auto && sp[:xaxis][:mirror] == true)
            toppad += 4mm
        else
            bottompad += 4mm
        end
    end
    # Add margin for y label
    if sp[:yaxis][:guide] != ""
        if sp[:yaxis][:guide_position] == :right || (sp[:yaxis][:guide_position] == :auto && sp[:yaxis][:mirror] == true)
            rightpad += 4mm
        else
            leftpad += 4mm
        end
    end
    if sp[:colorbar_title] != ""
        rightpad += 4mm
    end
    sp.minpad = Tuple(dpi * [leftpad, toppad, rightpad, bottompad])
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)
    # get data limits
    data_lims = gr_xy_axislims(sp)
    xy_lims = data_lims

    ratio = sp[:aspect_ratio]
    if ratio != :none
        if ratio == :equal
            ratio = 1
        end
        viewport_ratio = (viewport_plotarea[2] - viewport_plotarea[1]) / (viewport_plotarea[4] - viewport_plotarea[3])
        window_ratio = (data_lims[2] - data_lims[1]) / (data_lims[4] - data_lims[3]) / 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

    # fill in the plot area background
    bg = plot_color(sp[:background_color_inside])
    gr_fill_viewport(viewport_plotarea, bg)

    # reduced from before... set some flags based on the series in this subplot
    # TODO: can these be generic flags?
    outside_ticks = false
    cmap = hascolorbar(sp)
    draw_axes = sp[:framestyle] != :none
    # axes_2d = true
    for series in series_list(sp)
        st = series[:seriestype]
        if st == :pie
            draw_axes = false
        end
        if st == :heatmap
            outside_ticks = true
            for ax in (sp[:xaxis], sp[:yaxis])
                v = series[ax[:letter]]
                if diff(collect(extrema(diff(v))))[1] > 1e-6*std(v)
                    @warn("GR: heatmap only supported with equally spaced data.")
                end
            end
            x, y = heatmap_edges(series[:x], sp[:xaxis][:scale]), heatmap_edges(series[:y], sp[:yaxis][:scale])
            xy_lims = x[1], x[end], y[1], y[end]
            expand_extrema!(sp[:xaxis], x)
            expand_extrema!(sp[:yaxis], y)
            data_lims = gr_xy_axislims(sp)
        end
    end

    # set our plot area view
    gr_set_viewport_plotarea()

    # these are the Axis objects, which hold scale, lims, etc
    xaxis = sp[:xaxis]
    yaxis = sp[:yaxis]
    zaxis = sp[:zaxis]

    # set the scale flags and window
    xmin, xmax, ymin, ymax = data_lims
    scaleop = 0
    xtick, ytick = 1, 1
    if xmax > xmin && ymax > ymin
        # NOTE: for log axes, the major_x and major_y - if non-zero (omit labels) - control the minor grid lines (1 = draw 9 minor grid lines, 2 = no minor grid lines)
        # NOTE: for log axes, the x_tick and y_tick - if non-zero (omit axes) - only affect the output appearance (1 = nomal, 2 = scientiic notation)
        xaxis[:scale] == :log10 && (scaleop |= GR.OPTION_X_LOG)
        yaxis[:scale] == :log10 && (scaleop |= GR.OPTION_Y_LOG)
        xaxis[:flip]            && (scaleop |= GR.OPTION_FLIP_X)
        yaxis[:flip]            && (scaleop |= GR.OPTION_FLIP_Y)
        if scaleop & GR.OPTION_X_LOG == 0
            majorx = 1 #5
            xtick = GR.tick(xmin, xmax) / majorx
        else
            # log axis
            xtick = 2  # scientific notation
            majorx = 2 # no minor grid lines
        end
        if scaleop & GR.OPTION_Y_LOG == 0
            majory = 1 #5
            ytick = GR.tick(ymin, ymax) / majory
        else
            # log axis
            ytick = 2  # scientific notation
            majory = 2 # no minor grid lines
        end

        # NOTE: setwindow sets the "data coordinate" limits of the current "viewport"
        GR.setwindow(xmin, xmax, ymin, ymax)
        GR.setscale(scaleop)
    end

    # draw the axes
    gr_set_font(tickfont(xaxis))
    GR.setlinewidth(sp.plt[:thickness_scaling])

    if is3d(sp)
        zmin, zmax = gr_lims(zaxis, true)
        clims = sp[:clims]
        if is_2tuple(clims)
            isfinite(clims[1]) && (zmin = clims[1])
            isfinite(clims[2]) && (zmax = clims[2])
        end
        GR.setspace(zmin, zmax, round.(Int, sp[:camera])...)
        xtick = GR.tick(xmin, xmax) / 2
        ytick = GR.tick(ymin, ymax) / 2
        ztick = GR.tick(zmin, zmax) / 2
        ticksize = 0.01 * (viewport_plotarea[2] - viewport_plotarea[1])

        if xaxis[:grid]
            gr_set_line(xaxis[:gridlinewidth], xaxis[:gridstyle], xaxis[:foreground_color_grid])
            gr_set_transparency(xaxis[:gridalpha])
            GR.grid3d(xtick, 0, 0, xmin, ymax, zmin, 2, 0, 0)
        end
        if yaxis[:grid]
            gr_set_line(yaxis[:gridlinewidth], yaxis[:gridstyle], yaxis[:foreground_color_grid])
            gr_set_transparency(yaxis[:gridalpha])
            GR.grid3d(0, ytick, 0, xmin, ymax, zmin, 0, 2, 0)
        end
        if zaxis[:grid]
            gr_set_line(zaxis[:gridlinewidth], zaxis[:gridstyle], zaxis[:foreground_color_grid])
            gr_set_transparency(zaxis[:gridalpha])
            GR.grid3d(0, 0, ztick, xmin, ymax, zmin, 0, 0, 2)
        end
        gr_set_line(1, :solid, xaxis[:foreground_color_axis])
        gr_set_transparency(1)
        GR.axes3d(xtick, 0, ztick, xmin, ymin, zmin, 2, 0, 2, -ticksize)
        GR.axes3d(0, ytick, 0, xmax, ymin, zmin, 0, 2, 0, ticksize)

    elseif ispolar(sp)
        r = gr_set_viewport_polar()
        #rmin, rmax = GR.adjustrange(ignorenan_minimum(r), ignorenan_maximum(r))
        rmin, rmax = axis_limits(sp[:yaxis])
        gr_polaraxes(rmin, rmax, sp)

    elseif draw_axes
        if xmax > xmin && ymax > ymin
            GR.setwindow(xmin, xmax, ymin, ymax)
        end

        xticks, yticks, xspine_segs, yspine_segs, xtick_segs, ytick_segs, xgrid_segs, ygrid_segs, xminorgrid_segs, yminorgrid_segs, xborder_segs, yborder_segs = axis_drawing_info(sp)
        # @show xticks yticks #spine_segs grid_segs

        # draw the grid lines
        if xaxis[:grid]
            # gr_set_linecolor(sp[:foreground_color_grid])
            # GR.grid(xtick, ytick, 0, 0, majorx, majory)
            gr_set_line(xaxis[:gridlinewidth], xaxis[:gridstyle], xaxis[:foreground_color_grid])
            gr_set_transparency(xaxis[:gridalpha])
            gr_polyline(coords(xgrid_segs)...)
        end
        if yaxis[:grid]
            gr_set_line(yaxis[:gridlinewidth], yaxis[:gridstyle], yaxis[:foreground_color_grid])
            gr_set_transparency(yaxis[:gridalpha])
            gr_polyline(coords(ygrid_segs)...)
        end
        if xaxis[:minorgrid]
            # gr_set_linecolor(sp[:foreground_color_grid])
            # GR.grid(xtick, ytick, 0, 0, majorx, majory)
            gr_set_line(xaxis[:minorgridlinewidth], xaxis[:minorgridstyle], xaxis[:foreground_color_minor_grid])
            gr_set_transparency(xaxis[:minorgridalpha])
            gr_polyline(coords(xminorgrid_segs)...)
        end
        if yaxis[:minorgrid]
            gr_set_line(yaxis[:minorgridlinewidth], yaxis[:minorgridstyle], yaxis[:foreground_color_minor_grid])
            gr_set_transparency(yaxis[:minorgridalpha])
            gr_polyline(coords(yminorgrid_segs)...)
        end
        gr_set_transparency(1.0)

        # axis lines
        if xaxis[:showaxis]
            gr_set_line(1, :solid, xaxis[:foreground_color_border])
            GR.setclip(0)
            gr_polyline(coords(xspine_segs)...)
        end
        if yaxis[:showaxis]
            gr_set_line(1, :solid, yaxis[:foreground_color_border])
            GR.setclip(0)
            gr_polyline(coords(yspine_segs)...)
        end
        GR.setclip(1)

        # axis ticks
        if xaxis[:showaxis]
            if sp[:framestyle] in (:zerolines, :grid)
                gr_set_line(1, :solid, xaxis[:foreground_color_grid])
                gr_set_transparency(xaxis[:gridalpha])
            else
                gr_set_line(1, :solid, xaxis[:foreground_color_axis])
            end
            GR.setclip(0)
            gr_polyline(coords(xtick_segs)...)
        end
        if  yaxis[:showaxis]
            if sp[:framestyle] in (:zerolines, :grid)
                gr_set_line(1, :solid, yaxis[:foreground_color_grid])
                gr_set_transparency(yaxis[:gridalpha])
            else
                gr_set_line(1, :solid, yaxis[:foreground_color_axis])
            end
            GR.setclip(0)
            gr_polyline(coords(ytick_segs)...)
        end
        GR.setclip(1)

        # tick marks
        if !(xticks in (:none, nothing, false)) && xaxis[:showaxis]
            # x labels
            flip, mirror = gr_set_xticks_font(sp)
            for (cv, dv) in zip(xticks...)
                # use xor ($) to get the right y coords
                xi, yi = GR.wctondc(cv, sp[:framestyle] == :origin ? 0 : xor(flip, mirror) ? ymax : ymin)
                # @show cv dv ymin xi yi flip mirror (flip $ mirror)
                if xaxis[:ticks] in (:auto, :native)
                    # ensure correct dispatch in gr_text for automatic log ticks
                    if xaxis[:scale] in _logScales
                        dv = string(dv, "\\ ")
                    elseif xaxis[:formatter] in (:scientific, :auto)
                        dv = convert_sci_unicode(dv)
                    end
                end
                gr_text(xi, yi + (mirror ? 1 : -1) * 5e-3 * (xaxis[:tick_direction] == :out ? 1.5 : 1.0), string(dv))
            end
        end

        if !(yticks in (:none, nothing, false)) && yaxis[:showaxis]
            # y labels
            flip, mirror = gr_set_yticks_font(sp)
            for (cv, dv) in zip(yticks...)
                # use xor ($) to get the right y coords
                xi, yi = GR.wctondc(sp[:framestyle] == :origin ? 0 : xor(flip, mirror) ? xmax : xmin, cv)
                # @show cv dv xmin xi yi
                if yaxis[:ticks] in (:auto, :native)
                    # ensure correct dispatch in gr_text for automatic log ticks
                    if yaxis[:scale] in _logScales
                        dv = string(dv, "\\ ")
                    elseif yaxis[:formatter] in (:scientific, :auto)
                        dv = convert_sci_unicode(dv)
                    end
                end
                gr_text(xi + (mirror ? 1 : -1) * 1e-2 * (yaxis[:tick_direction] == :out ? 1.5 : 1.0), yi, string(dv))
            end
        end

        # border
        intensity = sp[:framestyle] == :semi ? 0.5 : 1.0
        if sp[:framestyle] in (:box, :semi)
            gr_set_line(intensity, :solid, xaxis[:foreground_color_border])
            gr_set_transparency(intensity)
            gr_polyline(coords(xborder_segs)...)
            gr_set_line(intensity, :solid, yaxis[:foreground_color_border])
            gr_set_transparency(intensity)
            gr_polyline(coords(yborder_segs)...)
        end
    end
    # end

    # add the guides
    GR.savestate()
    if sp[:title] != ""
        gr_set_font(titlefont(sp))
        loc = sp[:title_location]
        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()
            halign = GR.TEXT_HALIGN_CENTER
        end
        GR.settextalign(halign, GR.TEXT_VALIGN_TOP)
        gr_text(xpos, viewport_subplot[4], sp[:title])
    end

    if xaxis[:guide] != ""
        gr_set_font(guidefont(xaxis))
        if xaxis[:guide_position] == :top || (xaxis[:guide_position] == :auto && xaxis[:mirror] == true)
            GR.settextalign(GR.TEXT_HALIGN_CENTER, GR.TEXT_VALIGN_TOP)
            gr_text(gr_view_xcenter(), viewport_subplot[4], xaxis[:guide])
        else
            GR.settextalign(GR.TEXT_HALIGN_CENTER, GR.TEXT_VALIGN_BOTTOM)
            gr_text(gr_view_xcenter(), viewport_subplot[3], xaxis[:guide])
        end
    end

    if yaxis[:guide] != ""
        gr_set_font(guidefont(yaxis))
        GR.setcharup(-1, 0)
        if yaxis[:guide_position] == :right || (yaxis[:guide_position] == :auto && yaxis[:mirror] == true)
            GR.settextalign(GR.TEXT_HALIGN_CENTER, GR.TEXT_VALIGN_BOTTOM)
            gr_text(viewport_subplot[2], gr_view_ycenter(), yaxis[:guide])
        else
            GR.settextalign(GR.TEXT_HALIGN_CENTER, GR.TEXT_VALIGN_TOP)
            gr_text(viewport_subplot[1], gr_view_ycenter(), yaxis[:guide])
        end
    end
    GR.restorestate()

    gr_set_font(tickfont(xaxis))

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

    # calculate the colorbar limits once for a subplot
    clims = get_clims(sp)

    for (idx, series) in enumerate(series_list(sp))
        st = series[:seriestype]

        # update the current stored gradient
        if st in (:contour, :surface, :wireframe, :heatmap)
            gr_set_gradient(series[:fillcolor]) #, series[:fillalpha])
        elseif series[:marker_z] != nothing
            series[:markercolor] = gr_set_gradient(series[:markercolor])
        elseif series[:line_z] !=  nothing
            series[:linecolor] = gr_set_gradient(series[:linecolor])
        elseif series[:fill_z] != nothing
            series[:fillcolor] = gr_set_gradient(series[:fillcolor])
        end

        GR.savestate()

        # update the bounding window
        if ispolar(sp)
            gr_set_viewport_polar()
        else
            xmin, xmax, ymin, ymax = data_lims
            if xmax > xmin && ymax > ymin
                GR.setwindow(xmin, xmax, ymin, ymax)
            end
        end

        x, y, z = series[:x], series[:y], series[:z]
        frng = series[:fillrange]

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

        # recompute data
        if typeof(z) <: Surface
            z = vec(transpose_z(series, z.surf, false))
        elseif ispolar(sp)
            if frng != nothing
                _, frng = convert_to_polar(x, frng, (rmin, rmax))
            end
            x, y = convert_to_polar(x, y, (rmin, rmax))
        end

        if st == :straightline
            x, y = straightline_data(series)
        end

        if st in (:path, :scatter, :straightline)
            if length(x) > 1
                lz = series[:line_z]
                segments = iter_segments(series)
                # do area fill
                if frng != nothing
                    GR.setfillintstyle(GR.INTSTYLE_SOLID)
                    fr_from, fr_to = (is_2tuple(frng) ? frng : (y, frng))
                    for (i, rng) in enumerate(segments)
                        gr_set_fillcolor(get_fillcolor(series, clims, i))
                        fx = _cycle(x, vcat(rng, reverse(rng)))
                        fy = vcat(_cycle(fr_from,rng), _cycle(fr_to,reverse(rng)))
                        gr_set_transparency(get_fillalpha(series, i))
                        GR.fillarea(fx, fy)
                    end
                end

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

            if series[:markershape] != :none
                gr_draw_markers(series, x, y, clims)
            end

        elseif st == :contour
            zmin, zmax = clims
            GR.setspace(zmin, zmax, 0, 90)
            if typeof(series[:levels]) <: AbstractArray
                h = series[:levels]
            else
                h = series[:levels] > 1 ? range(zmin, stop=zmax, length=series[:levels]) : [(zmin + zmax) / 2]
            end
            if series[:fillrange] != nothing
                GR.surface(x, y, z, GR.OPTION_CELL_ARRAY)
            else
                GR.setlinetype(gr_linetype[get_linestyle(series)])
                GR.setlinewidth(max(0, get_linewidth(series) / (sum(gr_plot_size) * 0.001)))
                if plot_color(series[:linecolor]) == [plot_color(:black)]
                    GR.contour(x, y, h, z, 0 + (series[:contour_labels] == true ? 1 : 0))
                else
                    GR.contour(x, y, h, z, 1000 + (series[:contour_labels] == true ? 1 : 0))
                end
            end

            # create the colorbar of contour levels
            if cmap
                gr_set_line(1, :solid, yaxis[:foreground_color_axis])
                gr_set_viewport_cmap(sp)
                l = (length(h) > 1) ? round.(Int32, 1000 .+ (h .- ignorenan_minimum(h)) ./ (ignorenan_maximum(h) - ignorenan_minimum(h)) .* 255) : Int32[1000, 1255]
                GR.setwindow(xmin, xmax, zmin, zmax)
                GR.cellarray(xmin, xmax, zmax, zmin, 1, length(l), l)
                ztick = 0.5 * GR.tick(zmin, zmax)
                GR.axes(0, ztick, xmax, zmin, 0, 1, 0.005)
                gr_set_viewport_plotarea()
            end

        elseif st in [:surface, :wireframe]
            if st == :surface
                if length(x) == length(y) == length(z)
                    GR.trisurface(x, y, z)
                else
                    GR.gr3.surface(x, y, z, GR.OPTION_COLORED_MESH)
                end
            else
                GR.setfillcolorind(0)
                GR.surface(x, y, z, GR.OPTION_FILLED_MESH)
            end

        elseif st == :heatmap
            xmin, xmax, ymin, ymax = xy_lims
            zmin, zmax = clims
            m, n = length(x), length(y)
            xinds = sort(1:m, rev = xaxis[:flip])
            yinds = sort(1:n, rev = yaxis[:flip])
            z = reshape(reshape(z, m, n)[xinds, yinds], m*n)
            GR.setspace(zmin, zmax, 0, 90)
            grad = isa(series[:fillcolor], ColorGradient) ? series[:fillcolor] : cgrad()
            colors = [plot_color(grad[clamp((zi-zmin) / (zmax-zmin), 0, 1)], series[:fillalpha]) for zi=z]
            rgba = map(c -> UInt32( round(UInt, alpha(c) * 255) << 24 +
                                    round(UInt,  blue(c) * 255) << 16 +
                                    round(UInt, green(c) * 255) << 8  +
                                    round(UInt,   red(c) * 255) ), colors)
            w, h = length(x), length(y)
            GR.drawimage(xmin, xmax, ymax, ymin, w, h, rgba)

        elseif st in (:path3d, :scatter3d)
            # draw path
            if st == :path3d
                if length(x) > 1
                    lz = series[:line_z]
                    segments = iter_segments(series)
                    for (i, rng) in enumerate(segments)
                        gr_set_line(get_linewidth(series, i), get_linestyle(series, i), get_linecolor(series, clims, i)) #, series[:linealpha])
                        gr_set_transparency(get_linealpha(series, i))
                        GR.polyline3d(x[rng], y[rng], z[rng])
                    end
                end
            end

            # draw markers
            if st == :scatter3d || series[:markershape] != :none
                x2, y2 = unzip(map(GR.wc3towc, x, y, z))
                gr_draw_markers(series, x2, y2, clims)
            end

        # TODO: replace with pie recipe
        elseif st == :pie
            GR.selntran(0)
            GR.setfillintstyle(GR.INTSTYLE_SOLID)
            xmin, xmax, ymin, ymax = viewport_plotarea
            ymax -= 0.1 * (xmax - xmin)
            xcenter = 0.5 * (xmin + xmax)
            ycenter = 0.5 * (ymin + ymax)
            if xmax - xmin > ymax - ymin
                r = 0.5 * (ymax - ymin)
                xmin, xmax = xcenter - r, xcenter + r
            else
                r = 0.5 * (xmax - xmin)
                ymin, ymax = ycenter - r, ycenter + r
            end
            labels = pie_labels(sp, series)
            slices = series[:y]
            numslices = length(slices)
            total = sum(slices)
            a1 = 0
            x = zeros(3)
            y = zeros(3)
            for i in 1:numslices
                a2 = round(Int, a1 + (slices[i] / total) * 360.0)
                GR.setfillcolorind(980 + (i-1) % 20)
                GR.fillarc(xmin, xmax, ymin, ymax, a1, a2)
                α = 0.5 * (a1 + a2)
                cosf = r * cos(α * pi / 180)
                sinf = r * sin(α * pi / 180)
                x[1] = xcenter + cosf
                y[1] = ycenter + sinf
                x[2] = x[1] + 0.1 * cosf
                y[2] = y[1] + 0.1 * sinf
                y[3] = y[2]
                if 90 <= α < 270
                    x[3] = x[2] - 0.05
                    GR.settextalign(GR.TEXT_HALIGN_RIGHT, GR.TEXT_VALIGN_HALF)
                    gr_text(x[3] - 0.01, y[3], string(labels[i]))
                else
                    x[3] = x[2] + 0.05
                    GR.settextalign(GR.TEXT_HALIGN_LEFT, GR.TEXT_VALIGN_HALF)
                    gr_text(x[3] + 0.01, y[3], string(labels[i]))
                end
                gr_polyline(x, y)
                a1 = a2
            end
            GR.selntran(1)

        elseif st == :shape
            x, y = shape_data(series)
            for (i,rng) in enumerate(iter_segments(x, y))
                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
                    gr_set_fill(get_fillcolor(series, clims, i))
                    gr_set_transparency(get_fillalpha(series, i))
                    GR.fillarea(xseg, yseg)

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

        elseif st == :image
            z = transpose_z(series, series[:z].surf, true)'
            w, h = length(x), length(y)
            xinds = sort(1:w, rev = xaxis[:flip])
            yinds = sort(1:h, rev = yaxis[:flip])
            z = z[xinds, yinds]
            xmin, xmax = ignorenan_extrema(series[:x]); ymin, ymax = ignorenan_extrema(series[:y])
            if eltype(z) <: Colors.AbstractGray
                grey = round.(UInt8, float(z) * 255)
                rgba = map(c -> UInt32( 0xff000000 + UInt(c)<<16 + UInt(c)<<8 + UInt(c) ), grey)
            else
                rgba = map(c -> UInt32( round(UInt, alpha(c) * 255) << 24 +
                                        round(UInt,  blue(c) * 255) << 16 +
                                        round(UInt, green(c) * 255) << 8  +
                                        round(UInt,   red(c) * 255) ), z)
            end
            GR.drawimage(xmin, xmax, ymax, ymin, w, h, rgba)
        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)
            gr_text(GR.wctondc(xi, yi)..., str)
        end

        GR.restorestate()
    end

    # draw the colorbar
    GR.savestate()
    # special colorbar with steps is drawn for contours
    if cmap && any(series[:seriestype] != :contour for series in series_list(sp))
        gr_set_line(1, :solid, yaxis[:foreground_color_axis])
        gr_set_transparency(1)
        gr_colorbar(sp, clims)
    end
    GR.restorestate()

    # add the legend
    if sp[:legend] != :none
        GR.savestate()
        GR.selntran(0)
        GR.setscale(0)
        gr_set_font(legendfont(sp))
        w = 0
        i = 0
        n = 0
        if sp[:legendtitle] != nothing
            tbx, tby = gr_inqtext(0, 0, string(sp[:legendtitle]))
            w = tbx[3] - tbx[1]
            n += 1
        end
        for series in series_list(sp)
            should_add_to_legend(series) || continue
            n += 1
            if typeof(series[:label]) <: Array
                i += 1
                lab = series[:label][i]
            else
                lab = series[:label]
            end
            tbx, tby = gr_inqtext(0, 0, string(lab))
            w = max(w, tbx[3] - tbx[1])
        end
        if w > 0
            dy = _gr_point_mult[1] * sp[:legendfontsize] * 1.75
            h = dy*n
            (xpos,ypos) = gr_legend_pos(sp[:legend],w,h)
            GR.setfillintstyle(GR.INTSTYLE_SOLID)
            gr_set_fillcolor(sp[:background_color_legend])
            GR.fillrect(xpos - 0.08, xpos + w + 0.02, ypos + dy, ypos - dy * n)
            gr_set_line(1, :solid, sp[:foreground_color_legend])
            GR.drawrect(xpos - 0.08, xpos + w + 0.02, ypos + dy, ypos - dy * n)
            i = 0
            if sp[:legendtitle] != nothing
                GR.settextalign(GR.TEXT_HALIGN_CENTER, GR.TEXT_VALIGN_HALF)
                gr_set_textcolor(sp[:legendfontcolor])
                gr_set_transparency(1)
                gr_text(xpos - 0.03 + 0.5*w, ypos, string(sp[:legendtitle]))
                ypos -= dy
            end
            for series in series_list(sp)
                should_add_to_legend(series) || continue
                st = series[:seriestype]
                gr_set_line(get_linewidth(series), get_linestyle(series), get_linecolor(series, clims)) #, series[:linealpha])

                if (st == :shape || series[:fillrange] != nothing) && series[:ribbon] == nothing
                    gr_set_fill(get_fillcolor(series, clims)) #, series[:fillalpha])
                    l, r = xpos-0.07, xpos-0.01
                    b, t = ypos-0.4dy, ypos+0.4dy
                    x = [l, r, r, l, l]
                    y = [b, b, t, t, b]
                    gr_set_transparency(get_fillalpha(series))
                    gr_polyline(x, y, GR.fillarea)
                    gr_set_transparency(get_linealpha(series))
                    gr_set_line(get_linewidth(series), get_linestyle(series), get_linecolor(series, clims))
                    st == :shape && gr_polyline(x, y)
                end

                if st in (:path, :straightline)
                    gr_set_transparency(get_linealpha(series))
                    if series[:fillrange] == nothing || series[:ribbon] != nothing
                        GR.polyline([xpos - 0.07, xpos - 0.01], [ypos, ypos])
                    else
                        GR.polyline([xpos - 0.07, xpos - 0.01], [ypos+0.4dy, ypos+0.4dy])
                    end
                end

                if series[:markershape] != :none
                    gr_draw_markers(series, xpos - .035, ypos, clims, 6)
                end

                if typeof(series[:label]) <: Array
                    i += 1
                    lab = series[:label][i]
                else
                    lab = series[:label]
                end
                GR.settextalign(GR.TEXT_HALIGN_LEFT, GR.TEXT_VALIGN_HALF)
                gr_set_textcolor(sp[:legendfontcolor])
                gr_text(xpos, ypos, string(lab))
                ypos -= dy
            end
        end
        GR.selntran(1)
        GR.restorestate()
    end

    # add annotations
    GR.savestate()
    # update the bounding window
    if ispolar(sp)
        gr_set_viewport_polar()
    else
        xmin, xmax, ymin, ymax = data_lims
        if xmax > xmin && ymax > ymin
            GR.setwindow(xmin, xmax, ymin, ymax)
        end
    end
    for ann in sp[:annotations]
        x, y, val = locate_annotation(sp, ann...)
        x, y = if is3d(sp)
            # GR.wc3towc(x, y, z)
        else
            GR.wctondc(x, y)
        end
        gr_set_font(val.font)
        gr_text(x, y, val.str)
    end
    GR.restorestate()
end

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

const _gr_mimeformats = Dict(
    "application/pdf"         => "pdf",
    "image/png"               => "png",
    "application/postscript"  => "ps",
    "image/svg+xml"           => "svg",
)

const _gr_wstype_default = @static if Sys.islinux()
    "x11"
    # "cairox11"
elseif Sys.isapple()
    "quartz"
else
    "use_default"
end

const _gr_wstype = Ref(get(ENV, "GKSwstype", _gr_wstype_default))
gr_set_output(wstype::String) = (_gr_wstype[] = wstype)

for (mime, fmt) in _gr_mimeformats
    @eval function _show(io::IO, ::MIME{Symbol($mime)}, plt::Plot{GRBackend})
        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})
    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
        if _gr_wstype[] != "use_default"
            ENV["GKSwstype"] = _gr_wstype[]
        end
        gr_display(plt)
    end
end

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


text:
module Revtok
export tokenize, detokenize, buildvocab, segment

include("tokenizer.jl")
include("subwords.jl")

end

ext: jl


text:
# Some Boiler-plate code for preconditioning
#
# Meaning of P:
#    P ≈ ∇²E, so the preconditioned gradient is P^{-1} ∇E
#    P can be an arbitrary type but at a minimum it MUST provide:
#        ldiv!(x, P, b) ->  x = P \ b
#        dot(x, P, b)  ->  x' P b
#
#    If `dot` is not provided, then Optim.jl will try to define it via
#        dot(x, P, b) = dot(x, mul!(similar(x), P, y))
#
#    finally the preconditioner can be updated after each x-update using
#        precondprep!
#    but this is passed as an argument at the moment!
#

#####################################################
#  [0] Defaults and aliases for easier reading of the code
#      these can also be over-written if necessary.

if VERSION < v"1.4.0-DEV.92"
# an inner product w.r.t. a metric P (=preconditioner)
    dot(x, P, y) = dot(x, mul!(similar(x), P, y))
end

# default preconditioner update
precondprep!(P, x) = nothing

#####################################################
#  [1] Empty preconditioner = Identity
#

# out =  P^{-1} * A
ldiv!(out, ::Nothing, A) = copyto!(out, A)
# A' * P B
dot(A, ::Nothing, B) = dot(A, B)

#####################################################
#  [2] Diagonal preconditioner
#      P = Diag(d)
#      unfortunately, Base does not implement
#      ldiv!(a, P, b) or mul! for this type, so we do it by hand
#      TODO: maybe implement this in Base
ldiv!(out::Array, P::Diagonal, A::Array) = copyto!(out, A ./ P.diag)
dot(A::Array, P::Diagonal, B::Array) = dot(A, P.diag .* B)

#####################################################
#  [3] Inverse Diagonal preconditioner
#      here, P is stored by the entries of its inverse
#      TODO: maybe implement this in Base?
mutable struct InverseDiagonal
   diag
end
ldiv!(out::Array, P::InverseDiagonal, A::Array) = copyto!(out, A .* P.diag)
dot(A::Array, P::InverseDiagonal, B::Vector) = dot(A, B ./ P.diag)

#####################################################
#  [4] Matrix Preconditioner
#  the assumption here is that P is given by its inverse, which is typical
#     > ldiv! is about to be moved to Base, so we need a temporary hack
#     > mul! is already in Base, which defines `dot`
#  nothing to do!
ldiv!(x, P::AbstractMatrix, b) = copyto!(x, P \ b)


text:
# A Benchmark object represents a single "benchmarkable" that has an associated
# name for logging purposes.
#
# Fields:
#
#    name::String: The name of the benchmark.
#
#    f!::Function: A function that implements the "benchmarkable" protocol.

immutable Benchmark
    name::Compat.String
    f!::Function
end

# A BenchmarkSuite object represents a suite of related benchmarks.
#
# Fields:
#
#     name::String: The name of the suite of benchmarks.
#
#     benchmarks::Vector{Benchmark}: A vector of the benchmarks in the suite.

immutable BenchmarkSuite
    name::Compat.String
    benchmarks::Vector{Benchmark}
end


text:
#Use only via ZMQSer package

function launch_client(endpoint::ASCIIString)
    zctx = ZMQContext()
    requester = ZMQSocket(zctx, ZMQ_REQ)
    connect(requester, endpoint)
    return zctx, requester
end
launch_client() = launch_client("tcp://localhost:5555")

# Run svd remotely:
#   A = randn(3,5)
#   U, S, V = zmqcall(requester, :svd, A)
function zmqcall(requester::ZMQSocket, func::Symbol, args...)
    ex = expr(:call, {func, args...})
    zmq_serialize(requester, ex)
    zmq_deserialize(requester)
end

# Remotely parse a string and execute it, e.g.,
#    str = "x = randn(7); sort(x)"
#    y = zmqparse(requester, str)
# Within Julia it may be better to use a quote block, but this will
# surely be the easy way from other languages
function zmqparse(requester::ZMQSocket, str::ASCIIString)
    ex = :(parse_eval($str))
    zmq_serialize(requester, ex)
    zmq_deserialize(requester)
end

# Set a variable in the remote session
function zmqsetvar(requester::ZMQSocket, var::Symbol, val)
    ex = :($var = $val)
    zmq_serialize(requester, ex)
    zmq_deserialize(requester)
end

# Get a variable in the remote session
function zmqgetvar(requester::ZMQSocket, var::Symbol)
    ex = :($var)
    zmq_serialize(requester, ex)
    zmq_deserialize(requester)
end


text:
using DashBootstrapComponents, DashHtmlComponents

badges = html_span([
    dbc_badge(
        "Primary",
        color = "white",
        text_color = "primary",
        className = "border me-1",
    ),
    dbc_badge(
        "Secondary",
        color = "white",
        text_color = "secondary",
        className = "border me-1",
    ),
    dbc_badge(
        "Success",
        color = "white",
        text_color = "success",
        className = "border me-1",
    ),
    dbc_badge(
        "Warning",
        color = "white",
        text_color = "warning",
        className = "border me-1",
    ),
    dbc_badge("Danger", color = "white", text_color = "danger", className = "border me-1"),
    dbc_badge("Info", color = "white", text_color = "info", className = "border me-1"),
    dbc_badge("Dark", color = "white", text_color = "dark", className = "border me-1"),
    dbc_badge("Black", color = "white", text_color = "black", className = "border me-1"),
    dbc_badge("Muted", color = "white", text_color = "muted", className = "border me-1"),
    dbc_badge(
        "Light",
        color = "secondary",
        text_color = "light",
        className = "border me-1",
    ),
    dbc_badge("White", color = "secondary", text_color = "white", className = "border"),
])

max_stars_repo_name: glsdown/dash-bootstrap-components

max_issues_repo_name: glsdown/dash-bootstrap-components

max_forks_repo_name: glsdown/dash-bootstrap-components


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using JuliaFEM, LinearAlgebra, Statistics, Test

mesh = aster_read_mesh(@__DIR__()*"/testdata/primitives.med", "CYLINDER_20_TET4")
problem_elements = create_elements(mesh, "CYLINDER")
update!(problem_elements, "thermal conductivity", 200.0)
outer_elements = create_elements(mesh, "FACE2", "FACE3")
update!(outer_elements, "external temperature", 20.0)
update!(outer_elements, "heat transfer coefficient", 1.0)
#midline = Problem(Heat, "midline of rod", 1)
#midline.elements = create_elements(mesh, "INNER_LINE")
boundary_elements = create_elements(mesh, "FACE1")
update!(boundary_elements, "temperature 1", 100.0)

problem = Problem(Heat, "rod of length 20", 1)
add_elements!(problem, problem_elements)
outer = Problem(Heat, "outer surface", 1)
add_elements!(outer, outer_elements)
boundary = Problem(Dirichlet, "homogeneous dirichlet boundary", 1, "temperature")
add_elements!(boundary, boundary_elements)

analysis = Analysis(Linear)
add_problems!(analysis, problem, outer, boundary)
run!(analysis)

# Analytical solution
L = 20
k = 200.0
Tu = 20.0
h = 1.0
P = 2*pi
A = pi
α = h
β = sqrt((h*P)/(k*A))
T0 = 100.0
C = [1.0 1.0; (α+k*β)*exp(β*L) (α-k*β)*exp(-β*L)] \ [T0-Tu, 0]

T_diff = []
for x in range(0, stop=20)
    T_FEM = problem("temperature", [x, 0.0, 0.0], 0.0)[1]
    T_ACC = dot(C, [exp(β*x), exp(-β*x)]) + Tu
    push!(T_diff, norm(T_FEM - T_ACC))
end
@debug("mean diff on heat problem", mean(T_diff))

@test mean(T_diff) < 1.2 # mean diff = 1.14


text:
using Test
using TestExtras
using Random
using TensorKit
using SUNRepresentations: SUNIrrep
using Combinatorics
using TensorKit
using TensorKit: ProductSector, fusiontensor
using TensorOperations
TensorOperations.disable_cache() # avoids memory overflow during CI?
using Base.Iterators: take, product
using SparseArrayKit: SparseArray
import LinearAlgebra

const TK = TensorKit

Random.seed!(1234)

smallset(::Type{I}) where {I<:Sector} = take(values(I), 5)
function smallset(::Type{ProductSector{Tuple{I1,I2}}}) where {I1,I2}
    iter = product(smallset(I1), smallset(I2))
    s = collect(i ⊠ j for (i,j) in iter if dim(i)*dim(j) <= 6)
    return length(s) > 6 ? rand(s, 6) : s
end
function smallset(::Type{ProductSector{Tuple{I1,I2,I3}}}) where {I1,I2,I3}
    iter = product(smallset(I1), smallset(I2), smallset(I3))
    s = collect(i ⊠ j ⊠ k for (i,j,k) in iter if dim(i)*dim(j)*dim(k) <= 6)
    return length(s) > 6 ? rand(s, 6) : s
end
function randsector(::Type{I}) where {I<:Sector}
    s = collect(smallset(I))
    a = rand(s)
    while a == one(a) # don't use trivial label
        a = rand(s)
    end
    return a
end

Ti = time()
module GenericTests
    using Test
    using TestExtras
    using Random
    using SUNRepresentations
    include("generic.jl")
end

sectorlist = (SUNIrrep{3}, SUNIrrep{4}, SUNIrrep{5}, SUNIrrep{3} ⊠ SUNIrrep{3})
include("sectors.jl")
sectorlist = (SUNIrrep{3}, SUNIrrep{4}, SUNIrrep{5})
include("fusiontrees.jl")

Tf = time()
printstyled("Finished all tests in ",
            string(round((Tf-Ti)/60; sigdigits=3)),
            " minutes."; bold = true, color = Base.info_color())
println()


text:
# using DataFrames
# using Serialization
# using BioSequences
# using MLJ
using MLJGLMInterface
using StableRNGs
# using CSV

#%%

include("mlj_functions.jl")

#%%

do_GLM = true
do_evaluate = true
do_shap = true
do_bases_included_accuracy = true
make_figures = true
save_figures = true

do_GLM = false
do_evaluate = false
do_shap = false
# do_bases_included_accuracy = false
make_figures = false
save_figures = false

#%%

# x = x

# filename_csv = "df__N_reads__100000.csv"
# df = DataFrame(CSV.File(filename_csv, drop = [1]));

filename = "./data/df.data"
N_rows = 1_000_000
# N_rows = 10_000

X, y = get_Xy(filename, N_rows);
train, test = partition(eachindex(y), 0.75; shuffle = true, rng = StableRNG(123));

X_train = selectrows(X, train);
y_train = selectrows(y, train);
X_test = selectrows(X, test);
y_test = selectrows(y, test);

#%%

# x = x

#%%

if make_figures
    println("Signal proportion, base-stratified")
    get_base_stratified_signal_proportion(X[test, :], y_test, mean, Float64)
    get_base_stratified_signal_proportion(X[test, :], y_test, sum, Int64)
end

#%%

if make_figures
    println("Base count fractions")
    mask_signal = (y .== 1)
    base_counts_signal = get_base_counts_pr_position(X[mask_signal, :], normalise = true)
    base_counts_background =
        get_base_counts_pr_position(X[.!mask_signal, :], normalise = true)
    f_base_fraction = make_base_fraction_plot(base_counts_signal, base_counts_background)
    if save_figures
        save("./figures/base_fraction__$(N_rows).pdf", f_base_fraction)
    end
end

#%%

# get_matching_models(X, y)

# #%%
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
#                     ██      ██████
#                     ██      ██   ██
#                     ██      ██████
#                     ██      ██   ██
#                     ███████ ██   ██
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

LR = @load LogisticClassifier pkg = MLJLinearModels

pipe_LR = Pipeline(OneHotEncoder, LR(penalty = :none, fit_intercept = false), cache = false);

filename_LR = "./data/mach_LR__$(N_rows).jlso"

if isfile(filename_LR)
    println("Loading LR: $filename_LR")
    mach_LR = machine(filename_LR)
else
    println("Fitting LR: $filename_LR")
    mach_LR = machine(pipe_LR, X_train, y_train)
    fit!(mach_LR, verbosity = 0)
    MLJ.save(filename_LR, mach_LR)

    df_LR_long = get_df_LR_long(mach_LR)
    df_LR_wide = get_df_LR_wide(df_LR_long)
    if make_figures
        f_LR_fit_coef = plot_LR_fit_coefficients(df_LR_wide)
        if save_figures
            save("./figures/LR_fit_coefficient__$(N_rows).pdf", f_LR_fit_coef)
        end
    end

end

yhat_LR = predict(mach_LR, X_test);
acc_LR = accuracy(yhat_LR, y_test)
println("Accuracy, LR, ", round(acc_LR * 100, digits = 2), "%")
confusion_matrix(yhat_LR, y_test)

if make_figures
    f_score_density_LR = plot_score_density(yhat_LR, y_test, "Density of scores for LR")
    if save_figures
        save("./figures/score_density__LR__$(N_rows).pdf", f_score_density_LR)
    end
end

println("AUC, base-stratified, Logistic Regression:")
base_strat_auc_LR = get_base_stratified_measure(X_test, yhat_LR, y_test, area_under_curve)

if do_evaluate
    eval_LR = evaluate(
        pipe_LR,
        X,
        y,
        resampling = CV(nfolds = 5, shuffle = true, rng = 1234),
        measures = [LogLoss(), Accuracy()],
        verbosity = 0,
    )

    print_performance(eval_LR, "Logistic Regression")
end

#%%
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
#                      ██████  ██      ███    ███
#                     ██       ██      ████  ████
#                     ██   ███ ██      ██ ████ ██
#                     ██    ██ ██      ██  ██  ██
#                      ██████  ███████ ██      ██
# # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

if do_GLM

    pipe_GLM = Pipeline(
        OneHotEncoder(drop_last = true),
        x -> table(Matrix(x)),
        LinearBinaryClassifier(fit_intercept = true),
        cache = false,
    )

    filename_GLM = "./data/mach_GLM__$(N_rows).jlso"

    if isfile(filename_GLM)
        println("Loading GLM: $filename_GLM")
        mach_GLM = machine(filename_GLM)
    else
        println("Fitting GLM: $filename_GLM")
        mach_GLM = machine(pipe_GLM, X_train, y_train)
        fit!(mach_GLM, verbosity = 0)
        MLJ.save(filename_GLM, mach_GLM)

        fitresult = get_glm_fitresult(mach_GLM)

        fitted_params(mach_GLM).linear_binary_classifier.coef
        r = report(mach_GLM).linear_binary_classifier
        r.stderror
        # r.deviance
        # r.dof_residual
        # r.vcov
    end

    yhat_GLM = predict(mach_GLM, X_test)

    acc_glm = accuracy(yhat_GLM, y_test)
    println("Accuracy, GLM, ", round(100 * acc_glm, digits = 2), "%")
    confusion_matrix(yhat_GLM, y_test)

end

#%% LGB categorical
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
#                     ██       ██████  ██████   ██████
#                     ██      ██       ██   ██ ██
#                     ██      ██   ███ ██████  ██
#                     ██      ██    ██ ██   ██ ██
#                     ███████  ██████  ██████   ██████
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

LGB = @load LGBMClassifier verbosity = 0
categorical_columns = collect(1:size(X)[2]);

pipe_LGB = Pipeline(
    x -> convert_type(x, Float64),
    LGB(
        objective = "binary",
        num_iterations = 100,
        learning_rate = 0.1,
        feature_fraction = 0.8,
        bagging_fraction = 0.9,
        bagging_freq = 1,
        num_leaves = 1000,
        metric = ["auc", "binary_logloss"],
        categorical_feature = copy(categorical_columns),
    ),
    cache = false,
);

filename_LGB = "./data/mach_LGB__$(N_rows).jlso"

if isfile(filename_LGB)
    println("Loading LGB: $filename_LGB")
    mach_LGB = machine(filename_LGB)
else
    println("Fitting LGB: $filename_LGB")
    mach_LGB = machine(pipe_LGB, X_train, y_train)
    fit!(mach_LGB, verbosity = 0)
    MLJ.save(filename_LGB, mach_LGB)
end

yhat_LGB = predict(mach_LGB, X_test);
acc_LGB = accuracy(yhat_LGB, y_test);
println("Accuracy, LGB, ", round(100 * acc_LGB, digits = 2), "%")
confusion_matrix(yhat_LGB, y_test)

if make_figures
    f_score_density_LGB = plot_score_density(yhat_LGB, y_test, "Density of scores for LGB")
    if save_figures
        save("./figures/score_density__LGB__$(N_rows).pdf", f_score_density_LGB)
    end
end

println("AUC, base-stratified, LightGBM:")
base_strat_auc_LGB = get_base_stratified_measure(X_test, yhat_LGB, y_test, area_under_curve)

if do_evaluate
    eval_LGB = evaluate(
        pipe_LGB,
        X,
        y,
        resampling = CV(nfolds = 5, shuffle = true, rng = 1234),
        measures = [LogLoss(), Accuracy()],
        verbosity = 0,
    )
    print_performance(eval_LGB, "LightGBM")
end

#%%
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
#                     ██   ██ ██████   ██████
#                     ██   ██ ██   ██ ██    ██
#                     ███████ ██████  ██    ██
#                     ██   ██ ██      ██    ██
#                     ██   ██ ██       ██████
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

# atom = LGB(); #initialised a model with default params
# ensemble = @pipeline(OneHotEncoder, atom);
# mach = machine(pipe_LGB2, X, y)
# boostrange = range(ensemble, :(LGBm_classifier.num_iterations), lower=2, upper=500)
# curve = learning_curve!(mach, resampling=CV(nfolds=5),
#                         range=boostrange,
#                         resolution=100,
#                         measure=LogLoss())

# mach = fit!(machine(OneHotEncoder(), X))
# Xt = MLJ.transform(mach, X)

# LGB = LGB(); #initialised a model with default params
# LGBm = machine(LGB, Xt, y)
# boostrange = range(LGB, :num_iterations, lower=2, upper=500)
# curve = learning_curve!(LGBm, resampling=CV(nfolds=5),
#                         range=boostrange, resolution=100,
#                         measure=LogLoss())

# x = x

#%%
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
#                     ██████   ██████   ██████
#                     ██   ██ ██    ██ ██
#                     ██████  ██    ██ ██
#                     ██   ██ ██    ██ ██
#                     ██   ██  ██████   ██████
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

if make_figures

    y_hats = (yhat_LR, yhat_LGB)
    names_roc = ("Logistic Regression", "LightGBM")
    f_roc = plot_roc(y_hats, y_test, names_roc)
    if save_figures
        save("./figures/ROC__$(N_rows).pdf", f_roc)
    end

end

# x=x

#%%
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
#                     ███████ ██   ██  █████  ██████
#                     ██      ██   ██ ██   ██ ██   ██
#                     ███████ ███████ ███████ ██████
#                          ██ ██   ██ ██   ██ ██
#                     ███████ ██   ██ ██   ██ ██
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

if do_shap

    CairoMakie.activate!()

    data_shap_LR = shap_compute(mach = mach_LR, X = X_test, N = 100)

    if make_figures
        plot_shap_global(data_shap_LR, savefig = true, filename = "shap_global_LR")
        plot_shap_variables(
            data_shap_LR,
            savefig = save_figures,
            filename = "shap_feature_LR",
        )
    end

    data_shap_LGB = shap_compute(mach = mach_LGB, X = X_test, N = 100)
    if make_figures
        plot_shap_global(
            data_shap_LGB,
            savefig = save_figures,
            filename = "shap_global_LGB",
        )
        plot_shap_variables(
            data_shap_LGB,
            savefig = save_figures,
            filename = "shap_feature_LGB",
        )
    end

end

#%%
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
#                     ██████   █████  ███████ ███████ ███████
#                     ██   ██ ██   ██ ██      ██      ██
#                     ██████  ███████ ███████ █████   ███████
#                     ██   ██ ██   ██      ██ ██           ██
#                     ██████  ██   ██ ███████ ███████ ███████
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

if do_bases_included_accuracy

    filename_accuracies = "./data/accuracies__$(N_rows).data"

    if isfile(filename_accuracies)
        accuracies = deserialize(filename_accuracies).accuracies
        accuracies_centered = deserialize(filename_accuracies).accuracies_centered
    else
        accuracies = get_accuracies_pr_base(X, y, train, test)
        accuracies_centered = get_accuracies_pr_base_centered(
            X,
            y,
            train,
            test;
            add_analytical = true,
            add_custom = true,
        )

        serialize(
            filename_accuracies,
            (accuracies = accuracies, accuracies_centered = accuracies_centered),
        )

    end

    if make_figures
        # GLMakie.activate!()
        CairoMakie.activate!()

        f_acc = plot_accuracy_function_of_bases(accuracies)
        if save_figures
            save("./figures/accuracies_base_dependent__$(N_rows).pdf", f_acc)
        end

        f_acc_centered = plot_accuracy_function_of_bases_centered(accuracies_centered)
        if save_figures
            save(
                "./figures/accuracies_base_dependent_centered__$(N_rows).pdf",
                f_acc_centered,
            )
        end
    end

    #%%

end

#%%

max_stars_repo_name: ChristianMichelsen/reference-free

max_issues_repo_name: ChristianMichelsen/reference-free

max_forks_repo_name: ChristianMichelsen/reference-free


text:
using Revise
using AdFem
using LinearAlgebra
using PyPlot

m = 50
n = 50
h = 1/n

bdedge = bcedge("top", m, n, h)
bdnode = bcnode("left|right|bottom", m, n, h)

F1 = eval_f_on_gauss_pts((x,y)->3.0, m, n, h)
F2 = eval_f_on_gauss_pts((x,y)->-1.0, m, n, h)
F = compute_fem_source_term(F1, F2, m, n, h)

t1 = eval_f_on_boundary_edge((x,y)->-x-y, bdedge, m, n, h)
t2 = eval_f_on_boundary_edge((x,y)->2y, bdedge, m, n, h)
T = compute_fem_traction_term([t1 t2], bdedge, m, n, h)

D = constant(diagm(0=>[1,1,0.5]))
K = compute_fem_stiffness_matrix(D, m, n, h)
rhs = T - F
bdval = [eval_f_on_boundary_node((x,y)->x^2+y^2, bdnode, m, n, h);
        eval_f_on_boundary_node((x,y)->x^2-y^2, bdnode, m, n, h)]
DOF = [bdnode;bdnode .+ (m+1)*(n+1)]
K, rhs = impose_Dirichlet_boundary_conditions(K, rhs, DOF, bdval)
u = K\rhs
sess = Session(); init(sess)
S = run(sess, u)

figure(figsize=[10,4])
subplot(121)
visualize_scalar_on_fem_points(S[1:(m+1)*(n+1)], m, n, h)
subplot(122)
visualize_scalar_on_fem_points(S[(m+1)*(n+1)+1:2(m+1)*(n+1)], m, n, h)
savefig("numerical.png")

figure(figsize=[10,4])
x = LinRange(0, 1, 50)
X = zeros(50, 50)
Y = zeros(50, 50)
for i = 1:50
    for j = 1:50
        X[i,j] = x[i]
        Y[i,j] = x[j]
    end
end
subplot(121)
XY = fem_nodes(m, n, h)
X = XY[:,1]
Y = XY[:,2]
U = (@. X^2+Y^2)
V = (@. X^2-Y^2)
subplot(121)
visualize_scalar_on_fem_points(U, m, n, h)
subplot(122)
visualize_scalar_on_fem_points(V, m, n, h)
savefig("exact.png")

figure(figsize=[10,4])
subplot(121)
visualize_scalar_on_fem_points(U-S[1:(m+1)*(n+1)], m, n, h)
subplot(122)
visualize_scalar_on_fem_points(V-S[(m+1)*(n+1)+1:2(m+1)*(n+1)], m, n, h)
savefig("difference.png")


text:
using DepositionMTK
using Documenter

DocMeta.setdocmeta!(DepositionMTK, :DocTestSetup, :(using DepositionMTK); recursive=true)

makedocs(;
    modules=[DepositionMTK],
    authors="Chris Tessum <ctessum@gmail.com> and contributors",
    repo="https://github.com/ctessum/DepositionMTK.jl/blob/{commit}{path}#{line}",
    sitename="DepositionMTK.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://ctessum.github.io/DepositionMTK.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/ctessum/DepositionMTK.jl",
)


text:
using DataFrames
using StatsBase
using Statistics

# Optional: CategoricalArrays
# using CategoricalArrays

# Optional: Embedded R within Julia
# using RCall
# @rlibrary ggplot2

# Optional: Revise.jl tracks files and updates methods upon file writes
# using Revise
# includet(YourJuliaSourceFile.jl)

# ------------------------------------------------------------------------------
# Arrays
# ------------------------------------------------------------------------------
# Arrays are mutable and passed by reference.
Array{Char}(undef, 2, 3, 4) # uninitialized 2x3x4 array of Chars
a = [1,2]
a[1] # first element of array (1-based indexing)
a[end] # get last element of array
firstindex(a)
lastindex(a)

ndims(a) # number of dimensions
eltype(a) # type of elements in array
length(a) # number of elements in array
size(a) # dimensions of array as a tuple

zeros(2) # Float64 zeros
ones(2) # Float64 ones
trues(2)
falses(2)

r = range(0, stop = 1, length = 11)
collect(r) # converts iterator to array
collect(1:2:9)

reshape(1:6, 2, 3) # reshapes array into a 2x3 matrix, filled column-wise
fill("a", 2, 2)

[x*y for x in 1:2, y in 1:3] # list comprehension

# start:step:end
for i in 0:0.1:1
    println(i)
end

# ------------------------------------------------------------------------------
# DataFrames.jl package
# ------------------------------------------------------------------------------
df = DataFrame(A = 1:10, B = rand(10))
describe(df)
names(df)
nrow(df)
ncol(df)
push!(df, (50, 40.1))
rename!(df, :A => :A_new)
select!(df, Not(:A_new)) # remove column A
delete!(df, 1:4) # remove rows 1 through 4

df = DataFrame(A = 1:10, B = rand(10))
sort!(df, [order(:B), order(:A)]) # sort by A, then by B

filter(row -> row.A > 3, df)
df[df.A .> df.B, :]
df[(df.A .> 3) .& (df.A .< 5), :]
[mean(col) for col in eachcol(df)]
mean.(eachcol(df))

# ------------------------------------------------------------------------------
# Dictionaries
# ------------------------------------------------------------------------------
Dict{Char, Int}()
D = Dict("a"=>1, "b"=>2)
haskey(D, "b")
keys(D)
values(D)
delete!(D, "b")
get(D, "c", "default_value_if_no_c")

# ------------------------------------------------------------------------------
# Strings
# ------------------------------------------------------------------------------
"Hello" * ", world!" # string concatenation
string("Answer to the universe is ", 42)
occursin("cd", "abcd")
chop("testing") # removes last character of string and returns a SubString

# Regular expressions
r = r"A|B"
occursin(r, "CD") # false
match(r, "ACBD")

# ------------------------------------------------------------------------------
# Types
# ------------------------------------------------------------------------------
# var::TypeName # type annotation

# Type declaration; structs are immutable!
struct Person
    name::String
    age::UInt16
end

bob = Person("Bob", 30) # type instantiation

# Mutable type declaration
mutable struct Programmer
    name::String
    age::UInt16
    drinks_coffee::Bool
end

# Subtype declaration
abstract type SuperType end
struct SubType <: SuperType
    name::String
end

# Parametric types
struct Point{T <: Real}
    x::T
    y::T
end
p = Point{Float64}(1,2)

Union{Int, String} # Union types

fieldnames(Person)
Person.types # field types

('a', false)::Tuple{Char, Bool} # tuple type assertion
namedtuple = (p=1, q=2, r=3)
typeof(namedtuple)


text:
module MYGENrun2

using Modia
@usingModiaPlot
using StaticArrays

include("../src/Utils.jl")
include("../src/SIM8.jl")

# logData = true

GenDataMap = ReadDataSet("GENdata.yaml")
ZoneDataMap = ReadDataSet("ZONEdata.yaml")
ST8CMap = ReadDataSet("ST8Cdata.yaml")

instModel = Main.MYGEN2.instModel

include("GENinit01.start.jl")

function mySimulation(;SW1="TESTSW1", reuse=false)
    printstyled("\nNEW simulation\n",color=:light_blue)

    simulate!(instModel,
#         IDA(),
        CVODE_BDF(),
        startTime=-0.05,
        stopTime = 1.0,
        interval=2.123e-4,
        tolerance=1e-6,
        log=false,
        logStates=false,
        logEvents=false,
        logTiming = false,
        logParameters    = false,
        logEvaluatedParameters   = false,

        merge=Map(
            transformer = ZoneDataMap["setGEN"],
            meterGEN    = ZoneDataMap["setGEN"],
            synUnit     = Map(exciter = ST8CMap[SW1], SM = GenDataMap["NONsat"]),
            varvolts    = Map(Tstart=1e-3),     # Event muss T>0.0

            # das GEHT NICHT !! J = Map(w=Var(init=38.0)),
            #  instModel.parameters[:J][:w]=39.0
        ),
    )
#     printstyled("total time for simulation: ", color = :yellow)

#      plot(instModel, [("synUnit.SM.voltage"), ( "synUnit.exciter.y"), ( "synUnit.SM.ifd")], figure=1)

     plot(instModel, [( "synUnit.exciter.y"),("synUnit.phiDEG")], prefix=ST8CMap[SW1][:Project] *"----", reuse=reuse, figure=1)
#      plot(instModel, [ ("synUnit.SM.tau")], figure=1)
#      plot(instModel, [("synUnit.SM.i"),("GridVolts.volts")], figure=2)
#       plot(instModel, [ ("GridVolts.amps")], figure=2)

end

# open("LOG.txt", "a") do io
#     redirect_stdout(io) do
        mySimulation()
        mySimulation(SW1="TESTSW0", reuse=true)
#     end
# end

Ergebnis(instModel, [
    "VEXC",
    "LOAD",
    "meterGRID.power",
    "meterGEN.power",
    "synUnit.SM.voltage",
    "synUnit.SM.ifd",
    "synUnit.exciter.y",
    "synUnit.phiDEG",
    "J.w"])

printstyled(length(instModel.result_x[1,:]),"\n", color=:red)

# Speichern(instModel, ["time","synUnit.exciter.y", "synUnit.SM.ifd", "synUnit.SM.tau", "meterGEN.voltage", "meterGEN.current", "meterGRID.voltage", "synUnit.SM.phiKR"])

end

max_stars_repo_name: johhell/EPSphasor

max_issues_repo_name: johhell/EPSphasor

max_forks_repo_name: johhell/EPSphasor


text:
function TriMesh_Interpolate(x,eConn,eAdjacency,fNodes,xInterp)
#= """
  TriMesh_Interpolate - interpolates finite element functions at given points.

  This routine exploits the fact that the incoming mesh is Delauney.  If it
  isn't, it is possible to land in an infinite loop.  See the functions
  TriMesh_TestDelauney and TriMesh_ElementAdjacency.

  Author: Jeff Borggaard, Virginia Tech
          part of FEMfunctions.jl

  Licensing:
     This code is distributed under the MIT license.

  Usage:
  ```julia
    fValues = TriMesh_Interpolate(x,eConn,eAdjacency,fNodes,xInterp)
  ```

  Arguments:
  - `x`: nodal coordinates of the triangular mesh
  - `eConn`: element connectivity
  - `eAdjacency`: element adjacency
  - `fNodes`: an array of nodal values for m different functions
     (size(fNodes) = nNodes x nFunctions)
  - `xInterp`: a set of interpolation points (must be in the interior)

  Output argument:
  - `fValues`: an array of interpolated values
""" =#

  #  Get problem dimensions (we don't check the inputs)
  nNodes = size(x,1)
  nElements = size(eConn,1)
  nPoints = size(xInterp,1)
  nFunctions = size(fNodes,2)

  elementList = Array{Int64,1}(undef,nPoints)
  fValues = zeros(Float64,nPoints,nFunctions)

  for i=1:nPoints
    elementList[i], baryCoord1, baryCoord2 =
            TriMesh_Search(x,eConn[:,1:3],eAdjacency,xInterp[i,:])

    nodesLocal = eConn[ elementList[i], : ]
    xLocal = x[nodesLocal,:]
    xg,wg,ϕ,ϕ_x,ϕ_y = twodShape(xLocal,[baryCoord1],[baryCoord2],[1.0])

    for j=1:nFunctions
      tmp = ϕ*fNodes[nodesLocal,j]
      fValues[i,j] = tmp[1]
    end
  end

  return fValues, elementList
end


text:
"""
    LmResp

Encapsulates the response for a linear model

# Members

- `mu`: current value of the mean response vector or fitted value
- `offset`: optional offset added to the linear predictor to form `mu`
- `wts`: optional vector of prior frequency (a.k.a. case) weights for observations
- `y`: observed response vector

Either or both `offset` and `wts` may be of length 0
"""
mutable struct LmResp{V<:FPVector} <: ModResp  # response in a linear model
    mu::V                                  # mean response
    offset::V                              # offset added to linear predictor (may have length 0)
    wts::V                                 # prior weights (may have length 0)
    y::V                                   # response
    function LmResp{V}(mu::V, off::V, wts::V, y::V) where V
        n = length(y)
        length(mu) == n || error("mismatched lengths of mu and y")
        ll = length(off)
        ll == 0 || ll == n || error("length of offset is $ll, must be $n or 0")
        ll = length(wts)
        ll == 0 || ll == n || error("length of wts is $ll, must be $n or 0")
        new{V}(mu, off, wts, y)
    end
end

LmResp(y::FPVector, wts::FPVector=similar(y, 0)) =
    LmResp{typeof(y)}(fill!(similar(y), 0), similar(y, 0), wts, y)

LmResp(y::AbstractVector{<:Real}, wts::AbstractVector{<:Real}=similar(y, 0)) =
    LmResp(float(y), float(wts))

function updateμ!(r::LmResp{V}, linPr::V) where V<:FPVector
    n = length(linPr)
    length(r.y) == n || error("length(linPr) is $n, should be $(length(r.y))")
    length(r.offset) == 0 ? copyto!(r.mu, linPr) : broadcast!(+, r.mu, linPr, r.offset)
    deviance(r)
end

updateμ!(r::LmResp{V}, linPr) where {V<:FPVector} = updateμ!(r, convert(V, vec(linPr)))

function deviance(r::LmResp)
    y = r.y
    mu = r.mu
    wts = r.wts
    v = zero(eltype(y)) + zero(eltype(y)) * zero(eltype(wts))
    if isempty(wts)
        @inbounds @simd for i = eachindex(y,mu)
            v += abs2(y[i] - mu[i])
        end
    else
        @inbounds @simd for i = eachindex(y,mu,wts)
            v += abs2(y[i] - mu[i])*wts[i]
        end
    end
    v
end

function nulldeviance(r::LmResp)
    y = r.y
    wts = r.wts
    if isempty(wts)
        m = mean(y)
    else
        m = mean(r.y, weights(r.wts))
    end

    v = zero(eltype(y))*zero(eltype(wts))
    if isempty(wts)
        @inbounds @simd for i = 1:length(y)
            v += abs2(y[i] - m)
        end
    else
        @inbounds @simd for i = 1:length(y)
            v += abs2(y[i] - m)*wts[i]
        end
    end
    v
end

function loglikelihood(r::LmResp)
    n = isempty(r.wts) ? length(r.y) : sum(r.wts)
    -n/2 * (log(2π * deviance(r)/n) + 1)
end

function nullloglikelihood(r::LmResp)
    n = isempty(r.wts) ? length(r.y) : sum(r.wts)
    -n/2 * (log(2π * nulldeviance(r)/n) + 1)
end

residuals(r::LmResp) = r.y - r.mu

"""
    LinearModel

A combination of a [`LmResp`](@ref) and a [`LinPred`](@ref)

# Members

- `rr`: a `LmResp` object
- `pp`: a `LinPred` object
"""
struct LinearModel{L<:LmResp,T<:LinPred} <: LinPredModel
    rr::L
    pp::T
end

LinearAlgebra.cholesky(x::LinearModel) = cholesky(x.pp)

function StatsBase.fit!(obj::LinearModel)
    if isempty(obj.rr.wts)
        delbeta!(obj.pp, obj.rr.y)
    else
        delbeta!(obj.pp, obj.rr.y, obj.rr.wts)
    end
    installbeta!(obj.pp)
    updateμ!(obj.rr, linpred(obj.pp, zero(eltype(obj.rr.y))))
    return obj
end

function fit(::Type{LinearModel}, X::AbstractMatrix{<:Real}, y::AbstractVector{<:Real},
             allowrankdeficient::Bool=false; wts::AbstractVector{<:Real}=similar(y, 0))
    fit!(LinearModel(LmResp(y, wts), cholpred(X, allowrankdeficient)))
end

"""
    lm(X, y, allowrankdeficient::Bool=false; wts=similar(y, 0))

An alias for `fit(LinearModel, X, y, allowrankdeficient)`

The arguments `X` and `y` can be a `Matrix` and a `Vector` or a `Formula` and a `DataFrame`.

The keyword argument `wts` can be a `Vector` specifying frequency weights for observations.
Such weights are equivalent to repeating each observation a number of times equal
to its weight. Do note that this interpretation gives equal point estimates but
different standard errors from analytical (a.k.a. inverse variance) weights and
from probability (a.k.a. sampling) weights which are the default in some other
software.
"""
lm(X, y, allowrankdeficient::Bool=false; kwargs...) =
    fit(LinearModel, X, y, allowrankdeficient; kwargs...)

dof(x::LinearModel) = length(coef(x)) + 1

dof(obj::LinearModel{<:LmResp,<:DensePredChol{<:Real,<:CholeskyPivoted}}) = obj.pp.chol.rank + 1

"""
    deviance(obj::LinearModel)

For linear models, the deviance is equal to the residual sum of squares (RSS).
"""
deviance(obj::LinearModel) = deviance(obj.rr)

"""
    nulldeviance(obj::LinearModel)

For linear models, the deviance of the null model is equal to the total sum of squares (TSS).
"""
nulldeviance(obj::LinearModel) = nulldeviance(obj.rr)
loglikelihood(obj::LinearModel) = loglikelihood(obj.rr)
nullloglikelihood(obj::LinearModel) = nullloglikelihood(obj.rr)

r2(obj::LinearModel) = 1 - deviance(obj)/nulldeviance(obj)

function adjr2(obj::LinearModel)
    n = nobs(obj)
    # dof() includes the dispersion parameter
    p = dof(obj) - 1
    1 - (1 - r²(obj))*(n-1)/(n-p)
end

function dispersion(x::LinearModel, sqr::Bool=false)
    ssqr = deviance(x.rr)/dof_residual(x)
    return sqr ? ssqr : sqrt(ssqr)
end

function coeftable(mm::LinearModel; level::Real=0.95)
    cc = coef(mm)
    se = stderror(mm)
    tt = cc ./ se
    p = ccdf.(Ref(FDist(1, dof_residual(mm))), abs2.(tt))
    ci = se*quantile(TDist(dof_residual(mm)), (1-level)/2)
    levstr = isinteger(level*100) ? string(Integer(level*100)) : string(level*100)
    CoefTable(hcat(cc,se,tt,p,cc+ci,cc-ci),
              ["Coef.","Std. Error","t","Pr(>|t|)","Lower $levstr%","Upper $levstr%"],
              ["x$i" for i = 1:size(mm.pp.X, 2)], 4, 3)
end

"""
    predict(mm::LinearModel, newx::AbstractMatrix;
            interval::Union{Symbol,Nothing} = nothing, level::Real = 0.95)

If `interval` is `nothing` (the default), return a vector with the predicted values
for model `mm` and new data `newx`.
Otherwise, return a 3-column matrix with the prediction and
the lower and upper confidence bounds for a given `level` (0.95 equates alpha = 0.05).
Valid values of `interval` are `:confidence` delimiting the  uncertainty of the
predicted relationship, and `:prediction` delimiting estimated bounds for new data points.
"""
function predict(mm::LinearModel, newx::AbstractMatrix;
                 interval::Union{Symbol,Nothing}=nothing, level::Real = 0.95)
    retmean = newx * coef(mm)
    if interval === :confint
        Base.depwarn("interval=:confint is deprecated in favor of interval=:confidence")
        interval = :confidence
    end
    if interval === nothing
        return retmean
    end
    length(mm.rr.wts) == 0 || error("prediction with confidence intervals not yet implemented for weighted regression")
    R = cholesky!(mm.pp).U #get the R matrix from the QR factorization
    residvar = (ones(size(newx,2),1) * deviance(mm)/dof_residual(mm))
    if interval == :confidence
        retvariance = (newx/R).^2 * residvar
    elseif interval == :prediction
        retvariance = (newx/R).^2 * residvar .+ deviance(mm)/dof_residual(mm)
    else
        error("only :confidence and :prediction intervals are defined")
    end
    retinterval = quantile(TDist(dof_residual(mm)), (1. - level)/2) * sqrt.(retvariance)
    (prediction = retmean, lower = retmean .+ retinterval, upper = retmean .- retinterval)
end

function confint(obj::LinearModel; level::Real=0.95)
    hcat(coef(obj),coef(obj)) + stderror(obj) *
    quantile(TDist(dof_residual(obj)), (1. - level)/2.) * [1. -1.]
end


text:
# Here we override broadcasting for banded matrices.
# The design is to to exploit the broadcast machinery so that
# banded matrices that conform to the banded matrix interface but are not
# <: AbstractBandedMatrix can get access to fast copyto!, lmul!, rmul!, axpy!, etc.
# using broadcast variants (B .= A, B .= 2.0 .* A, etc.)

abstract type AbstractBlockStyle{N} <: AbstractArrayStyle{N} end
struct BlockStyle{N} <: AbstractBlockStyle{N} end
struct PseudoBlockStyle{N} <: AbstractBlockStyle{N} end

BlockStyle(::Val{N}) where {N} = BlockStyle{N}()
PseudoBlockStyle(::Val{N}) where {N} = PseudoBlockStyle{N}()
BlockStyle{M}(::Val{N}) where {N,M} = BlockStyle{N}()
PseudoBlockStyle{M}(::Val{N}) where {N,M} = PseudoBlockStyle{N}()
BroadcastStyle(::Type{<:BlockArray{<:Any,N}}) where {N} = BlockStyle{N}()
BroadcastStyle(::Type{<:PseudoBlockArray{<:Any,N}}) where {N} = PseudoBlockStyle{N}()
BroadcastStyle(::Type{<:AdjOrTrans{<:Any,<:BlockArray{<:Any,N}}}) where {N} = BlockStyle{2}()
BroadcastStyle(::Type{<:AdjOrTrans{<:Any,<:PseudoBlockArray{<:Any,N}}}) where {N} = PseudoBlockStyle{2}()

BroadcastStyle(::DefaultArrayStyle{N}, b::AbstractBlockStyle{M}) where {M,N} = typeof(b)(Val(max(M,N)))
BroadcastStyle(a::AbstractBlockStyle{N}, ::DefaultArrayStyle{M}) where {M,N} = typeof(a)(Val(max(M,N)))
BroadcastStyle(::StructuredMatrixStyle, b::AbstractBlockStyle{M}) where {M} = typeof(b)(Val(max(M,2)))
BroadcastStyle(a::AbstractBlockStyle{M}, ::StructuredMatrixStyle) where {M} = typeof(a)(Val(max(M,2)))
BroadcastStyle(::BlockStyle{M}, ::PseudoBlockStyle{N}) where {M,N} = BlockStyle(Val(max(M,N)))
BroadcastStyle(::PseudoBlockStyle{M}, ::BlockStyle{N}) where {M,N} = BlockStyle(Val(max(M,N)))

# sortedunion can assume inputs are already sorted so this could be improved
sortedunion(a,b) = sort!(union(a,b))
sortedunion(a::Base.OneTo, b::Base.OneTo) = Base.OneTo(max(last(a),last(b)))
sortedunion(a::AbstractUnitRange, b::AbstractUnitRange) = min(first(a),first(b)):max(last(a),last(b))
combine_blockaxes(a, b) = _BlockedUnitRange(sortedunion(blocklasts(a), blocklasts(b)))

Base.Broadcast.axistype(a::T, b::T) where T<:BlockedUnitRange = length(b) == 1 ? a : combine_blockaxes(a, b)
Base.Broadcast.axistype(a::BlockedUnitRange, b::BlockedUnitRange) = length(b) == 1 ? a : combine_blockaxes(a, b)
Base.Broadcast.axistype(a::BlockedUnitRange, b) = length(b) == 1 ? a : combine_blockaxes(a, b)
Base.Broadcast.axistype(a, b::BlockedUnitRange) = length(b) == 1 ? a : combine_blockaxes(a, b)

similar(bc::Broadcasted{<:AbstractBlockStyle{N}}, ::Type{T}) where {T,N} =
    BlockArray{T,N}(undef, axes(bc))

similar(bc::Broadcasted{PseudoBlockStyle{N}}, ::Type{T}) where {T,N} =
    PseudoBlockArray{T,N}(undef, axes(bc))

"""
    SubBlockIterator(subblock_lasts::Vector{Int}, block_lasts::Vector{Int})
    SubBlockIterator(A::AbstractArray, bs::NTuple{N,AbstractUnitRange{Int}} where N, dim::Integer)

An iterator for iterating `BlockIndexRange` of the blocks specified by
`subblock_lasts`.  The `Block` index part of `BlockIndexRange` is
determined by `subblock_lasts`.  That is to say, the `Block` index first
specifies one of the block represented by `subblock_lasts` and then the
inner-block index range specifies the region within the block.  Each
such block corresponds to a block specified by `blocklasts`.

Note that the invariance `subblock_lasts ⊂ block_lasts` must hold and must
be ensured by the caller.

# Examples
```jldoctest
julia> using BlockArrays

julia> import BlockArrays: SubBlockIterator, BlockIndexRange

julia> A = BlockArray(1:6, 1:3);

julia> subblock_lasts = axes(A, 1).lasts;

julia> @assert subblock_lasts == [1, 3, 6];

julia> block_lasts = [1, 3, 4, 6];

julia> for idx in SubBlockIterator(subblock_lasts, block_lasts)
           B = @show view(A, idx)
           @assert !(parent(B) isa BlockArray)
           idx :: BlockIndexRange
           idx.block :: Block{1}
           idx.indices :: Tuple{UnitRange}
       end
view(A, idx) = 1:1
view(A, idx) = 2:3
view(A, idx) = 4:4
view(A, idx) = 5:6

julia> [idx.block.n[1] for idx in SubBlockIterator(subblock_lasts, block_lasts)]
4-element Vector{Int64}:
 1
 2
 3
 3

julia> [idx.indices[1] for idx in SubBlockIterator(subblock_lasts, block_lasts)]
4-element Vector{UnitRange{Int64}}:
 1:1
 1:2
 1:1
 2:3
```
"""
struct SubBlockIterator
    subblock_lasts::Vector{Int}
    block_lasts::Vector{Int}
end

Base.IteratorEltype(::Type{<:SubBlockIterator}) = Base.HasEltype()
Base.eltype(::Type{<:SubBlockIterator}) = BlockIndexRange{1,Tuple{UnitRange{Int64}}}

Base.IteratorSize(::Type{<:SubBlockIterator}) = Base.HasLength()
Base.length(it::SubBlockIterator) = length(it.block_lasts)

SubBlockIterator(arr::AbstractArray, bs::NTuple{N,AbstractUnitRange{Int}}, dim::Integer) where N =
    SubBlockIterator(blocklasts(axes(arr, dim)), blocklasts(bs[dim]))

function Base.iterate(it::SubBlockIterator, state=nothing)
    if state === nothing
        i,j = 1,1
    else
        i, j = state
    end
    length(it.block_lasts)+1 == i && return nothing
    idx = i == 1 ? (1:it.block_lasts[i]) : (it.block_lasts[i-1]+1:it.block_lasts[i])

    bir = Block(j)[j == 1 ? idx : idx .- it.subblock_lasts[j-1]]
    if it.subblock_lasts[j] == it.block_lasts[i]
        j += 1
    end
    return (bir, (i + 1, j))
end

subblocks(::Any, bs::NTuple{N,AbstractUnitRange{Int}}, dim::Integer) where N =
    (nothing for _ in blockaxes(bs[dim], 1))

function subblocks(arr::AbstractArray, bs::NTuple{N,AbstractUnitRange{Int}}, dim::Integer) where N
    if size(arr, dim) == 1
        return (BlockIndexRange(Block(1), 1:1) for _ in blockaxes(bs[dim], 1))
    end
    return SubBlockIterator(arr, bs, dim)
end

@inline _bview(arg, ::Vararg) = arg
@inline _bview(A::AbstractArray, I...) = view(A, I...)

@inline function Base.Broadcast.materialize!(dest, bc::Broadcasted{BS}) where {Style,BS<:AbstractBlockStyle}
    return copyto!(dest, Base.Broadcast.instantiate(Base.Broadcast.Broadcasted{BS}(bc.f, bc.args, combine_blockaxes.(axes(dest),axes(bc)))))
end

@generated function _generic_blockbroadcast_copyto!(dest::AbstractArray,
                            bc::Broadcasted{<:AbstractBlockStyle{NDims}, <:Any, <:Any, Args}) where {NDims, Args <: Tuple}

    NArgs = length(Args.parameters)

    # `bvar(0, dim)` is a variable for BlockIndexRange of `dim`-th dimension
    # of `dest` array.  `bvar(i, dim)` is a similar variable of `i`-th
    # argument in `bc.args`.
    bvar(i, dim) = Symbol("blockindexrange_", i, "_", dim)

    function forloop(dim)
        if dim > 0
            quote
                for ($(bvar(0, dim)), $(bvar.(1:NArgs, dim)...),) in zip(
                        subblocks(dest, bs, $dim),
                        subblocks.(bc.args, Ref(bs), Ref($dim))...)
                    $(forloop(dim - 1))
                end
            end
        else
            bview(a, i) = :(_bview($a, $([bvar(i, d) for d in 1:NDims]...)))
            destview = bview(:dest, 0)
            argblocks = [bview(:(bc.args[$i]), i) for i in 1:NArgs]
            quote
                broadcast!(bc.f, $destview, $(argblocks...))
            end
        end
    end

    quote
        bs = axes(bc)
        if !blockisequal(axes(dest), bs)
            copyto!(PseudoBlockArray(dest, bs), bc)
            return dest
        end

        $(forloop(NDims))
        return dest
    end
end

copyto!(dest::AbstractArray,
        bc::Broadcasted{<:AbstractBlockStyle{NDims}, <:Any, <:Any, Args}) where {NDims, Args <: Tuple} =
    _generic_blockbroadcast_copyto!(dest, bc)

# type-stable version of _bview.(args, K)
__bview(args::Tuple{}, K) = ()
__bview(args::Tuple, K) = tuple(_bview(args[1],K), __bview(tail(args), K)...)

function _fast_blockbradcast_copyto!(dest, bc)
    @inbounds for K in blockaxes(bc)[1]
        broadcast!(bc.f, view(dest,K), __bview(bc.args, K)...)
    end
    dest
end

_hasscalarlikevec() = false
_hasscalarlikevec(a, b...) = _hasscalarlikevec(b...)
_hasscalarlikevec(a::AbstractVector, b...) = size(a,1) == 1 || _hasscalarlikevec(b...)

blockisequalorscalar(ax, ::Number) = true
blockisequalorscalar(ax, a) = blockisequal(ax, Base.axes1(a))

function copyto!(dest::AbstractVector,
        bc::Broadcasted{<:AbstractBlockStyle{1}, <:Any, <:Any, Args}) where {Args <: Tuple}
    _hasscalarlikevec(bc.args...) && return _generic_blockbroadcast_copyto!(dest, bc)
    ax = axes(dest,1)
    for a in bc.args
        blockisequalorscalar(ax, a) || return _generic_blockbroadcast_copyto!(dest, bc)
    end
    return _fast_blockbradcast_copyto!(dest, bc)
end
@inline function Broadcast.instantiate(bc::Broadcasted{Style}) where {Style <:BlockStyle}
    bcf = Broadcast.instantiate(Broadcast.flatten(Broadcasted{Nothing}(bc.f, bc.args, bc.axes)))
    return Broadcasted{Style}(bcf.f, bcf.args, bcf.axes)
end

_removeblocks(a::Broadcasted) = broadcasted(a.f, map(_removeblocks,a.args)...)
_removeblocks(a::PseudoBlockArray) = a.blocks
_removeblocks(a::BlockSlice) = a.indices
_removeblocks(a::Adjoint) = _removeblocks(parent(a))'
_removeblocks(a::Transpose) = transpose(_removeblocks(parent(a)))
_removeblocks(a::SubArray{<:Any,N,<:PseudoBlockArray}) where N = view(_removeblocks(parent(a)), map(_removeblocks, parentindices(a))...)
_removeblocks(a) = a
copy(bc::Broadcasted{PseudoBlockStyle{N}}) where N = PseudoBlockArray(materialize(_removeblocks(bc)), axes(bc))

for op in (:+, :-, :*)
    @eval function copy(bc::Broadcasted{BlockStyle{N},<:Any,typeof($op),<:Tuple{<:AbstractArray{<:Number,N}}}) where N
        (A,) = bc.args
        _BlockArray(broadcast(a -> broadcast($op, a), blocks(A)), axes(A))
    end
end

for op in (:+, :-, :*, :/, :\)
    @eval begin
        function copy(bc::Broadcasted{BlockStyle{N},<:Any,typeof($op),<:Tuple{<:Number,<:AbstractArray{<:Number,N}}}) where N
            x,A = bc.args
            _BlockArray(broadcast((x,a) -> broadcast($op, x, a), x, blocks(A)), axes(A))
        end
        function copy(bc::Broadcasted{BlockStyle{N},<:Any,typeof($op),<:Tuple{<:AbstractArray{<:Number,N},<:Number}}) where N
            A,x = bc.args
            _BlockArray(broadcast((a,x) -> broadcast($op, a, x), blocks(A),x), axes(A))
        end
    end
end

# exploit special cases for *, for example, *(::Number, ::Diagonal)
for op in (:*, :/)
    @eval @inline $op(A::BlockArray, x::Number) = _BlockArray($op(blocks(A),x), axes(A))
end
for op in (:*, :\)
    @eval @inline $op(x::Number, A::BlockArray) = _BlockArray($op(x,blocks(A)), axes(A))
end

###
# SubViews
###

_blocktype(::Type{<:BlockArray{<:Any,N,<:AbstractArray{R,N}}}) where {N,R} = R

BroadcastStyle(::Type{<:SubArray{T,N,Arr,<:NTuple{N,BlockSlice1},false}}) where {T,N,Arr<:BlockArray} =
    BroadcastStyle(_blocktype(Arr))

# special cases for SubArrays which we want to broadcast by Block
BroadcastStyle(::Type{<:SubArray{<:Any,N,<:Any,I}}) where {N,I<:Tuple{BlockSlice{<:Any,<:BlockedUnitRange},Vararg{Any}}} = BlockStyle{N}()
BroadcastStyle(::Type{<:SubArray{<:Any,N,<:Any,I}}) where {N,I<:Tuple{BlockSlice{<:Any,<:BlockedUnitRange},BlockSlice{<:Any,<:BlockedUnitRange},Vararg{Any}}} = BlockStyle{N}()
BroadcastStyle(::Type{<:SubArray{<:Any,N,<:Any,I}}) where {N,I<:Tuple{Any,BlockSlice{<:Any,<:BlockedUnitRange},Vararg{Any}}} = BlockStyle{N}()

BroadcastStyle(::Type{<:SubArray{<:Any,N,<:PseudoBlockArray,I}}) where {N,I<:Tuple{BlockSlice{<:Any,<:BlockedUnitRange},Vararg{Any}}} = PseudoBlockStyle{N}()
BroadcastStyle(::Type{<:SubArray{<:Any,N,<:PseudoBlockArray,I}}) where {N,I<:Tuple{BlockSlice{<:Any,<:BlockedUnitRange},BlockSlice{<:Any,<:BlockedUnitRange},Vararg{Any}}} = PseudoBlockStyle{N}()
BroadcastStyle(::Type{<:SubArray{<:Any,N,<:PseudoBlockArray,I}}) where {N,I<:Tuple{Any,BlockSlice{<:Any,<:BlockedUnitRange},Vararg{Any}}} = PseudoBlockStyle{N}()

###
# Fill
###

for op in (:*, :/)
    @eval begin
        broadcasted(::AbstractBlockStyle, ::typeof($op), a::Zeros, b::AbstractArray) = FillArrays._broadcasted_zeros(a, b)
        broadcasted(::AbstractBlockStyle, ::typeof($op), a::Ones{T}, b::AbstractArray{V}) where {T,V} = LinearAlgebra.copy_oftype(b, Base.promote_op(*, T, V))
    end
end

for op in (:*, :\)
    @eval begin
        broadcasted(::AbstractBlockStyle, ::typeof($op), a::AbstractArray, b::Zeros) = FillArrays._broadcasted_zeros(a, b)
        broadcasted(::AbstractBlockStyle, ::typeof($op), a::AbstractArray{T}, b::Ones{V}) where {T,V} = LinearAlgebra.copy_oftype(a, Base.promote_op(*, T, V))
    end
end