text:
const P2 = GeometryBasics.Point2{Float64}
const P3 = GeometryBasics.Point3{Float64}

const _haligns = :hcenter, :left, :right
const _valigns = :vcenter, :top, :bottom

nanpush!(a::AVec{P2}, b) = (push!(a, P2(NaN, NaN)); push!(a, b))
nanappend!(a::AVec{P2}, b) = (push!(a, P2(NaN, NaN)); append!(a, b))
nanpush!(a::AVec{P3}, b) = (push!(a, P3(NaN, NaN, NaN)); push!(a, b))
nanappend!(a::AVec{P3}, b) = (push!(a, P3(NaN, NaN, NaN)); append!(a, b))
compute_angle(v::P2) = (angle = atan(v[2], v[1]); angle < 0 ? 2π - angle : angle)

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

struct Shape{X<:Number,Y<:Number}
    x::Vector{X}
    y::Vector{Y}
    # function Shape(x::AVec, y::AVec)
    #     # if x[1] != x[end] || y[1] != y[end]
    #     #     new(vcat(x, x[1]), vcat(y, y[1]))
    #     # else
    #         new(x, y)
    #     end
    # end
end

"""
    Shape(x, y)
    Shape(vertices)

Construct a polygon to be plotted
"""
Shape(verts::AVec) = Shape(RecipesPipeline.unzip(verts)...)
Shape(s::Shape) = deepcopy(s)

get_xs(shape::Shape) = shape.x
get_ys(shape::Shape) = shape.y
vertices(shape::Shape) = collect(zip(shape.x, shape.y))

#deprecated
@deprecate shape_coords coords

"return the vertex points from a Shape or Segments object"
coords(shape::Shape) = shape.x, shape.y

coords(shapes::AVec{<:Shape}) = unzip(map(coords, shapes))

"get an array of tuples of points on a circle with radius `r`"
partialcircle(start_θ, end_θ, n = 20, r = 1) =
    [(r * cos(u), r * sin(u)) for u in range(start_θ, stop = end_θ, length = n)]

"interleave 2 vectors into each other (like a zipper's teeth)"
function weave(x, y; ordering = Vector[x, y])
    ret = eltype(x)[]
    done = false
    while !done
        for o in ordering
            try
                push!(ret, popfirst!(o))
            catch
            end
        end
        done = isempty(x) && isempty(y)
    end
    ret
end

"create a star by weaving together points from an outer and inner circle.  `n` is the number of arms"
function makestar(n; offset = -0.5, radius = 1.0)
    z1 = offset * π
    z2 = z1 + π / (n)
    outercircle = partialcircle(z1, z1 + 2π, n + 1, radius)
    innercircle = partialcircle(z2, z2 + 2π, n + 1, 0.4radius)
    Shape(weave(outercircle, innercircle))
end

"create a shape by picking points around the unit circle.  `n` is the number of point/sides, `offset` is the starting angle"
makeshape(n; offset = -0.5, radius = 1.0) =
    Shape(partialcircle(offset * π, offset * π + 2π, n + 1, radius))

function makecross(; offset = -0.5, radius = 1.0)
    z2 = offset * π
    z1 = z2 - π / 8
    outercircle = partialcircle(z1, z1 + 2π, 9, radius)
    innercircle = partialcircle(z2, z2 + 2π, 5, 0.5radius)
    Shape(
        weave(
            outercircle,
            innercircle,
            ordering = Vector[outercircle, innercircle, outercircle],
        ),
    )
end

from_polar(angle, dist) = P2(dist * cos(angle), dist * sin(angle))

makearrowhead(angle; h = 2.0, w = 0.4, tip = from_polar(angle, h)) = Shape(
    P2[
        (0, 0),
        from_polar(angle - 0.5π, w) - tip,
        from_polar(angle + 0.5π, w) - tip,
        (0, 0),
    ],
)

const _shapes = KW(
    :circle    => makeshape(20),
    :rect      => makeshape(4, offset = -0.25),
    :diamond   => makeshape(4),
    :utriangle => makeshape(3, offset = 0.5),
    :dtriangle => makeshape(3, offset = -0.5),
    :rtriangle => makeshape(3, offset = 0.0),
    :ltriangle => makeshape(3, offset = 1.0),
    :pentagon  => makeshape(5),
    :hexagon   => makeshape(6),
    :heptagon  => makeshape(7),
    :octagon   => makeshape(8),
    :cross     => makecross(offset = -0.25),
    :xcross    => makecross(),
    :vline     => Shape([(0, 1), (0, -1)]),
    :hline     => Shape([(1, 0), (-1, 0)]),
)

for n in 4:8
    _shapes[Symbol("star$n")] = makestar(n)
end

Shape(k::Symbol) = deepcopy(_shapes[k])

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

# uses the centroid calculation from https://en.wikipedia.org/wiki/Centroid#Centroid_of_polygon
"return the centroid of a Shape"
function center(shape::Shape)
    x, y = coords(shape)
    n = length(x)
    A, Cx, Cy = 0, 0, 0
    for i in 1:n
        ip1 = i == n ? 1 : i + 1
        A += x[i] * y[ip1] - x[ip1] * y[i]
    end
    A *= 0.5
    for i in 1:n
        ip1 = i == n ? 1 : i + 1
        m = (x[i] * y[ip1] - x[ip1] * y[i])
        Cx += (x[i] + x[ip1]) * m
        Cy += (y[i] + y[ip1]) * m
    end
    Cx / 6A, Cy / 6A
end

function scale!(shape::Shape, x::Real, y::Real = x, c = center(shape))
    sx, sy = coords(shape)
    cx, cy = c
    for i in eachindex(sx)
        sx[i] = (sx[i] - cx) * x + cx
        sy[i] = (sy[i] - cy) * y + cy
    end
    shape
end

"""
    scale(shape, x, y = x, c = center(shape))
    scale!(shape, x, y = x, c = center(shape))

Scale shape by a factor.
"""
scale(shape::Shape, x::Real, y::Real = x, c = center(shape)) =
    scale!(deepcopy(shape), x, y, c)

function translate!(shape::Shape, x::Real, y::Real = x)
    sx, sy = coords(shape)
    for i in eachindex(sx)
        sx[i] += x
        sy[i] += y
    end
    shape
end

"""
    translate(shape, x, y = x)
    translate!(shape, x, y = x)

Translate a Shape in space.
"""
translate(shape::Shape, x::Real, y::Real = x) = translate!(deepcopy(shape), x, y)

rotate_x(x::Real, y::Real, θ::Real, centerx::Real, centery::Real) =
    ((x - centerx) * cos(θ) - (y - centery) * sin(θ) + centerx)

rotate_y(x::Real, y::Real, θ::Real, centerx::Real, centery::Real) =
    ((y - centery) * cos(θ) + (x - centerx) * sin(θ) + centery)

rotate(x::Real, y::Real, θ::Real, c) = (rotate_x(x, y, θ, c...), rotate_y(x, y, θ, c...))

function rotate!(shape::Shape, θ::Real, c = center(shape))
    x, y = coords(shape)
    for i in eachindex(x)
        xi = rotate_x(x[i], y[i], θ, c...)
        yi = rotate_y(x[i], y[i], θ, c...)
        x[i], y[i] = xi, yi
    end
    shape
end

"rotate an object in space"
function rotate(shape::Shape, θ::Real, c = center(shape))
    x, y = coords(shape)
    x_new = rotate_x.(x, y, θ, c...)
    y_new = rotate_y.(x, y, θ, c...)
    Shape(x_new, y_new)
end

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

mutable struct Font
    family::AbstractString
    pointsize::Int
    halign::Symbol
    valign::Symbol
    rotation::Float64
    color::Colorant
end

"""
    font(args...)
Create a Font from a list of features. Values may be specified either as
arguments (which are distinguished by type/value) or as keyword arguments.
# Arguments
- `family`: AbstractString. "serif" or "sans-serif" or "monospace"
- `pointsize`: Integer. Size of font in points
- `halign`: Symbol. Horizontal alignment (:hcenter, :left, or :right)
- `valign`: Symbol. Vertical aligment (:vcenter, :top, or :bottom)
- `rotation`: Real. Angle of rotation for text in degrees (use a non-integer type)
- `color`: Colorant or Symbol
# Examples
```julia-repl
julia> font(8)
julia> font(family="serif", halign=:center, rotation=45.0)
```
"""
function font(args...; kw...)
    # defaults
    family = "sans-serif"
    pointsize = 14
    halign = :hcenter
    valign = :vcenter
    rotation = 0
    color = colorant"black"

    for arg in args
        T = typeof(arg)
        @assert arg !== :match

        if T == Font
            family = arg.family
            pointsize = arg.pointsize
            halign = arg.halign
            valign = arg.valign
            rotation = arg.rotation
            color = arg.color
        elseif arg == :center
            halign = :hcenter
            valign = :vcenter
        elseif arg ∈ _haligns
            halign = arg
        elseif arg ∈ _valigns
            valign = arg
        elseif T <: Colorant
            color = arg
        elseif T <: Symbol || T <: AbstractString
            try
                color = parse(Colorant, string(arg))
            catch
                family = string(arg)
            end
        elseif T <: Integer
            pointsize = arg
        elseif T <: Real
            rotation = convert(Float64, arg)
        else
            @warn "Unused font arg: $arg ($T)"
        end
    end

    for sym in keys(kw)
        if sym == :family
            family = string(kw[sym])
        elseif sym == :pointsize
            pointsize = kw[sym]
        elseif sym == :halign
            halign = kw[sym]
            halign == :center && (halign = :hcenter)
            @assert halign ∈ _haligns
        elseif sym == :valign
            valign = kw[sym]
            valign == :center && (valign = :vcenter)
            @assert valign ∈ _valigns
        elseif sym == :rotation
            rotation = kw[sym]
        elseif sym == :color
            color = parse(Colorant, kw[sym])
        else
            @warn "Unused font kwarg: $sym"
        end
    end

    Font(family, pointsize, halign, valign, rotation, color)
end

function scalefontsize(k::Symbol, factor::Number)
    f = default(k)
    f = round(Int, factor * f)
    default(k, f)
end

"""
    scalefontsizes(factor::Number)

Scales all **current** font sizes by `factor`. For example `scalefontsizes(1.1)` increases all current font sizes by 10%. To reset to initial sizes, use `scalefontsizes()`
"""
function scalefontsizes(factor::Number)
    for k in keys(merge(_initial_plt_fontsizes, _initial_sp_fontsizes))
        scalefontsize(k, factor)
    end

    for letter in (:x, :y, :z)
        for k in keys(_initial_ax_fontsizes)
            scalefontsize(get_attr_symbol(letter, k), factor)
        end
    end
end

"""
    scalefontsizes()

Resets font sizes to initial default values.
"""
function scalefontsizes()
    for k in keys(merge(_initial_plt_fontsizes, _initial_sp_fontsizes))
        f = default(k)
        if k in keys(_initial_fontsizes)
            factor = f / _initial_fontsizes[k]
            scalefontsize(k, 1.0 / factor)
        end
    end

    for letter in (:x, :y, :z)
        for k in keys(_initial_ax_fontsizes)
            if k in keys(_initial_fontsizes)
                f = default(get_attr_symbol(letter, k))
                factor = f / _initial_fontsizes[k]
                scalefontsize(get_attr_symbol(letter, k), 1.0 / factor)
            end
        end
    end
end

resetfontsizes() = scalefontsizes()

"Wrap a string with font info"
struct PlotText
    str::AbstractString
    font::Font
end
PlotText(str) = PlotText(string(str), font())

"""
    text(string, args...; kw...)

Create a PlotText object wrapping a string with font info, for plot annotations.
`args` and `kw` are passed to `font`.
"""
text(t::PlotText) = t
text(t::PlotText, font::Font) = PlotText(t.str, font)
text(str::AbstractString, f::Font) = PlotText(str, f)
text(str, args...; kw...) = PlotText(string(str), font(args...; kw...))

Base.length(t::PlotText) = length(t.str)

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

struct Stroke
    width
    color
    alpha
    style
end

"""
    stroke(args...; alpha = nothing)

Define the properties of the stroke used in plotting lines
"""
function stroke(args...; alpha = nothing)
    width = 1
    color = :black
    style = :solid

    for arg in args
        T = typeof(arg)

        # if arg in _allStyles
        if allStyles(arg)
            style = arg
        elseif T <: Colorant
            color = arg
        elseif T <: Symbol || T <: AbstractString
            try
                color = parse(Colorant, string(arg))
            catch
            end
        elseif allAlphas(arg)
            alpha = arg
        elseif allReals(arg)
            width = arg
        else
            @warn "Unused stroke arg: $arg ($(typeof(arg)))"
        end
    end

    Stroke(width, color, alpha, style)
end

struct Brush
    size  # fillrange, markersize, or any other sizey attribute
    color
    alpha
end

function brush(args...; alpha = nothing)
    size = 1
    color = :black

    for arg in args
        T = typeof(arg)

        if T <: Colorant
            color = arg
        elseif T <: Symbol || T <: AbstractString
            try
                color = parse(Colorant, string(arg))
            catch
            end
        elseif allAlphas(arg)
            alpha = arg
        elseif allReals(arg)
            size = arg
        else
            @warn "Unused brush arg: $arg ($(typeof(arg)))"
        end
    end

    Brush(size, color, alpha)
end

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

mutable struct SeriesAnnotations
    strs::AVec  # the labels/names
    font::Font
    baseshape::Union{Shape,AVec{Shape},Nothing}
    scalefactor::Tuple
end

_text_label(lab::Tuple, font) = text(lab[1], font, lab[2:end]...)
_text_label(lab::PlotText, font) = lab
_text_label(lab, font) = text(lab, font)

series_annotations(scalar) = series_annotations([scalar])
series_annotations(anns::SeriesAnnotations) = anns
series_annotations(::Nothing) = nothing

function series_annotations(anns::AMat{SeriesAnnotations})
    @assert size(anns, 1) == 1 "matrix of SeriesAnnotations must be a row vector"
    return anns
end

function series_annotations(anns::AMat, outer_args...)
    # Types that represent annotations for an entire series
    whole_series = Union{AVec,Tuple{AVec,Vararg{Any}}}

    # whole_series types can only be in a row vector
    if size(anns, 1) > 1
        for ann in Iterators.filter(ann -> ann isa whole_series, anns)
            (throw ∘ ArgumentError)(
                "Given series annotation must be the only element in its column:\n$ann",
            )
        end
    end

    ann_vec = map(eachcol(anns)) do col
        ann = first(col) isa whole_series ? first(col) : col

        # Override arguments from outer tuple with args from inner tuple
        strs, inner_args = Iterators.peel(wraptuple(ann))
        series_annotations(strs, outer_args..., inner_args...)
    end

    return permutedims(ann_vec)
end

function series_annotations(strs::AVec, args...)
    fnt = font()
    shp = nothing
    scalefactor = 1, 1
    for arg in args
        if isa(arg, Shape) || (isa(arg, AVec) && eltype(arg) == Shape)
            shp = arg
        elseif isa(arg, Font)
            fnt = arg
        elseif isa(arg, Symbol) && haskey(_shapes, arg)
            shp = _shapes[arg]
        elseif isa(arg, Number)
            scalefactor = arg, arg
        elseif is_2tuple(arg)
            scalefactor = arg
        elseif isa(arg, AVec)
            strs = collect(zip(strs, arg))
        else
            @warn "Unused SeriesAnnotations arg: $arg ($(typeof(arg)))"
        end
    end
    # if scalefactor != 1
    #     for s in get(shp)
    #         scale!(s, scalefactor, scalefactor, (0, 0))
    #     end
    # end
    SeriesAnnotations([_text_label(s, fnt) for s in strs], fnt, shp, scalefactor)
end

function series_annotations_shapes!(series::Series, scaletype::Symbol = :pixels)
    anns = series[:series_annotations]
    # msw, msh = anns.scalefactor
    # ms = series[:markersize]
    # msw, msh = if isa(ms, AVec)
    #     1, 1
    # elseif is_2tuple(ms)
    #     ms
    # else
    #     ms, ms
    # end

    # @show msw msh
    if anns !== nothing && anns.baseshape !== nothing
        # we use baseshape to overwrite the markershape attribute
        # with a list of custom shapes for each
        msw, msh = anns.scalefactor
        msize = Float64[]
        shapes = Vector{Shape}(undef, length(anns.strs))
        for i in eachindex(anns.strs)
            str = _cycle(anns.strs, i)

            # get the width and height of the string (in mm)
            sw, sh = text_size(str, anns.font.pointsize)

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

            # we save the size of the larger direction to the markersize list,
            # and then re-scale a copy of baseshape to match the w/h ratio
            maxscale = max(xscale, yscale)
            push!(msize, maxscale)
            baseshape = _cycle(anns.baseshape, i)
            shapes[i] =
                scale(baseshape, msw * xscale / maxscale, msh * yscale / maxscale, (0, 0))
        end
        series[:markershape] = shapes
        series[:markersize] = msize
    end
    return
end

mutable struct EachAnn
    anns
    x
    y
end

function Base.iterate(ea::EachAnn, i = 1)
    if ea.anns === nothing || isempty(ea.anns.strs) || i > length(ea.y)
        return
    end

    tmp = _cycle(ea.anns.strs, i)
    str, fnt = if isa(tmp, PlotText)
        tmp.str, tmp.font
    else
        tmp, ea.anns.font
    end
    ((_cycle(ea.x, i), _cycle(ea.y, i), str, fnt), i + 1)
end

# -----------------------------------------------------------------------
annotations(anns::AMat) = map(annotations, anns)
annotations(sa::SeriesAnnotations) = sa
annotations(anns::AVec) = anns
annotations(anns) = Any[anns]
annotations(::Nothing) = []

_annotationfont(sp::Subplot) = Plots.font(;
    family = sp[:annotationfontfamily],
    pointsize = sp[:annotationfontsize],
    halign = sp[:annotationhalign],
    valign = sp[:annotationvalign],
    rotation = sp[:annotationrotation],
    color = sp[:annotationcolor],
)

_annotation(sp::Subplot, font, lab, pos...; alphabet = "abcdefghijklmnopqrstuvwxyz") = (
    pos...,
    lab == :auto ? text("($(alphabet[sp[:subplot_index]]))", font) : _text_label(lab, font),
)

# Expand arrays of coordinates, positions and labels into individual annotations
# and make sure labels are of type PlotText
function process_annotation(sp::Subplot, xs, ys, labs, font = _annotationfont(sp))
    anns = []
    labs = makevec(labs)
    xlength = length(methods(length, (typeof(xs),))) == 0 ? 1 : length(xs)
    ylength = length(methods(length, (typeof(ys),))) == 0 ? 1 : length(ys)
    for i in 1:max(xlength, ylength, length(labs))
        x, y, lab = _cycle(xs, i), _cycle(ys, i), _cycle(labs, i)
        x = typeof(x) <: TimeType ? Dates.value(x) : x
        y = typeof(y) <: TimeType ? Dates.value(y) : y
        push!(anns, _annotation(sp, font, lab, x, y))
    end
    anns
end

function process_annotation(
    sp::Subplot,
    positions::Union{AVec{Symbol},Symbol,Tuple},
    labs,
    font = _annotationfont(sp),
)
    anns = []
    positions, labs = makevec(positions), makevec(labs)
    for i in 1:max(length(positions), length(labs))
        pos, lab = _cycle(positions, i), _cycle(labs, i)
        push!(anns, _annotation(sp, font, lab, get(_positionAliases, pos, pos)))
    end
    anns
end

_relative_position(xmin, xmax, pos::Length{:pct}) = xmin + pos.value * (xmax - xmin)

# Give each annotation coordinates based on specified position
function locate_annotation(
    sp::Subplot,
    pos::Symbol,
    label::PlotText;
    position_multiplier = Dict{Symbol,Tuple{Float64,Float64}}(
        :topleft      => (0.1pct, 0.9pct),
        :topcenter    => (0.5pct, 0.9pct),
        :topright     => (0.9pct, 0.9pct),
        :bottomleft   => (0.1pct, 0.1pct),
        :bottomcenter => (0.5pct, 0.1pct),
        :bottomright  => (0.9pct, 0.1pct),
    ),
)
    x, y = position_multiplier[pos]
    (
        _relative_position(axis_limits(sp, :x)..., x),
        _relative_position(axis_limits(sp, :y)..., y),
        label,
    )
end
locate_annotation(sp::Subplot, x, y, label::PlotText) = (x, y, label)
locate_annotation(sp::Subplot, x, y, z, label::PlotText) = (x, y, z, label)

locate_annotation(sp::Subplot, rel::NTuple{2,<:Number}, label::PlotText) = (
    _relative_position(axis_limits(sp, :x)..., rel[1] * Plots.pct),
    _relative_position(axis_limits(sp, :y)..., rel[2] * Plots.pct),
    label,
)
locate_annotation(sp::Subplot, rel::NTuple{3,<:Number}, label::PlotText) = (
    _relative_position(axis_limits(sp, :x)..., rel[1] * Plots.pct),
    _relative_position(axis_limits(sp, :y)..., rel[2] * Plots.pct),
    _relative_position(axis_limits(sp, :z)..., rel[3] * Plots.pct),
    label,
)
# -----------------------------------------------------------------------

"type which represents z-values for colors and sizes (and anything else that might come up)"
struct ZValues
    values::Vector{Float64}
    zrange::Tuple{Float64,Float64}
end

function zvalues(
    values::AVec{T},
    zrange::Tuple{T,T} = (ignorenan_minimum(values), ignorenan_maximum(values)),
) where {T<:Real}
    ZValues(collect(float(values)), map(Float64, zrange))
end

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

function expand_extrema!(a::Axis, surf::Surface)
    ex = a[:extrema]
    for vi in surf.surf
        expand_extrema!(ex, vi)
    end
    ex
end

"For the case of representing a surface as a function of x/y... can possibly avoid allocations."
struct SurfaceFunction <: AbstractSurface
    f::Function
end

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

# # I don't want to clash with ValidatedNumerics, but this would be nice:
# ..(a::T, b::T) = (a, b)

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

# style is :open or :closed (for now)
struct Arrow
    style::Symbol
    side::Symbol  # :head (default), :tail, or :both
    headlength::Float64
    headwidth::Float64
end

"""
    arrow(args...)

Define arrowheads to apply to lines - args are `style` (`:open` or `:closed`),
`side` (`:head`, `:tail` or `:both`), `headlength` and `headwidth`
"""
function arrow(args...)
    style = :simple
    side = :head
    headlength = headwidth = 0.3
    setlength = false
    for arg in args
        T = typeof(arg)
        if T == Symbol
            if arg in (:head, :tail, :both)
                side = arg
            else
                style = arg
            end
        elseif T <: Number
            # first we apply to both, but if there's more, then only change width after the first number
            headwidth = Float64(arg)
            if !setlength
                headlength = headwidth
            end
            setlength = true
        elseif T <: Tuple && length(arg) == 2
            headlength, headwidth = Float64(arg[1]), Float64(arg[2])
        else
            @warn "Skipped arrow arg $arg"
        end
    end
    Arrow(style, side, headlength, headwidth)
end

# allow for do-block notation which gets called on every valid start/end pair which
# we need to draw an arrow
function add_arrows(func::Function, x::AVec, y::AVec)
    for i in 2:length(x)
        xyprev = (x[i - 1], y[i - 1])
        xy = (x[i], y[i])
        if ok(xyprev) && ok(xy)
            if i == length(x) || !ok(x[i + 1], y[i + 1])
                # add the arrow from xyprev to xy
                func(xyprev, xy)
            end
        end
    end
end

# -----------------------------------------------------------------------
"create a BezierCurve for plotting"
mutable struct BezierCurve{T<:GeometryBasics.Point}
    control_points::Vector{T}
end

function (bc::BezierCurve)(t::Real)
    p = zero(P2)
    n = length(bc.control_points) - 1
    for i in 0:n
        p += bc.control_points[i + 1] * binomial(n, i) * (1 - t)^(n - i) * t^i
    end
    p
end

@deprecate curve_points coords

coords(curve::BezierCurve, n::Integer = 30; range = [0, 1]) =
    map(curve, Base.range(first(range), stop = last(range), length = n))

function extrema_plus_buffer(v, buffmult = 0.2)
    vmin, vmax = ignorenan_extrema(v)
    vdiff = vmax - vmin
    buffer = vdiff * buffmult
    vmin - buffer, vmax + buffer
end

### Legend

@add_attributes subplot struct Legend
    background_color = :match
    foreground_color = :match
    position = :best
    title = nothing
    font::Font = font(8)
    title_font::Font = font(11)
    column = 1
end :match = (
    :legend_font_family,
    :legend_font_color,
    :legend_title_font_family,
    :legend_title_font_color,
)


text:
@testset "glaser_test.jl" begin

    epl375 = LaserLab.PulsedLaser(375.0nm, 140.0mW, 1.0MHz, 65.0ps)
    lmu_40x_nuv = LaserLab.Objective("LMU-40X-NUV", 5.3mm, 5.1mm, 40.0)
    w0 = lmu_40x_nuv.d /2.0
    gepl375 = LaserLab.GaussianLaser(epl375, w0)
    g2epl375 = LaserLab.propagate_paralell_beam(epl375, lmu_40x_nuv)

    @test gepl375.w0 ≈ w0
    @test gepl375.laser.P ≈ epl375.P
    @test g2epl375.w0 ≈ epl375.λ/(π * lmu_40x_nuv.NA)
    @test LaserLab.spot_size(g2epl375) ≈ 2 * g2epl375.w0
    @test LaserLab.angular_divergence(g2epl375) ≈ 2 * g2epl375.θ0
    @test LaserLab.depth_of_focus(g2epl375)  ≈ 2 * g2epl375.z0

    fI = LaserLab.I(g2epl375)
    I0 = fI(0.0*μm, 0.0*μm )
    @test g2epl375.I0 ≈ I0

    aspot = π * (LaserLab.spot_size(g2epl375) / 2.0)^2
    @test aspot ≈ π * g2epl375.w0^2

    P0 = I0 * aspot
    @test LaserLab.n_photons(g2epl375.laser.λ, P0) ≈ LaserLab.n_photons(g2epl375.laser.λ, g2epl375.I0 * aspot)

    nγ = uconvert(Hz, epl375.P / LaserLab.photon_energy(epl375.λ))
    @test 2 * nγ / (π * gepl375.w0^2) == gepl375.γ0
end

ext: jl

max_stars_repo_name: proscrite/LaserLab

max_issues_repo_name: proscrite/LaserLab

max_forks_repo_name: proscrite/LaserLab


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

using Markdown
using InteractiveUtils

# ╔═╡ bcaefa4e-85c6-11eb-3316-0db4881ebd30
begin
	using PlutoUI
	using NativeSVG
end

# ╔═╡ 13d62b06-85c7-11eb-00b9-0399467b1ab5
md"""# Dictionaries

This chapter presents another built-in type called a dictionary."""

# ╔═╡ 1f1c885e-85c7-11eb-1452-13a8943ebc90
md"""## A Dictionary Is a Mapping

A *dictionary* is like an array, but more general. In an array, the indices have to be integers; in a dictionary they can be (almost) any type.

A dictionary contains a collection of indices, which are called *keys*, and a collection of *values*. Each key is associated with a single value. The association of a key and a value is called a *key-value* pair, or sometimes an *item*.

In mathematical language, a dictionary represents a *mapping* from keys to values, so you can also say that each key “maps to” a value. As an example, we’ll build a dictionary that maps from English to Spanish words, so the keys and the values are all strings.

The function `Dict` creates a new dictionary with no items (because Dict is the name of a built-in function, you should avoid using it as a variable name):

```julia
julia> eng2sp = Dict()
Dict{Any,Any} with 0 entries
```

The types of the keys and values in the dictionary are specified in curly braces: here, both are of type `Any`.

The dictionary is empty. To add items to the dictionary, you can use square brackets:

```julia
julia> eng2sp["one"] = "uno";
```

This line creates an item that maps from the key `"one"` to the value `"uno"`. If we print the dictionary again, we see a key-value pair with an arrow `=>` between the key and value:

```julia
julia> eng2sp
Dict{Any,Any} with 1 entry:
  "one" => "uno"
```

This output format is also an input format. For example, you can create a new dictionary with three items as follows:

```julia
julia> eng2sp = Dict("one" => "uno", "two" => "dos", "three" => "tres") Dict{String,String} with 3 entries:
  "two"   => "dos"
  "one"   => "uno"
  "three" => "tres"
```

Here all the initial keys and values are strings, so a `Dict{String,String}` is created.

The order of the items in a dictionary is unpredictable. If you type the same example
on your computer, you might get a different result.

But that’s not a problem because the elements of a dictionary are never indexed with integer indices. Instead, you use the keys to look up the corresponding values:

```julia
julia> eng2sp["two"]
"dos"
```

The key `"two"` always maps to the value `"dos"`, so the order of the items doesn’t matter.

If the key isn’t in the dictionary, you get an exception:

```julia
julia> eng2sp["four"]
ERROR: KeyError: key "four" not found
```

The `length` function works on dictionaries; it returns the number of key-value pairs:

```julia
julia> length(eng2sp)
3
```

The function keys returns a collection with the keys of the dictionary:

```julia
julia> ks = keys(eng2sp);

julia> print(ks)
["two", "one", "three"]
```

Now you can use the `∈` operator to see whether something appears as a key in the dictionary:

```julia
julia> "one" ∈ ks
true
julia> "uno" ∈ ks
false
```

To see whether something appears as a value in a dictionary, you can use the function `values`, which returns a collection of values, and then use the `∈` operator:

```julia
julia> vs = values(eng2sp);

julia> "uno" ∈ vs
true
```

The `∈` operator uses different algorithms for arrays and dictionaries. For arrays, it searches the elements of the array in order, as described in “Searching”. As the array gets longer, the search time gets longer in direct proportion.

For dictionaries, Julia uses an algorithm called a *hash table* that has a remarkable property: the `∈` operator takes about the same amount of time no matter how many items are in the dictionary.
"""

# ╔═╡ 27e4cf0c-85c8-11eb-0d95-ab9ddc845387
md"""## Dictionaries as Collections of Counters

Suppose you are given a string and you want to count how many times each letter appears. There are several ways you could do it:

* You could create 26 variables, one for each letter of the alphabet. Then you could traverse the string and, for each character, increment the corresponding counter, probably using a chained conditional.
* You could create an array with 26 elements. Then you could convert each character to a number (using the built-in function Int), use the number as an index into the array, and increment the appropriate counter.
* You could create a dictionary with characters as keys and counters as the corresponding values. The first time you see a character, you would add an item to the dictionary. After that you would increment the value of an existing item.

Each of these options performs the same computation, but each of them implements that computation in a different way.

An *implementation* is a way of performing a computation. Some implementations are better than others. For example, an advantage of the dictionary implementation is that we don’t have to know ahead of time which letters appear in the string and we only have to make room for the letters that do appear.

Here is what the code might look like:

```julia
function histogram(s)
	d = Dict()
	for c in s
		if c ∉ keys(d)
			d[c] = 1
		else
			d[c] += 1
		end
	end
	d
end
```

The name of the function is `histogram`, which is a statistical term for a collection of counters (or frequencies).

The first line of the function creates an empty dictionary. The `for` loop traverses the string. Each time through the loop, if the character `c` is not in the dictionary, we create a new item with key `c` and the initial value `1` (since we have seen this letter once). If `c` is already in the dictionary we increment `d[c]`.

Here’s how it works:

```julia
julia> h = histogram("brontosaurus") Dict{Any,Any} with 8 entries:
  'n' => 1
  's' => 2
  'a' => 1
  'r' => 2
  't' => 1
  'o' => 2
  'u' => 2
  'b' => 1
```

The histogram indicates that the letters *a* and *b* appear once, *o* appears twice, and so on.

Dictionaries have a function called `get` that takes a key and a default value. If the key appears in the dictionary, `get` returns the corresponding value; otherwise, it returns the default value. For example:

```julia
julia> h = histogram("a")
Dict{Any,Any} with 1 entry:
  'a' => 1
julia> get(h, 'a', 0)
1
julia> get(h, 'b', 0)
0
```
"""

# ╔═╡ c8dff15e-85c8-11eb-1230-bd39702cf460
md"""#### Exercise 11-1

Use `get` to write `histogram` more concisely. You should be able to eliminate the `if` statement. What does the function `get!` do? Can you use it to simplify `histogram`?
"""

# ╔═╡ 0b2e93b0-85c9-11eb-0dc3-61b1cbfb022b
md"""## Looping and Dictionaries

You can traverse the keys of a dictionary in a for statement. For example, `printhist` prints each key and the corresponding value:

```julia
function printhist(h)
	for c in keys(h)
		println(c, " ", h[c])
	end
end
```

Here’s what the output looks like:

```julia
julia> h = histogram("parrot");

julia> printhist(h)
a 1
r 2
p 1
o 1
t 1
```

Again, the keys are in no particular order. To traverse the keys in sorted order, you can combine `sort` and `collect`:

```julia
julia> for c in sort(collect(keys(h)))
	       println(c, " ", h[c])
       end
a 1
o 1
p 1
r 2
t 1
```
"""

# ╔═╡ 82a2f936-85c9-11eb-1d33-0d9dc98fdac0
md"""## Reverse Lookup

Given a dictionary `d` and a key `k`, it is easy to find the corresponding value `v = d[k]`. This operation is called a *lookup*.

But what if you have `v` and you want to find `k`? You have two problems. First, there might be more than one key that maps to the value `v`. Depending on the application, you might be able to pick one, or you might have to make an array that contains all of them. Second, there is no simple syntax to do a *reverse lookup*; you have to search.

!!! danger
    A reverse lookup is much slower than a forward lookup; if you have to do it often, or if the dictionary gets big, the performance of your program will suffer.

Here is a function that takes a value and returns the first key that maps to that value:

```julia
function reverselookup(d, v)
	for k in keys(d)
		if d[k]==v
			return k
		end
	end
	error("LookupError")
end
```

This function is yet another example of the search pattern, but it uses a function we haven’t seen before: `error`. The `error` function is used to produce an `ErrorException` that interrupts the normal flow of control. In this case it has the message `"LookupError"`, indicating that a key does not exist.

If we get to the end of the loop, that means `v` doesn’t appear in the dictionary as a value, so we throw an exception.

Here is an example of a successful reverse lookup:

```julia
julia> h = histogram("parrot");

julia> key = reverselookup(h, 2)
'r': ASCII/Unicode U+0072 (category Ll: Letter, lowercase)
```

And an unsuccessful one:

```julia
julia> key = reverselookup(h, 3)
ERROR: LookupError
```

The effect when you generate an exception is the same as when Julia throws one: it prints a stacktrace and an error message.

!!! tip
    Julia provides an optimized way to do a reverse lookup: findall(isequal(3), h).
"""

# ╔═╡ 6d87ad02-85ca-11eb-2a3d-999d28a7c4b7
md"""## Dictionaries and Arrays

Arrays can appear as values in a dictionary. For example, if you are given a dictionary that maps from letters to frequencies, you might want to invert it—that is, create a dictionary that maps from frequencies to letters. Since there might be several letters with the same frequency, each value in the inverted dictionary should be an array of letters.

Here is a function that inverts a dictionary:

```julia
function invertdict(d)
	inverse = Dict()
	for key in keys(d)
		val = d[key]
		if val ∉ keys(inverse)
			inverse[val] = [key]
		else
			push!(inverse[val], key)
		end
	end
	inverse
end
```

Each time through the loop, `key` gets a key from `d` and `val` gets the corresponding value. If `val` is not in `inverse`, that means we haven’t seen it before, so we create a new item and initialize it with a *singleton* (an array that contains a single element). Otherwise, we have seen this value before, so we append the corresponding key to the array.

Here is an example:

```julia
julia> hist = histogram("parrot");

julia> inverse = invertdict(hist) Dict{Any,Any} with 2 entries:
  2 => ['r']
  1 => ['a', 'p', 'o', 't']
```

Figure 11-1 is a state diagram showing `hist` and `inverse`. A dictionary is represented as a box with the key-value pairs inside. If the values are integers, floats, or strings, I draw them inside the box, but I usually draw arrays outside the box, just to keep the diagram simple.
"""

# ╔═╡ 4fb90040-85cb-11eb-23dd-ad9ea62b94c3
Drawing(width=720, height=150) do
	defs() do
        marker(id="arrow", markerWidth="10", markerHeight="10", refX="0", refY="3", orient="auto", markerUnits="strokeWidth") do
      		path(d="M0,0 L0,6 L9,3 z", fill="black")
		end
	end
	text(x=105, y=80, font_family="JuliaMono, monospace", text_anchor="end", font_size="0.85rem", font_weight=600) do
		str("hist ->")
	end
	rect(x=120, y=20, width=150, height=110, fill="rgb(242, 242, 242)", stroke="black")
	text(x=140, y=40, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("'o'")
	end
	line(x1=170, y1=35, x2=220, y2=35, stroke="black", marker_end="url(#arrow)")
	text(x=240, y=40, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("1")
	end
	text(x=140, y=60, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("'a'")
	end
	line(x1=170, y1=55, x2=220, y2=55, stroke="black", marker_end="url(#arrow)")
	text(x=240, y=60, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("1")
	end
	text(x=140, y=80, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("'p'")
	end
	text(x=140, y=100, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("'t'")
	end
	text(x=140, y=120, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("'r'")
	end
	line(x1=170, y1=115, x2=220, y2=115, stroke="black", marker_end="url(#arrow)")
	text(x=240, y=120, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("2")
	end
	line(x1=170, y1=95, x2=220, y2=95, stroke="black", marker_end="url(#arrow)")
	text(x=240, y=100, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("1")
	end
	line(x1=170, y1=75, x2=220, y2=75, stroke="black", marker_end="url(#arrow)")
	text(x=240, y=80, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("1")
	end
	text(x=410, y=80, font_family="JuliaMono, monospace", text_anchor="end", font_size="0.85rem", font_weight=600) do
		str("inverse ->")
	end
    rect(x=520, y=10, width=150, height=90, fill="rgb(242, 242, 242)", stroke="black")
	text(x=540, y=30, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("1")
	end
	line(x1=560, y1=25, x2=610, y2=25, stroke="black", marker_end="url(#arrow)")
	text(x=630, y=30, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("'o'")
	end
	text(x=540, y=50, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("2")
	end
	line(x1=560, y1=45, x2=610, y2=45, stroke="black", marker_end="url(#arrow)")
	text(x=630, y=50, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("'a'")
	end
	text(x=540, y=70, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("3")
	end
	line(x1=560, y1=65, x2=610, y2=65, stroke="black", marker_end="url(#arrow)")
	text(x=630, y=70, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("'p'")
	end
	text(x=540, y=90, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("4")
	end
	line(x1=560, y1=85, x2=610, y2=85, stroke="black", marker_end="url(#arrow)")
	text(x=630, y=90, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("'t'")
	end
	rect(x=520, y=110, width=150, height=30, fill="rgb(242, 242, 242)", stroke="black")
	text(x=540, y=130, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("1")
	end
	line(x1=560, y1=125, x2=610, y2=125, stroke="black", marker_end="url(#arrow)")
	text(x=630, y=130, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("'r'")
	end
	rect(x=420, y=10, width=50, height=130, fill="rgb(242, 242, 242)", stroke="black")
	text(x=440, y=60, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("1")
	end
	line(x1=460, y1=55, x2=510, y2=55, stroke="black", marker_end="url(#arrow)")
	text(x=440, y=130, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("2")
	end
	line(x1=460, y1=125, x2=510, y2=125, stroke="black", marker_end="url(#arrow)")
end

# ╔═╡ 26c0d704-85cf-11eb-2a13-293b42210e70
md"*Figure 11-1. State diagram.*"

# ╔═╡ 35e8ee92-85cf-11eb-215c-21912c5fc92b
md"""!!! note
    I mentioned earlier that a dictionary is implemented using a hash table. That means that the keys have to be *hashable*.

    A *hash* is a function that takes a value (of any kind) and returns an integer. Dictionaries use these integers, called hash values, to store and look up key-value pairs.
"""

# ╔═╡ 8316e7a0-85cf-11eb-2ecb-df63279cadc6
md"""## Memos

If you played with the `fibonacci` function from “One More Example”, you might have noticed that the bigger the argument you provide, the longer the function takes to run. Furthermore, the runtime increases quickly.

To understand why, consider Figure 11-2, which shows the *call graph* for `fibonacci` with `n = 4`.
"""

# ╔═╡ c1e519e8-85cf-11eb-2ce1-4bcd389f56bf
Drawing(width=720, height=310) do
	defs() do
        marker(id="arrow", markerWidth="10", markerHeight="10", refX="0", refY="3", orient="auto", markerUnits="strokeWidth") do
      		path(d="M0,0 L0,6 L9,3 z", fill="black")
		end
	end
	rect(x=310, y=10, width=100, height=50, fill="rgb(242, 242, 242)", stroke="black")
	text(x=325, y=30, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("fibonacci")
	end
	text(x=325, y=50, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("n")
	end
	line(x1=340, y1=45, x2=375, y2=45, stroke="black", marker_end="url(#arrow)")
	text(x=390, y=50, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("4")
	end
	line(x1=310, y1=60, x2=260, y2=85, stroke="black", marker_end="url(#arrow)")
	line(x1=410, y1=60, x2=460, y2=85, stroke="black", marker_end="url(#arrow)")
	rect(x=150, y=90, width=100, height=50, fill="rgb(242, 242, 242)", stroke="black")
	text(x=165, y=110, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("fibonacci")
	end
	text(x=165, y=130, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("n")
	end
	line(x1=180, y1=125, x2=215, y2=125, stroke="black", marker_end="url(#arrow)")
	text(x=230, y=130, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("3")
	end
	line(x1=150, y1=140, x2=130, y2=162, stroke="black", marker_end="url(#arrow)")
	rect(x=75, y=170, width=100, height=50, fill="rgb(242, 242, 242)", stroke="black")
	text(x=90, y=190, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("fibonacci")
	end
	text(x=90, y=210, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("n")
	end
	line(x1=105, y1=205, x2=140, y2=205, stroke="black", marker_end="url(#arrow)")
	text(x=155, y=210, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("2")
	end
	line(x1=250, y1=140, x2=270, y2=162, stroke="black", marker_end="url(#arrow)")
	rect(x=225, y=170, width=100, height=50, fill="rgb(242, 242, 242)", stroke="black")
	text(x=240, y=190, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("fibonacci")
	end
	text(x=240, y=210, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("n")
	end
	line(x1=255, y1=205, x2=290, y2=205, stroke="black", marker_end="url(#arrow)")
	text(x=305, y=210, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("1")
	end
	line(x1=75, y1=220, x2=53, y2=242, stroke="black", marker_end="url(#arrow)")
	rect(x=0, y=250, width=100, height=50, fill="rgb(242, 242, 242)", stroke="black")
	text(x=15, y=270, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("fibonacci")
	end
	text(x=15, y=290, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("n")
	end
	line(x1=30, y1=285, x2=65, y2=285, stroke="black", marker_end="url(#arrow)")
	text(x=80, y=290, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("1")
	end
	line(x1=175, y1=220, x2=197, y2=242, stroke="black", marker_end="url(#arrow)")
	rect(x=150, y=250, width=100, height=50, fill="rgb(242, 242, 242)", stroke="black")
	text(x=150+15, y=270, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("fibonacci")
	end
	text(x=150+15, y=290, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("n")
	end
	line(x1=150+30, y1=285, x2=150+65, y2=285, stroke="black", marker_end="url(#arrow)")
	text(x=150+80, y=290, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("0")
	end
	rect(x=470, y=90, width=100, height=50, fill="rgb(242, 242, 242)", stroke="black")
	text(x=485, y=110, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("fibonacci")
	end
	text(x=485, y=130, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("n")
	end
	line(x1=500, y1=125, x2=535, y2=125, stroke="black", marker_end="url(#arrow)")
	text(x=550, y=130, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("2")
	end
	line(x1=570, y1=140, x2=590, y2=162, stroke="black", marker_end="url(#arrow)")
	rect(x=720-175, y=170, width=100, height=50, fill="rgb(242, 242, 242)", stroke="black")
	text(x=720-175+15, y=190, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("fibonacci")
	end
	text(x=720-175+15, y=210, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("n")
	end
	line(x1=720-175+30, y1=205, x2=720-175+65, y2=205, stroke="black", marker_end="url(#arrow)")
	text(x=720-175+80, y=210, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("0")
	end
	line(x1=470, y1=140, x2=450, y2=162, stroke="black", marker_end="url(#arrow)")
	rect(x=720-325, y=170, width=100, height=50, fill="rgb(242, 242, 242)", stroke="black")
	text(x=720-325+15, y=190, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("fibonacci")
	end
	text(x=720-325+15, y=210, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("n")
	end
	line(x1=720-325+30, y1=205, x2=720-325+65, y2=205, stroke="black", marker_end="url(#arrow)")
	text(x=720-325+80, y=210, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("1")
	end
end

# ╔═╡ 2a7646fe-85d9-11eb-3bbe-95f975f3901f
md"*Figure 11-2. Call graph.*"

# ╔═╡ 388a8d90-85d9-11eb-08f1-6b622cfb2c11
md"""A call graph shows a set of function frames, with lines connecting each frame to the frames of the functions it calls. At the top of the graph, `fibonacci` with `n = 4` calls `fibonacci` with `n = 3` and `n = 2`. In turn, `fibonacci` with `n = 3` calls `fibonacci` with `n = 2`and `n = 1`, and so on.

Count how many times `fibonacci(0)` and `fibonacci(1)` are called. This is an inefficient solution to the problem, and it gets worse as the argument gets bigger.

One solution is to keep track of values that have already been computed by storing them in a dictionary. A previously computed value that is stored for later use is called a *memo*. Here is a “memoized” version of `fibonacci`:

```julia
const known = Dict(0=>0, 1=>1)

function fibonacci(n)
	if n ∈ keys(known)
		return known[n]
	end
    res = fibonacci(n-1) + fibonacci(n-2)
	known[n] = res
	res
end
```

`known` is a constant dictionary that keeps track of the Fibonacci numbers we already know. It starts with two items: `0` maps to `0` and `1` maps to `1`.

Whenever `fibonacci` is called, it checks `known`. If the result is already there, it can return immediately. Otherwise, it has to compute the new value, add it to the dictionary, and return it.

If you run this version of `fibonacci` and compare it with the original, you will find that it is much faster.
"""

# ╔═╡ 170c40ea-85da-11eb-2079-a912e41da021
md"""## Global Variables

In the previous example, `known` is created outside the function, so it belongs to the special frame called `Main`. Variables in `Main` are sometimes called *global* because they can be accessed from any function. Unlike local variables, which disappear when their function ends, global variables persist from one function call to the next.

For performance reasons, you should declare a global variable *constant* with the keyword `const`. You can no longer reassign the variable, but if it refers to a mutable value, you can modify the value.

!!! danger
    Global variables can be useful, but if you have a lot of them, and you modify them frequently, they can make programs hard to debug and perform badly.
"""

# ╔═╡ 698eb48a-85da-11eb-272f-157bb3fec731
md"""## Debugging

As you work with bigger datasets, it can become unwieldy to debug by printing and checking the output by hand. Here are some suggestions for debugging large datasets:
* Scale down the input.

  If possible, reduce the size of the dataset. For example, if the program reads a text file, start with just the first 10 lines, or with the smallest example you can find that errors. You should not edit the files themselves, but rather modify the pro‐ gram so it reads only the first n lines.

  If there is an error, you can reduce n to the smallest value that manifests the error, and then increase it gradually as you find and correct errors.

* Check summaries and types.

  Instead of printing and checking the entire dataset, consider printing summaries of the data: for example, the number of items in a dictionary or the total of an array of numbers.

  A common cause of runtime errors is a value that is not the right type. For debugging this kind of error, it is often enough to print the type of a value.

* Write self-checks.

  Sometimes you can write code to check for errors automatically. For example, if you are computing the average of an array of numbers, you could check that the result is not greater than the largest element in the array or less than the smallest. This is called a “sanity check.”

  Another kind of check compares the results of two different computations to see if they are consistent. This is called a “consistency check.”

* Format the output.

  Formatting debugging output can make it easier to spot an error. We saw an example in “Debugging”.

  Again, time you spend building scaffolding can reduce the time you spend debugging.
"""

# ╔═╡ 057b7444-85db-11eb-1dc3-179e19732c1c
md"""## Glossary

*dictionary*:
A mapping from keys to their corresponding values.

*key*:
An object that appears in a dictionary as the first part of a key-value pair.

*value*:
An object that appears in a dictionary as the second part of a key-value pair. This is more specific than our previous use of the word “value.”

*key-value pair*:
The representation of the mapping from a key to a value.

*item*:
In a dictionary, another name for a key-value pair.

*mapping*:
A relationship in which each element of one set corresponds to an element of another set.

*hash table*:
The algorithm used to implement Julia dictionaries.

*implementation*:
A way of performing a computation.

*lookup*:
A dictionary operation that takes a key and finds the corresponding value.

*reverse lookup*:
A dictionary operation that takes a value and finds one or more keys that map to it.

*singleton*:
An array (or other sequence) with a single element.

*hashable*:
A type that has a hash function.

*hash function*:
A function used by a hash table to compute the location for a key.

*call graph*:
A diagram that shows every frame created during the execution of a program, with an arrow from each caller to each callee.

*memo*:
A computed value stored to avoid unnecessary future computation.

*global variable*:
A variable defined outside a function. Global variables can be accessed from any function.

*constant global variable*
A global variable that cannot be reassigned.
"""

# ╔═╡ 64397616-85db-11eb-0150-fd07faf0ee9b
md"""## Exercises

#### Exercise 11-2

Write a function that reads the words in words.txt and stores them as keys in a dictionary. It doesn’t matter what the values are. Then you can use the `∈` operator as a fast way to check whether a string is in the dictionar.

If you did “Exercise 10-10”, you can compare the speed of this implementation with the array `∈` operator and the bisection search.
"""

# ╔═╡ 9513aed2-85db-11eb-19ee-dd547267503b
md"""#### Exercise 11-3

Read the documentation of the dictionary function `get!` and use it to write a more concise version of `invertdict`.
"""

# ╔═╡ a5cd8176-85db-11eb-193c-75864c42871e
md"""#### Exercise 11-4

Memoize the Ackermann function from “Exercise 6-5” and see if memoization makes it possible to evaluate the function with bigger arguments.
"""

# ╔═╡ b8989c78-85db-11eb-3a4a-c3bb2dc3a988
md"""#### Exercise 11-5

If you did “Exercise 10-7”, you already have a function named `hasduplicates` that takes an array as a parameter and returns true if there is any object that appears more than once in the array.

Use a dictionary to write a faster, simpler version of `hasduplicates`.
"""

# ╔═╡ d4c39484-85db-11eb-0ed1-a3268690454e
md"""#### Exercise 11-6

Two words are “rotate pairs” if you can rotate one of them and get the other (see rotateword in “Exercise 8-11”).

Write a program that reads a word array and finds all the rotate pairs.
"""

# ╔═╡ e388e7f8-85db-11eb-0451-df6af9de2adc
md"""#### Exercise 11-7

Here’s another Puzzler from Car Talk:

> [A contributor] came upon a common one-syllable, five-letter word recently that has the following unique property. When you remove the first letter, the remaining letters form a homophone of the original word, that is a word that sounds exactly the same. Replace the first letter, that is, put it back and remove the second letter and the result is yet another homophone of the original word. And the question is, what’s the word?
>
> Now I’m going to give you an example that doesn’t work. Let’s look at the five-letter word, ‘wrack.’ W-R-A-C-K, you know like to ‘wrack with pain.’ If I remove the first letter, I am left with a four-letter word, ‘R-A-C-K.’ As in, ‘Holy cow, did you see the rack on that buck! It must have been a nine-pointer!’ It’s a perfect homophone. If you put the ‘w’ back, and remove the ‘r,’ instead, you’re left with the word, ‘wack,’ which is a real word, it’s just not a homophone of the other two words.
>
> But there is, however, at least one word that [I] know of, which will yield two homophones if you remove either of the first two letters to make two, new four-letter words. The question is, what’s the word?

You can use the dictionary from “Exercise 11-2” to check whether a string is in the word array.
!!! tip
    To check whether two words are homophones, you can use the Carnegie Mellon University Pronouncing Dictionary.

Write a program that lists all the words that solve the Puzzler.
"""

# ╔═╡ Cell order:
# ╟─bcaefa4e-85c6-11eb-3316-0db4881ebd30
# ╟─13d62b06-85c7-11eb-00b9-0399467b1ab5
# ╟─1f1c885e-85c7-11eb-1452-13a8943ebc90
# ╟─27e4cf0c-85c8-11eb-0d95-ab9ddc845387
# ╟─c8dff15e-85c8-11eb-1230-bd39702cf460
# ╟─0b2e93b0-85c9-11eb-0dc3-61b1cbfb022b
# ╟─82a2f936-85c9-11eb-1d33-0d9dc98fdac0
# ╟─6d87ad02-85ca-11eb-2a3d-999d28a7c4b7
# ╟─4fb90040-85cb-11eb-23dd-ad9ea62b94c3
# ╟─26c0d704-85cf-11eb-2a13-293b42210e70
# ╟─35e8ee92-85cf-11eb-215c-21912c5fc92b
# ╟─8316e7a0-85cf-11eb-2ecb-df63279cadc6
# ╟─c1e519e8-85cf-11eb-2ce1-4bcd389f56bf
# ╟─2a7646fe-85d9-11eb-3bbe-95f975f3901f
# ╟─388a8d90-85d9-11eb-08f1-6b622cfb2c11
# ╟─170c40ea-85da-11eb-2079-a912e41da021
# ╟─698eb48a-85da-11eb-272f-157bb3fec731
# ╟─057b7444-85db-11eb-1dc3-179e19732c1c
# ╟─64397616-85db-11eb-0150-fd07faf0ee9b
# ╟─9513aed2-85db-11eb-19ee-dd547267503b
# ╟─a5cd8176-85db-11eb-193c-75864c42871e
# ╟─b8989c78-85db-11eb-3a4a-c3bb2dc3a988
# ╟─d4c39484-85db-11eb-0ed1-a3268690454e
# ╟─e388e7f8-85db-11eb-0451-df6af9de2adc


text:
# GUI event loops and toolkit integration for Python, most importantly
# to support plotting in a window without blocking.  Currently, we
# support wxWidgets, Qt4, and GTK+.

############################################################################

# global variable to specify default GUI toolkit to use
gui = :default # one of :default, :wx, :qt, :tk, or :gtk

pyexists(mod) = try
    pyimport(mod)
    true
  catch
    false
  end

# Tkinter was renamed to tkinter in Python 3
function tkinter_name()
    return pyversion.major < 3 ? "Tkinter" : "tkinter"
end

pygui_works(gui::Symbol) = gui == :default ||
    ((gui == :wx && pyexists("wx")) ||
     (gui == :gtk && pyexists("gtk")) ||
     (gui == :gtk3 && pyexists("gi")) ||
     (gui == :tk && pyexists(tkinter_name())) ||
     (gui == :qt_pyqt4 && pyexists("PyQt4")) ||
     (gui == :qt_pyside && pyexists("PySide")) ||
     (gui == :qt && (pyexists("PyQt4") || pyexists("PySide"))))

# get or set the default GUI; doesn't affect running GUI
"""
    pygui()

Return the current GUI toolkit as a symbol.
"""
function pygui()
    global gui
    if gui::Symbol != :default && pygui_works(gui::Symbol)
        return gui::Symbol
    else
        for g in (:tk, :qt, :wx, :gtk, :gtk3)
            if pygui_works(g)
                gui = g
                return gui::Symbol
            end
        end
        error("No supported Python GUI toolkit is installed.")
    end
end
function pygui(g::Symbol)
    global gui
    if g != gui::Symbol
        if !(g in (:wx,:gtk,:gtk3,:tk,:qt,:qt_pyqt4,:qt_pyside,:default))
            throw(ArgumentError("invalid gui $g"))
        elseif !pygui_works(g)
            error("Python GUI toolkit for $g is not installed.")
        end
        gui = g
    end
    return g::Symbol
end

############################################################################
# Event loops for various toolkits.

# call doevent(status) every sec seconds
function install_doevent(doevent::Function, sec::Real)
    timeout = Base.Timer(doevent,sec,sec)
    return timeout
end

macro doevent(body)
    Expr(:function, :($(esc(:doevent))(async)), body)
end

# For PyPlot issue #181: recent pygobject releases emit a warning
# if we don't specify which version we want:
function gtk_requireversion(gtkmodule::AbstractString, vers::VersionNumber=v"3.0")
    if startswith(gtkmodule, "gi.")
        gi = pyimport("gi")
        if gi[:get_required_version]("Gtk") === nothing
            gi[:require_version]("Gtk", string(vers))
        end
    end
end

# GTK (either GTK2 or GTK3, depending on gtkmodule):
function gtk_eventloop(gtkmodule::AbstractString, sec::Real=50e-3)
    gtk_requireversion(gtkmodule)
    gtk = pyimport(gtkmodule)
    events_pending = gtk["events_pending"]
    main_iteration = gtk["main_iteration"]
    @doevent begin
        # handle all pending
        while pycall(events_pending, Bool)
            pycall(main_iteration, PyObject)
        end
    end
    install_doevent(doevent, sec)
end

# Qt4: (PyQt4 or PySide module)
function qt_eventloop(QtModule="PyQt4", sec::Real=50e-3)
    QtCore = pyimport("$QtModule.QtCore")
    instance = QtCore["QCoreApplication"]["instance"]
    AllEvents = QtCore["QEventLoop"]["AllEvents"]
    processEvents = QtCore["QCoreApplication"]["processEvents"]
    maxtime = PyObject(50)
    @doevent begin
        app = pycall(instance, PyObject)
        if app.o != pynothing[]
            app["_in_event_loop"] = true
            pycall(processEvents, PyObject, AllEvents, maxtime)
        end
    end
    install_doevent(doevent, sec)
end

# wx:  (based on IPython/lib/inputhookwx.py, which is 3-clause BSD-licensed)
function wx_eventloop(sec::Real=50e-3)
    wx = pyimport("wx")
    GetApp = wx["GetApp"]
    EventLoop = wx["EventLoop"]
    EventLoopActivator = wx["EventLoopActivator"]
    @doevent begin
        app = pycall(GetApp, PyObject)
        if app.o != pynothing[]
            app["_in_event_loop"] = true
            evtloop = pycall(EventLoop, PyObject)
            ea = pycall(EventLoopActivator, PyObject, evtloop)
            Pending = evtloop["Pending"]
            Dispatch = evtloop["Dispatch"]
            while pycall(Pending, Bool)
                pycall(Dispatch, PyObject)
            end
            pydecref(ea) # deactivate event loop
            pycall(app["ProcessIdle"], PyObject)
        end
    end
    install_doevent(doevent, sec)
end

# Tk: (Tkinter/tkinter module)
function Tk_eventloop(sec::Real=50e-3)
    Tk = pyimport(tkinter_name())
    @doevent begin
        root = Tk["_default_root"]
        if root.o != pynothing[]
            pycall(root["update"], PyObject)
        end
    end
    install_doevent(doevent, sec)
end
# cache running event loops (so that we don't start any more than once)
const eventloops = Dict{Symbol,Timer}()

"""
    pygui_start(gui::Symbol = pygui())

Start the event loop of a certain toolkit.

The argument `gui` defaults to the current default GUI, but it could be `:wx`, `:gtk`, `:gtk3`, `:tk`, or `:qt`.

"""
function pygui_start(gui::Symbol=pygui(), sec::Real=50e-3)
    pygui(gui)
    if !haskey(eventloops, gui)
        if gui == :wx
            eventloops[gui] = wx_eventloop(sec)
        elseif gui == :gtk
            eventloops[gui] = gtk_eventloop("gtk", sec)
        elseif gui == :gtk3
            eventloops[gui] = gtk_eventloop("gi.repository.Gtk", sec)
        elseif gui == :tk
            eventloops[gui] = Tk_eventloop(sec)
        elseif gui == :qt_pyqt4
            eventloops[gui] = qt_eventloop("PyQt4", sec)
        elseif gui == :qt_pyside
            eventloops[gui] = qt_eventloop("PySide", sec)
        elseif gui == :qt
            try
                eventloops[gui] = qt_eventloop("PyQt4", sec)
            catch
                eventloops[gui] = qt_eventloop("PySide", sec)
            end
        else
            throw(ArgumentError("unsupported GUI type $gui"))
        end
    end
    gui
end

"""
    pygui_stop(gui::Symbol = pygui())

Stop any running event loop for gui. The `gui` argument defaults to current default GUI.

"""
function pygui_stop(gui::Symbol=pygui())
    if haskey(eventloops, gui)
        Base.close(pop!(eventloops, gui))
        true
    else
        false
    end
end

pygui_stop_all() = for gui in keys(eventloops); pygui_stop(gui); end

############################################################################


text:
using PowerModels
using Ipopt
using Gurobi
using JuMP
using .PP2PM

import JSON

function read_time_series(json_path)
    time_series = Dict()
    open(json_path, "r") do f
        time_series = JSON.parse(f)  # parse and transform data
    end
    return time_series
end

function set_pq_values_from_timeseries(mn, time_series)
    # This function iterates over multinetwork entries and sets p, q values
    # of loads and "sgens" (which are loads with negative P and Q values)

    # iterate over networks (which represent the time steps)
    for (t, network) in mn["nw"]
        t_j = string(parse(Int64,t) - 1)
        # iterate over all loads for this network
        for (i, load) in network["load"]
            # update variables from time series here
#             print("\nload before: ")
#             print(load["pd"])
            load["pd"] = time_series[t_j][parse(Int64,i)] / mn["baseMVA"]
#             print("\nload after: ")
#             print(load["pd"])
        end
    end

    return mn
end

function run_powermodels(json_path)
    # load converted pandapower network
    pm = PP2PM.load_pm_from_json(json_path)
    # copy network n_time_steps time step times
    n_time_steps = pm["n_time_steps"]
    mn = PowerModels.replicate(pm, pm["n_time_steps"])
    mn["time_elapsed"] = pm["time_elapsed"]
    # set P, Q values of loads and generators from time series
    if isfile("/tmp/timeseries.json")
        time_series = read_time_series("/tmp/timeseries.json")
        mn = set_pq_values_from_timeseries(mn, time_series)
    else
        print("Running storage without time series")
    end

    ipopt_solver = JuMP.with_optimizer(Ipopt.Optimizer, print_level = 0)

    # run multinetwork storage opf
    result = PowerModels._run_mn_strg_opf(mn, PowerModels.ACPPowerModel, ipopt_solver)
    print_summary(result)
    return result
end

max_stars_repo_name: julffers/pandapower

max_issues_repo_name: shravanchipli/pandapower

max_forks_repo_name: gdgarcia/pandapower


text:
using Test, Match

input = joinpath(@__DIR__, "input")
program = readlines(input)

struct Op
    code::String
    value::Int
end

function parse_program(p)
    return map(parse_op, p)
end

function parse_op(op)
    code = op[1:3]
    val = parse(Int, op[5:end])
    return Op(code, val)
end

function run(program, ptr = 1, visited = Set{Int}())
    acc = 0
    v = copy(visited)
    bound = length(program)
    while !(ptr in v) && ptr > 0 && ptr <= bound
        op = program[ptr]
        code = op.code
        val = op.value
        push!(v, ptr)
        @match code begin
            "acc" => begin acc += val; ptr += 1 end
            "jmp" => begin ptr += val end
            "nop" => begin ptr += 1 end
        end
    end
    complete = ptr == length(program) + 1
    return [complete, acc, v]
end

function mutate(program, visited)
    winning = Set{Int}()
    losing = Set{Int}()
    for i in 1:length(program)
        if !(i in union(winning, losing))
            res = run(program, i)
            res[1] ? union!(winning, res[3]) : union!(losing, res[3])
        end
    end
    for i in visited
        op = program[i]
        if op.code == "jmp"
            if i + 1 in winning
                new = Op("nop", op.value)
                dup = copy(program)
                dup[i] = new
                return dup
            end
        elseif op.code == "nop"
            if i + op.value in winning
                new = Op("jmp", op.value)
                dup = copy(program)
                dup[i] = new
                return dup
            end
        end
    end
end

prog = parse_program(program)
r1 = run(prog)

p1 = r1[2]

mutant = mutate(prog, r1[3])
r2 = run(mutant)
p2 = r2[2]

@assert(p1 == 1489)
@assert(p2 == 1539)

sample = [
    "nop +0",
    "acc +1",
    "jmp +4",
    "acc +3",
    "jmp -3",
    "acc -99",
    "acc +1",
    "jmp -4",
    "acc +6"
]

t = parse_program(sample)
r = run(t)
m = mutate(t, r[3])
v = run(m)

@test r[2] == 5
@test r[1] == false

@test v[2] == 8
@test v[1] == true

println("-----------------------------------------------------------------------")
println("handheld halting -- part one :: $p1")
println("handheld halting -- part two :: $p2")
println("-----------------------------------------------------------------------")


text:
u1 = (L = 10.0,)
p1 = (k_L = 1.0, k_R = 1/20, m = 2.0)
# make sure that all states go before parameters
popt = SLVector(L = 10.1, k_L = 1.1, k_R = 1/20.1)

# use allow_missing_opt = false for type stability
ps = @inferred ProblemParSetter_sym(keys(u1),keys(p1),keys(popt), Val(false))

@testset "warning on missing symbols" begin
    state_syms = keys(u1)
    par_syms = keys(p1)
    popt_syms = (:L, :k_L, :M1, :M2)
    NO = length(popt_syms)
    psw = @test_logs (:warn,r"missing optimization parameters") ProblemParSetter_sym(state_syms, par_syms, popt_syms)
end;

@testset "MethodError if missing symbols are not allowed" begin
    @test_throws MethodError ps1 = ProblemParSetter_sym((:x, :RHS), (:τ,), (:RHS, :τ, :M), Val(false))
end;

@testset "access symbols and counts" begin
    @test (@inferred symbols_state(ps)) == keys(u1)
    @test (@inferred symbols_par(ps)) == keys(p1)
    @test (@inferred symbols_paropt(ps)) == keys(popt)
    #
    @test (@inferred count_state(ps)) == length(u1)
    @test (@inferred count_par(ps)) == length(p1)
    @test (@inferred count_paropt(ps)) == length(popt)
end;

# ps.statemap
# ps.optinfo

@testset "label Tuples and SVectors" begin
    # not type stable, because names are not parameters of ps
    @test label_paropt(ps, Tuple(popt)) == NamedTuple(popt)
    @test label_paropt(ps, SVector(popt)) == popt
    @test label_paropt(ps, collect(popt)) == LArray(popt)
    #
    @test label_state(ps, Tuple(u1)) == u1
    @test label_par(ps, Tuple(p1)) == p1
    @test label_state(ps, SVector(Tuple(u1))) == SLVector(u1)
    @test label_par(ps, SVector(Tuple(p1))) == SLVector(p1)
    @test label_state(ps, collect(u1)) == LVector(u1)
    @test label_par(ps, collect(p1)) == LVector(p1)
end;

@testset "update_statepar and get_paropt" begin
    u0o, po = @inferred update_statepar(ps, popt, u1, p1)
    #@btime update_statepar($ps, $popt, $u1, $p1) # zero allocations
    @test collect(u0o) == [10.1]
    @test collect(po) == [1.1,1/20.1,2.0]
    #
    # retrieve popt again
    # using Cthulhu
    # @descend_code_warntype get_paropt(ps, u0o, po)
    popt2 = @inferred get_paropt(ps, u0o, po)
    @test all(popt2 .== popt)
    #using BenchmarkTools
    #@btime get_paropt($ps, $u0o, $po)
    #
    # no @inferred because labelling causes type instability
    popt2n = get_paropt_labeled(ps, u0o, po)
    @test popt2n == NamedTuple(popt)
    #
    popt2m = get_paropt_labeled(ps, collect(u0o), collect(po))
    @test popt2m == LArray(popt)
end;

@testset "update_statepar and get_paropt for Vector" begin
    u1vec = collect(u1)
    p1vec = collect(p1)
    u0o, po = @inferred update_statepar(ps, popt, u1vec, p1vec)
    #@btime update_statepar($ps, $popt, $u1, $p1) # zero allocations
    @test collect(u0o) == [10.1]
    @test collect(po) == [1.1,1/20.1,2.0]
    #
    # retrieve popt again
    #@code_warntype get_paropt(ps, u0o, po)
    popt2 = @inferred get_paropt(ps, u0o, po)
    @test all(popt2 .== popt)
    #using BenchmarkTools
    #@btime get_paropt($ps, $u0o, $po)
    #
    # no @inferred because labelling ist not type stable
    popt2n = get_paropt_labeled(ps, u0o, po)
    @test popt2n == LArray(popt)
end;

# vcat AxisArray does not recturn AxisArray
# @testset "upate_statepar and get_paropt for AxisArray" begin
#     # test with AbstractVector different from Vector
#     u1vec = AxisArray(collect(u1); row = keys(u1))
#     p1vec = collect(p1)
#     u0o, po = @inferred update_statepar(ps, popt, u1vec, p1vec)
#     #@btime update_statepar($ps, $popt, $u1, $p1) # zero allocations
#     @test collect(u0o) == [10.1]
#     @test collect(po) == [1.1,1/20.1,2.0]
#     #
#     # retrieve popt again
#     #@code_warntype get_paropt(ps, u0o, po)
#     popt2 = @inferred get_paropt(ps, u0o, po)
#     @test all(popt2 .== popt)
#     #using BenchmarkTools
#     #@btime get_paropt($ps, $u0o, $po)
#     #
#     # no @inferred because labelling ist not type stable
#     popt2n = get_paropt_labeled(ps, u0o, po)
#     @test popt2n == LArray(popt)
# end;

@testset "update_statepar and get_paropt for p only AbstractVector" begin
    p1vec = collect(p1)
    u0o, po = @inferred update_statepar(ps, popt, u1, p1vec)
    #@btime update_statepar($ps, $popt, $u1, $p1) # zero allocations
    @test collect(u0o) == [10.1]
    @test collect(po) == [1.1,1/20.1,2.0]
    #
    # retrieve popt again
    #@code_warntype get_paropt(ps, u0o, po)
    popt2 = @inferred get_paropt(ps, u0o, po)
    @test all(popt2 .== popt)
    #using BenchmarkTools
    #@btime get_paropt($ps, $u0o, $po)
    #
    # no @inferred because labelling ist not type stable
    popt2n = get_paropt_labeled(ps, u0o, po)
    @test popt2n == NamedTuple(popt) #u0o is NamedTuple
end;

@testset "merge: create a mofified popt AbstractVector" begin
    popt3 = @inferred merge(popt, (k_L = 1.2,))
    @test length(popt3) == length(popt)
    @test popt3.k_L == 1.2
    @test popt3.L == popt.L
    @test popt3.k_R == popt.k_R
end;

@testset "merge: create a mofified popt AbstractVector on LVector" begin
    poptlv = LArray(popt)
    popt3 = @inferred merge(poptlv, (k_L = 1.2,))
    @test length(popt3) == length(popt)
    @test popt3.k_L == 1.2
    @test popt3.L == popt.L
    @test popt3.k_R == popt.k_R
end;

@testset "update ODEProblem optimized parameters" begin
    f = (u,p,t) -> p[1]*u
    u0 = (u1=1/2,)
    p = (p1=1.1,p2=2)
    tspan = (0.0,1.0)
    prob = ODEProblem(f,SVector(Tuple(u0)),tspan,SVector(Tuple(p)))
    #sol = solve(prob)
    #sol[end]
    ps = ProblemParSetter_sym(keys(u0),keys(p),(:u1,:p2))
    popt = (u1=1/4, p2=1.2)
    #
    # update_statepar
    prob2 = update_statepar(ps, popt, prob)
    @test prob2.u0[1] == popt.u1
    @test prob2.p[1] == p.p1 # not updated
    @test prob2.p[2] == popt.p2
    #
    # get_paropt
    @test get_paropt(ps, prob2) == collect(popt) #Tuple(popt) # Vector because provided p as SVector <: AbstractVector
    @test get_paropt_labeled(ps, prob2) == LVector(popt) # LVector because provided p as SVector <: AbstractVector
    @test label_par(ps, prob.p) == LVector(p)
end;

@testset "gradient with Tuple" begin
    fcost = (popt) -> begin
        u0, p = update_statepar(ps, popt, u1, p1)
        d = sum(get_paropt(ps, u0, p))
        d*d
    end
    fcost(popt)
    @test typeof(ForwardDiff.gradient(fcost, popt)) == typeof(popt)
end;

@testset "gradient with SVector" begin
    fcost = (popt) -> begin
        u0, p = update_statepar(ps, popt, u1, p1)
        d = sum(get_paropt(ps, u0, p))
        d*d
    end
    poptsv = SVector(popt)
    fcost(poptsv)
    @test typeof(ForwardDiff.gradient(fcost, poptsv)) == typeof(poptsv)
end;

@testset "gradient with AbstractVector" begin
    fcost = (popt) -> begin
        u0, p = update_statepar(ps, popt, u1, p1)
        d = sum(get_paropt(ps, u0, p))
        d*d
    end
    poptv = collect(popt)
    fcost(poptv)
    @test typeof(ForwardDiff.gradient(fcost, poptv)) == typeof(poptv)
end;

@testset "construct from ODESystem" begin
    @named m = samplesystem()
    popt_names = (:RHS, :τ)
    ps1 = ProblemParSetter_sym(m, popt_names)
    @test symbols_state(ps1) == (:x, :RHS)
    @test symbols_par(ps1) == (:τ,)
    @test symbols_paropt(ps1) == popt_names
    em = embed_system(m;name=:m, simplify=false)
    ps1 = ProblemParSetter_sym(em, popt_names; strip=true)
    @test symbols_state(ps1) == (:x, :RHS)
    @test symbols_par(ps1) == (:τ,)
    @test symbols_paropt(ps1) == popt_names
end;

@testset "error message on state not in front of parameters" begin
    @named m = samplesystem()
    popt_names = (:τ, :RHS) # note :RHS is a state and should be in fron to :τ
    @test_throws ErrorException ps1 = ProblemParSetter_sym(m, popt_names)
    # @test symbols_paropt(ps1) == popt_names # fails too
    #
    # in addition to missing parameters
    popt_names = (:τ, :RHS, :M) # note :RHS is a state and should be in fron to :τ
    @test_throws ErrorException ps1 = ProblemParSetter_sym(m, popt_names)
end;

@testset "name_paropt" begin
    xn = @inferred name_paropt(ps, collect(1:count_paropt(ps)))
    @test names(xn)[1] == collect(symbols_paropt(ps))
    @test xn[:k_R] == 3
    xn = @inferred name_state(ps, collect(1:count_state(ps)))
    @test names(xn)[1] == collect(symbols_state(ps))
    xn = @inferred name_par(ps, collect(1:count_par(ps)))
    @test names(xn)[1] == collect(symbols_par(ps))
    #
    frandsym = () -> begin
        syms_arr = rand([:L,:k_L,:k_R],2)
        ntuple(i -> syms_arr[i], 2)
        rand() > 0.5 ? (:L,:k_L,:k_R) : (:L,:k_L)
    end
    #frandsym()
    psr = @inferred ProblemParSetter_sym(keys(u1),keys(p1),frandsym(), Val(false))
    # cannot be inferred, because labels are not known at construction time
    #xl = @inferred label_paropt(psr, collect(1:count_paropt(psr)))
    # but NamedVector is ok
    xn = @inferred name_paropt(psr, collect(1:count_paropt(psr)))
    ftmp = () -> begin
        psr = ProblemParSetter_sym(keys(u1),keys(p1),frandsym(), Val(false))
        xn = name_paropt(psr, collect(1:count_paropt(psr)))
    end
    # @code_warntype ftmp()  # all red!! number of parameters not known
    # need different ProblemParSetter_sym that does not store integers in type signature
    # and hence will not support LabelledArrays
    _ftmp = (namesopt) -> begin
        local psr = ProblemParSetter_sym(keys(u1),keys(p1),namesopt, Val(false))
        #xn = @inferred label_paropt(psr, collect(1:count_paropt(psr)))
        xn = @inferred name_paropt(psr, collect(1:count_paropt(psr)))
    end
    ftmp = () -> begin
        local namesopt = frandsym()
        xn = @inferred _ftmp(namesopt)
        @show typeof(xn)
        # return value of entire function is not type stable, # because namesopt is not typestable
        # but inside this function typeof(xn) is known
    end
    xn = ftmp()

end;


text:
@testset "Three letters name" begin

    @test three2residue("ALA") == Residue('A')
    @test three2residue("Ala") == Residue('A')
    @test three2residue("ala") == Residue('A')

    @test_throws ErrorException three2residue("AL")
    @test_throws ErrorException three2residue("Alanine")

    @test residue2three(Residue('A')) == "ALA"

    @test_throws ErrorException residue2three(GAP)
    @test_throws ErrorException residue2three(reinterpret(Residue, -30))

    @test three2residue("XAA") == XAA
    @test three2residue("UNK") == XAA

    # Example for three letter amino acid sequence
    seq = "AlaCysAspGluPheGlyHisIleLysAsxXaaGlx"
    @test [ three2residue(seq[i:i+2]) for i in 1:3:length(seq) ] == res"ACDEFGHIKXXX"
end


text:
#http://julia-programming-language.2336112.n4.nabble.com/Re-Any-method-to-save-the-variables-in-workspace-to-file-td3123.html

# load("filename")
function loadvars(filename)
  f = open(filename, "r")
  try
    eval(parse(readall(f)))
  finally
    close(f)
  end
end

#@savevars("filename",var)
macro savevars(filename, vars...)
  printexprs = map(vars) do var
    :(print(f, ";", $(string(var)), " = "); showall(f, $(esc(var))))
  end
  quote
    local f = open($(esc(filename)), "w")
    try
      $(Expr(:block, printexprs...))
    finally
      close(f)
    end
  end
end

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
using NeuralPDE, Flux, ModelingToolkit, GalacticOptim, Optim, DiffEqFlux
using Quadrature,Cubature,Cuba
using Plots

function solve(opt)
    strategy = QuadratureTraining()

    @parameters x t
    @variables u(..)
    Dt = Differential(t)
    Dxx = Differential(x)^2

    eq = Dt(u(x,t)) - Dxx(u(x,t)) ~ -exp(-t) * (sin(pi * x) - pi^2 * sin(pi * x))

    bcs = [u(x,0) ~ sin(pi*x),
           u(-1,t) ~ 0.,
           u(1,t) ~ 0.]

    domains = [x ∈ IntervalDomain(-1.0,1.0),
               t ∈ IntervalDomain(0.0,1.0)]

    chain = FastChain(FastDense(2,18,Flux.σ),FastDense(18,18,Flux.σ),FastDense(18,1))

    discretization = PhysicsInformedNN(chain,strategy)

    indvars = [x, t]   #phisically independent variables
    depvars = [u]       #dependent (target) variable

    loss = []
    initial_time = 0

    times = []

    cb = function (p,l)
        if initial_time == 0
            initial_time = time()
        end
        push!(times, time() - initial_time)
        #println("Current loss for $opt is: $l")
        push!(loss, l)
        return false
    end

    @named pde_system = PDESystem(eq, bcs, domains, indvars, depvars)
    prob = discretize(pde_system, discretization)

    if opt == "both"
        res = GalacticOptim.solve(prob, ADAM(); cb = cb, maxiters=50)
        prob = remake(prob,u0=res.minimizer)
        res = GalacticOptim.solve(prob, BFGS(); cb = cb, maxiters=150)
    else
        res = GalacticOptim.solve(prob, opt; cb = cb, maxiters=200)
    end

    times[1] = 0.001

    return loss, times #add numeric solution
end

opt1 = ADAM()
opt2 = ADAM(0.005)
opt3 = ADAM(0.05)
opt4 = RMSProp()
opt5 = RMSProp(0.005)
opt6 = RMSProp(0.05)
opt7 = Optim.BFGS()
opt8 = Optim.LBFGS()

loss_1, times_1 = solve(opt1)
loss_2, times_2 = solve(opt2)
loss_3, times_3 = solve(opt3)
loss_4, times_4 = solve(opt4)
loss_5, times_5 = solve(opt5)
loss_6, times_6 = solve(opt6)
loss_7, times_7 = solve(opt7)
loss_8, times_8 = solve(opt8)
loss_9, times_9 = solve("both")

p = plot([times_1, times_2, times_3, times_4, times_5, times_6, times_7, times_8, times_9], [loss_1, loss_2, loss_3, loss_4, loss_5, loss_6, loss_7, loss_8, loss_9],xlabel="time (s)", ylabel="loss", xscale=:log10, yscale=:log10, labels=["ADAM(0.001)" "ADAM(0.005)" "ADAM(0.05)" "RMSProp(0.001)" "RMSProp(0.005)" "RMSProp(0.05)" "BFGS()" "LBFGS()" "ADAM + BFGS"], legend=:bottomleft, linecolor=["#2660A4" "#4CD0F4" "#FEC32F" "#F763CD" "#44BD79" "#831894" "#A6ED18" "#980000" "#FF912B"])

p = plot(1:201, [loss_1, loss_2, loss_3, loss_4, loss_5, loss_6, loss_7, loss_8, loss_9[2:end]], xlabel="iterations", ylabel="loss", yscale=:log10, labels=["ADAM(0.001)" "ADAM(0.005)" "ADAM(0.05)" "RMSProp(0.001)" "RMSProp(0.005)" "RMSProp(0.05)" "BFGS()" "LBFGS()" "ADAM + BFGS"], legend=:bottomleft, linecolor=["#2660A4" "#4CD0F4" "#FEC32F" "#F763CD" "#44BD79" "#831894" "#A6ED18" "#980000" "#FF912B"])

@show loss_1[end], loss_2[end], loss_3[end], loss_4[end], loss_5[end], loss_6[end], loss_7[end], loss_8[end], loss_9[end]

using SciMLBenchmarks
SciMLBenchmarks.bench_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])

max_stars_repo_name: SciML/SciMLBenchmarksOutput

max_issues_repo_name: SciML/SciMLBenchmarksOutput

max_forks_repo_name: SciML/SciMLBenchmarksOutput


text:
using AutoLandmarking
using Flux
using CuArrays
import ImageView

push!(LOAD_PATH, Base.@__DIR__)

# Import data with custom import functions
image_data, imname = load_imgs("C:/Users/immanueldiener/Desktop/Master/scripts/no_enamel", (128,128,128), false)
landmark_data = read_landmarks("C:/Users/immanueldiener/Desktop/Master/scripts/no_enamel", 10, "@1")

# X and Y coordinates are swapped in avizo files, wo I'll swap them back here.
#this step is essential for all the data augmentation to work.
landmark_data = swap_xy(landmark_data)

#some images have a higher value in the fill voxels than in the actual object voxels. This will fix the issue.
# The highest value must always be the filler value!
image_data[:,:,:,[109,130]] = image_data[:,:,:,[109,130]] .- 1
image_data[:,:,:,[109,130]] = image_data[:,:,:,[109,130]] .* -1

image_data2, imname2 = load_imgs("C:/Users/immanueldiener/Desktop/Master/master_data\\additional_teeth_volumes\\", (128,128,128), true)
landmark_data2 = read_landmarks("C:/Users/immanueldiener/Desktop/Master/master_data\\additional_teeth_volumes\\", 10, "@1")

#Swap X and Y again
landmark_data2 = swap_xy(landmark_data2)

#all these images have reversed fill and object values. Tha highest value must always be the filler value!
image_data2 = (image_data2 .- 1) .* -1

images = cat(image_data, image_data2, dims=4)
lms = cat(landmark_data, landmark_data2, dims=2)

#put all the landmarks that are not yet on the surface of an object onto the landmark
lms_on_surface = landmark_to_surface(images, lms, 10)

#define training and testing sets
X_train = images[:,:,:,setdiff(1:243,86:125)]
X_test = images[:,:,:,86:125]
y_train = lms_on_surface[:, setdiff(1:243,86:125)]
y_test = lms_on_surface[:, 86:125]

#augment the training data for better learning and less overfitting
X_train, y_train = mirror_vol(X_train, y_train)
flip1, lm_flipped1 = flip_volume_front(X_train, y_train)
flip2, lm_flipped2 = flip_volume_side(X_train, y_train)
flip3, lm_flipped3 = flip_volume_front(flip1, lm_flipped1)
flip4, lm_flipped4 = flip_volume_side(flip3, lm_flipped3)

# we'll be using flip2 and flip4 and rotate them by 30/60 and 20/70 degrees respectively
flip2_rot1, lm_flip2_rot1 = rotate_volumes(flip2, lm_flipped2, 30)
flip2_rot2, lm_flip2_rot2 = rotate_volumes(flip2, lm_flipped2, 60)

flip4_rot1, lm_flip4_rot1 = rotate_volumes(flip4, lm_flipped4, 20)
flip4_rot2, lm_flip4_rot2 = rotate_volumes(flip4, lm_flipped4, 70)

# lets also rotate the original data (including the mirror images)
og_rot1, lm_og_rot1 = rotate_volumes(X_train, y_train, 50)
og_rot2, lm_og_rot2 = rotate_volumes(X_train, y_train, 80)

# jittering around of the object inside the volume. Jittering is landmark based and
# a 10 Voxel padding will be added after the minimum/maximum landmark coordinate in
# any dimension.

jit1, lm_jit1 = jitter_3D(flip2_rot1, lm_flip2_rot1, 10)
jit2, lm_jit2 = jitter_3D(flip4_rot2, lm_flip4_rot2, 10)
jit3, lm_jit3 = jitter_3D(og_rot2, lm_og_rot2, 10)

#lets also use aligned versions of the original dataset
aligned, lm_aligend, reco = align_principal(X_train, y_train, 192)
resized, lm_resized = resize_relevant(aligned, lm_aligend, 128)
aligned = nothing

jit4, lm_jit4 = jitter_3D(X_train, y_train, 10)

#generate all the depthmaps for input into the network
complete1 = depth_map_all_sides(flip2_rot1)
complete2 = depth_map_all_sides(flip2_rot2)
complete3 = depth_map_all_sides(flip4_rot1)
complete4 = depth_map_all_sides(flip4_rot2)
complete5 = depth_map_all_sides(og_rot1)
complete6 = depth_map_all_sides(og_rot2)
complete7 = depth_map_all_sides(jit1)
complete8 = depth_map_all_sides(jit2)
complete9 = depth_map_all_sides(jit3)
complete10 = depth_map_all_sides(resized)
complete11 = depth_map_all_sides(flip1)
complete12 = depth_map_all_sides(flip3)
complete13 = depth_map_all_sides(jit4)

X_train = cat(complete1, complete2, complete3, complete4, complete5, complete6, complete7, complete8,
    complete9, complete10, complete11, complete12, complete13, dims=4)
y_train = cat(lm_flip2_rot1, lm_flip2_rot2, lm_flip4_rot1, lm_flip4_rot2, lm_og_rot1, lm_og_rot2,
    lm_jit1, lm_jit2, lm_jit3, lm_resized, lm_flipped1, lm_flipped3, lm_jit4, dims=2)

# make depthmaps of the testing data
X_test = depth_map_all_sides(X_test)

# define the cost function
cost(x, y) = sum((model(x)-y).^2)|>gpu

#define the model
model = Flux.mapleaves(cu, AutoLandmarking.vgg19)

# define the trainingrate and optimiser
opt = Flux.ADAM(0.000015)

# redefinition of the dropout function to work with testmode!() since this
# functionality is not working at the moment. This step only needs to be done
# if you wish to perform uncertainty estimation using the dropout method (AutoLM.response_distribution())
# after training.
using Random
function Flux.dropout(x, p; dims = :)
    q = 1 - p
    y = rand!(similar(x, Flux._dropout_shape(x, dims)))
    y .= Flux._dropout_kernel.(y, p, q)
    x .* y
end

import Zygote
Zygote.@adjoint function Flux.dropout(x, p; dims = :)
   q = 1 - p
   y = rand!(similar(x, Flux._dropout_shape(x, dims)))
   y .= Flux._dropout_kernel.(y, p, q)
   return x .* y, Δ -> (Δ .* y, nothing)
end
Zygote.refresh()

function run_model(modell, X, y)
  train_data = Flux.mapleaves(cu, X)
  train_labels = Flux.mapleaves(cu, y)
  dataset = Flux.Data.DataLoader(train_data, train_labels, batchsize = 4, shuffle=true)
  Flux.train!(cost, params(model), dataset, opt)
  testmode!(model)
  cosima = cost_whole_data_2D(train_data, train_labels, cost)
  testmode!(model, false)
  return cosima
end

# create lists to store metrics
accs = []
costs = []

# run the training for 300 epochs. batches of 128 will be loaded onto the GPU,
# which will be further subdevided into minibatches of 4 (as defined in AutoML.run_model).
for i in 1:300
  costr = 0
  for j in 1:128:size(X_train, 4)
    if j+127 > size(X_train, 4)
      train_data = X_train[:,:,:,j:end]
      train_labels = y_train[:,j:end]
    else
      train_data = X_train[:,:,:,j:j+127]
      train_labels = y_train[:,j:j+127]
    end
    cosima = run_model(model, train_data, train_labels)
    costr = costr + cosima
    train_data = nothing
    train_labels = nothing
    # GC.gc()
  end
  push!(costs, costr)
  println("epoch ", i, " finished")
  println(costr)
  testmode!(model)
  acc2, max1, min1 = AutoLandmarking.avg_accuracy_per_point(model, gpu(X_test[:,:,:,:]), gpu(y_test[:,:]), 3)
  testmode!(model, false)
  println("median deviation per point on testing dataset: ", acc2, "with maximum", max1)
  println("median sum is ", sum(acc2))
  push!(accs, sum(acc2))
  if sum(acc2)[1]<7.5
    break
  end
end

#save model
using BSON
model_cpu = cpu(model)
BSON.@save "C:/Users/immanueldiener/Desktop/Master/master_data/model_demo.bson" model_cpu

#load model back in
BSON.@load "C:/Users/immanueldiener/Desktop/Master/master_data/model_demo.bson" model_cpu
model = gpu(model_cpu)

# ================== After training ===================

#solve issue with arrays inside vector
aks = []
for i in accs
  push!(aks, i[1])
end

# check what the best accuracy is
minimum(aks)

#plot the development of cost and accuracy
import Plots
Plots.plot(aks, legend=:topright, label="sum of deviations per point", color= :red, xlabel="epochs", ylabel = "acuracy")
plt = Plots.plot!(Plots.twinx(), costs, label="training loss", legend=:bottomright, ylabel = "loss")

# check for outliers (greatest distance from the procrustes mean).
# teh idea is to eliminate out all the predicted landmark sets that
# deviate a lot from the procrustes mean (reference), since these are
# most likely bad predictions. do get the procrustes mean, we can use training
# testing data together to get a more accureate mean. This is an ideal approch
# because otherwise we would not realize if all the predictions are biased in some
# way that is congruent over all predictions but still wrong.


text:
module Mapper2

const SRCDIR = @__DIR__
const PKGDIR = dirname(SRCDIR)
using Reexport

using Logging

# Aux modules
include("Helper.jl")
using .Helper

include("MapperGraphs.jl")

# Central data types.
include("MapperCore.jl")

# Placement modules
include("Place/SA/SA.jl")

# Routing modules
include("Route/Route.jl")

# exports from Helper.
# TODO: DO we really need these exports?
export Address, dim_max, dim_min, place!

# Use submodules to make exports visible.
@reexport using .MapperCore

#############
# PLACEMENT #
#############
# Placement Algorithms
@reexport using .SA

# Default Placement Algorithm
place!(map::Map; kwargs...) = SA.place!(map; kwargs...)

###########
# ROUTING #
###########
@reexport using .Routing

end #module Mapper2


text:
struct DockerExecutor <: SandboxExecutor
end

Base.show(io::IO, exe::DockerExecutor) = write(io, "Docker Executor")

function executor_available(::Type{DockerExecutor}; verbose::Bool = false)
    # don't even try to exec if it doesn't exist
    if Sys.which("docker") === nothing
        if verbose
            @info("No `docker` command found; docker unavailable")
        end
        return false
    end

    # Return true if we can `docker ps`; if we can't, then there's probably a permissions issue
    if !success(`docker ps`)
        if verbose
            @warn("Unable to run `docker ps`; perhaps you're not in the `docker` group?")
        end
        return false
    end
    return probe_executor(DockerExecutor(); test_read_only_map=true, test_read_write_map=true, verbose)
end

timestamps_path() = joinpath(@get_scratch!("docker_timestamp_hashes"), "path_timestamps.toml")
function load_timestamps()
    path = timestamps_path()
    if !isfile(path)
        return Dict()
    end
    try
        return TOML.parsefile(path)
    catch e
        @error("couldn't load $(path)", exception=e)
        return Dict()
    end
end

function save_timestamp(image_name::String, timestamp::Float64)
    open(timestamps_path(), read=true, write=true, create=true) do io
        timestamps = try
            TOML.parse(io)
        catch
            return Dict()
        end
        seekstart(io)
        timestamps[image_name] = timestamp
        TOML.print(io, timestamps)
    end
end

docker_image_name(root_path::String) = "sandbox_rootfs:$(string(Base._crc32c(root_path), base=16))"
function should_build_docker_image(root_path::String)
    # If this image has been built before, compare its historical timestamp to the current one
    curr_ctime = max_directory_ctime(root_path)
    prev_ctime = get(load_timestamps(), docker_image_name(root_path), 0.0)
    return curr_ctime != prev_ctime
end

"""
    build_docker_image(root_path::String)

Docker doesn't like volume mounts within volume mounts, like we do with `sandbox`.
So we do things "the docker way", where we construct a rootfs docker image, then mount
things on top of that, with no recursive mounting.  We cut down on unnecessary work
somewhat by quick-scanning the directory for changes and only rebuilding if changes
are detected.
"""
function build_docker_image(root_path::String; verbose::Bool = false)
    image_name = docker_image_name(root_path)
    if should_build_docker_image(root_path)
        max_ctime = max_directory_ctime(root_path)
        if verbose
            @info("Building docker image $(image_name) with max timestamp $(max_ctime)")
        end

        # Build the docker image
        open(`docker import - $(image_name)`, "w", verbose ? stdout : devnull) do io
            Tar.create(root_path, io)
        end

        # Record that we built it
        save_timestamp(image_name, max_ctime)
    end

    return image_name
end

function build_executor_command(exe::DockerExecutor, config::SandboxConfig, user_cmd::Cmd)
    # Begin building `docker` args
    cmd_string = String["docker", "run", "--rm", "--privileged", "-i"]

    # If we're doing a fully-interactive session, tell it to allocate a psuedo-TTY
    if all(isa.((config.stdin, config.stdout, config.stderr), Base.TTY))
        push!(cmd_string, "-t")
    end

    # Build the docker image that corresponds to this rootfs
    image_name = build_docker_image(config.read_only_maps["/"])

    # Start in the right directory
    append!(cmd_string, ["-w", config.pwd])

    # Add in read-only mappings (skipping the rootfs)
    for (dst, src) in config.read_only_maps
        if dst == "/"
            continue
        end
        append!(cmd_string, ["-v", "$(src):$(dst):ro"])
    end

    # Add in read-write mappings
    for (dst, src) in config.read_write_maps
        append!(cmd_string, ["-v", "$(src):$(dst)"])
    end

    # Apply environment mappings, first from `config`, next from `user_cmd`.
    for (k, v) in config.env
        append!(cmd_string, ["-e", "$(k)=$(v)"])
    end
    if user_cmd.env !== nothing
        for pair in user_cmd.env
            append!(cmd_string, ["-e", pair])
        end
    end

    # Finally, append the docker image name user-requested command string
    push!(cmd_string, image_name)
    append!(cmd_string, user_cmd.exec)

    docker_cmd = Cmd(cmd_string)

    # If the user has asked that this command be allowed to fail silently, pass that on
    if user_cmd.ignorestatus
        docker_cmd = ignorestatus(docker_cmd)
    end

    return docker_cmd
end


text:
struct SpringForce <: AbstractForce
    A :: Real
end

function potential(f::SpringForce, x::Real, t::Real)::Real
    - f.A * x^2
end

function virtual_workrate(f::SpringForce, x::Real,v::Real, t::Real)::Real
    total_force(f, x, t) * v
end

function erased(f::SpringForce, x::Real)::Bool
    x < 0
end