text:
StaggeredGridCalculus# The following without the qualifier StaggeredGridCalculus cannot find the Grid constructor
# of StaggeredGridCalculus:
# Grid(lprim::AbsVecReal, Npml::NTuple{2,Int}, isbloch::Bool) =
#     StaggeredGridCalculus.Grid((lprim,), ([Npml[nN]], [Npml[nP]]), (isbloch,))

# Calculate ghost points from l, L, and isbloch.
lghost(l::NTuple{2,NTuple{K,AbstractVector{<:Real}}},  # grid point locations
       L::SVec{K,Float64},  # domain size
       isbloch::SVec{K,Bool}  # boundary condition
       ) where {K} =
    (map((lprimₖ,Lₖ) -> lprimₖ[1]+Lₖ, SVec(l[nPR]), L),  # lg[PRIM]
     map((isblochₖ,lprimₖ,ldualₖ,Lₖ) -> (isblochₖ ? ldualₖ[end]-Lₖ : 2lprimₖ[1]-ldualₖ[1]), isbloch, SVec(l[nPR]), SVec(l[nDL]), L))  # lg[DUAL]

@testset "grid" begin

@test isa(Grid(tuple(1:10), [true]), Any)

@testset "Grid{1}, Bloch boundary" begin
    isbloch = [true]
    N = 22
    ∆ldual = rand(N)
    L = sum(∆ldual)
    l₀ = L / 2
    lprim = cumsum([-l₀; ∆ldual])

    g1 = Grid(tuple(lprim), isbloch)

    @test g1.N == [N]
    @test g1.L ≈ [L]
    @test all(SVec.(lprim) .∈ Ref(g1))  # `∈` supports only SVec, not Vector
    ldual = StaggeredGridCalculus.movingavg(lprim)
    pop!(lprim)
    @test g1.l ≈ ((lprim,), (ldual,))
    ∆lprim = [ldual[1]+L-ldual[end]; diff(ldual)]
    @test g1.∆l ≈ ((∆lprim,), (∆ldual,))
    @test g1.isbloch == isbloch
    # @test g1.Npml == ([Npmln], [Npmlp])
    # @test g1.Lpml ≈ ([sum(∆ldual[1:Npmln])], [sum(∆ldual[end-Npmlp+1:end])])
    # @test g1.lpml ≈ ([sum(∆ldual[1:Npmln])-l₀], [sum(∆ldual)-sum(∆ldual[end-Npmlp+1:end])-l₀])
    @test g1.bounds ≈ ([lprim[1]], [lprim[end]+∆ldual[end]])
    @test SVec(-l₀) ∈ g1  # `∈` supports only SVec
    lg = lghost(((lprim,),(ldual,)), SVec(L), SVec{1}(isbloch))
    lprim_g = g1.ghosted.l[nPR][1]
    ldual_g = g1.ghosted.l[nDL][1]
    @test lprim_g[1:end-1]≈lprim && lprim_g[end]≈lg[nPR][1] && ldual_g[2:end]≈ldual && ldual_g[1]≈lg[nDL][1]
    τlprim_g = g1.ghosted.τl[nPR][1]
    τldual_g = g1.ghosted.τl[nDL][1]
    @test τlprim_g[1:end-1]≈lprim && τlprim_g[end]≈lg[nPR][1]-L && τldual_g[2:end]≈ldual && τldual_g[1]≈lg[nDL][1]+L
    τindprim_g = g1.ghosted.τind[nPR][1]
    τinddual_g = g1.ghosted.τind[nDL][1]
    @test all(τindprim_g[1:end-1].==1:N) && τindprim_g[end]==1 && all(τinddual_g[2:end].==2:N+1) && τinddual_g[1]==N+1
    ∆τprim_g = g1.ghosted.∆τ[nPR][1]
    ∆τdual_g = g1.ghosted.∆τ[nDL][1]
    @test iszero(∆τprim_g[1:end-1]) && ∆τprim_g[end]≈-L && iszero(∆τdual_g[2:end]) && ∆τdual_g[1]≈L
end  # @testset "Grid{1}, primal boundary"

@testset "Grid{3}" begin
    N = [29, 23, 6]
    ∆ldual = ntuple(d->rand(N[d]), length(N))
    isbloch = [false, false, true]

    L = SVec(sum.(∆ldual))  # SFloat{3}
    l₀ = L ./ 2  # SFloat{3}
    lprim = map((v,s)->v.-s, map(x->[0; cumsum(x)], ∆ldual), (l₀...,))  # tuple of vectors

    g3 = Grid(lprim, isbloch)

    @test g3.N == N
    @test g3.L ≈ L
    ldual = StaggeredGridCalculus.movingavg.(lprim)
    pop!.(lprim)
    @test g3.l ≈ (lprim, ldual)
    ∆lprim = diff.(ldual)
    prepend!(∆lprim[1], 2(ldual[1][1]-lprim[1][1]))
    prepend!(∆lprim[2], 2(ldual[2][1]-lprim[2][1]))
    prepend!(∆lprim[3], ldual[3][1]+L[3]-ldual[3][end])
    @test g3.∆l ≈ (∆lprim, ∆ldual)
    @test g3.isbloch == isbloch
    # @test g3.Npml == Npml
    # @test g3.Lpml ≈ (
    #     [sum(∆ldual[d][1:Npml[nN][d]]) for d = 1:3],
    #     [sum(∆ldual[d][end-Npml[nP][d]+1:end]) for d = 1:3],
    # )
    # @test g3.lpml ≈ (
    #     [sum(∆ldual[d][1:Npml[nN][d]]) - l₀[d] for d = 1:3],
    #     [sum(∆ldual[d]) - sum(∆ldual[d][end-Npml[nP][d]+1:end]) - l₀[d] for d = 1:3]
    # )
    @test g3.bounds ≈ (
        [lprim[d][1] for d = 1:3],
        [lprim[d][end] + ∆ldual[d][end] for d = 1:3]
    )
    @test -l₀ ∈ g3
    @test all(g3.bounds .∈ Ref(g3))
    lg = lghost((lprim,ldual), L, g3.isbloch)
    lprim_g = g3.ghosted.l[nPR]
    ldual_g = g3.ghosted.l[nDL]
    @test pop!.(lprim_g)≈lg[nPR].data && popfirst!.(ldual_g)≈lg[nDL].data && lprim_g≈lprim && ldual_g≈ldual
    τlprim_g = g3.ghosted.τl[nPR]
    τldual_g = g3.ghosted.τl[nDL]
    @test pop!.(τlprim_g)≈lg[nPR].data.-(0,0,L[3]) && popfirst!.(τldual_g)≈(ldual[1][1],ldual[2][1],lg[nDL][3]+L[3]) && τlprim_g≈lprim && τldual_g≈ldual
    τindprim_g = g3.ghosted.τind[nPR]
    τinddual_g = g3.ghosted.τind[nDL]
    @test pop!.(τindprim_g)==(N[1]+1,N[2]+1,1) && popfirst!.(τinddual_g)==(2,2,N[3]+1) && τindprim_g==map(m->collect(1:m), tuple(N...)) && τinddual_g==map(m->collect(2:m+1), tuple(N...))
    ∆τprim_g = g3.ghosted.∆τ[nPR]
    ∆τdual_g = g3.ghosted.∆τ[nDL]
    @test pop!.(∆τprim_g)≈(0,0,-L[3]) && popfirst!.(∆τdual_g)≈(0,0,L[3]) && all(iszero.(∆τprim_g)) && all(iszero.(∆τdual_g))
end  # @testset "Grid{3}"

# @testset "Grid{2}" begin
#     isbloch = (false, false, true)
#     Npml = ((10,9), (5,6), (1,4))
#     N = (10, 12, 1)
#     M = (+).(sum.(Npml), N)
#     ∆ldual = rand.(M)  # tuple of arrays
#     L = sum.(∆ldual)
#     l₀ = (/).(L, 2)
#     lprim = (-).(map(x->[0; cumsum(x)], ∆ldual), l₀)
#
#     g3 = Grid3D(lprim, Npml, isbloch)
#     normal_axis = Ŷ
#     h, v = Int.(next3(normal_axis))
#     g2 = Grid2D(g3, normal_axis)
#
#     lprim = (lprim[h], lprim[v])
#     ∆ldual = (∆ldual[h], ∆ldual[v])
#     M = (M[h], M[v])
#     L = (L[h], L[v])
#     isbloch = (isbloch[h], isbloch[v])
#     Npml = (Npml[h], Npml[v])
#     l₀ = (l₀[h], l₀[v])
#
#     @test g2.normal_axis == Ŷ
#     @test g2.N == M
#     @test g2.L ≈ L
#     @test bounds_(g2) ≈ ((lprim[nHRZ][1],lprim[nHRZ][end]), (lprim[nVRT][1],lprim[nVRT][end]))
#     zdual = (lprim[nHRZ][1:end-1] + lprim[nHRZ][2:end]) / 2
#     xdual = (lprim[nVRT][1:end-1] + lprim[nVRT][2:end]) / 2
#     ldual = (zdual, xdual)
#     lprim = (nw->lprim[nw][1:end-1]).(nPD)
#     @test l_(g2) ≈ (lprim, ldual)
#     ∆lprim = (
#         [ldual[nHRZ][1]+L[nHRZ]-ldual[1][end]; diff(zdual)],
#         [2(ldual[nVRT][1]-lprim[nVRT][1]); diff(xdual)]
#     )
#     @test ∆l_(g2) ≈ (∆lprim, ∆ldual)
#     @test isbloch_(g2) == isbloch
#     @test Npml_(g2) == Npml
#     @test Lpml_(g2) ≈ (
#         (sum(∆ldual[nHRZ][1:Npml[nHRZ][nN]]), sum(∆ldual[nHRZ][end-Npml[nHRZ][nP]+1:end])),
#         (sum(∆ldual[nVRT][1:Npml[nVRT][nN]]), sum(∆ldual[nVRT][end-Npml[nVRT][nP]+1:end]))
#     )
#     @test lpml_(g2) ≈ (
#         (sum(∆ldual[nHRZ][1:Npml[nHRZ][nN]]) - l₀[nHRZ], sum(∆ldual[nHRZ]) - sum(∆ldual[nHRZ][end-Npml[nHRZ][nP]+1:end]) - l₀[nHRZ]),
#         (sum(∆ldual[nVRT][1:Npml[nVRT][nN]]) - l₀[nVRT], sum(∆ldual[nVRT]) - sum(∆ldual[nVRT][end-Npml[nVRT][nP]+1:end]) - l₀[nVRT])
#     )
#     @test contains(g2, map(-,l₀))
#     @test all(.!contains.(g2, [bounds_(g2,NEG), bounds_(g2,POS)], false))
#     @test all(contains.(g2, [bounds_(g2,NEG), bounds_(g2,POS)]))
#     # @test pbounds(g2, withpml=true) ≈ ((-l₀[1],sum(∆ldual[1])-l₀[1]), (-l₀[2],sum(∆ldual[2])-l₀[2]))
#     # @test pbounds(g2, withpml=false) ≈ lpml_(g2)
#     # @test lplot(g2, withinterp=true, withpml=true) ≈ (
#     #     (lprim[1],[lprim[1][1];ldual[1][2:end-1];lprim[1][end]]),
#     #     (lprim[2],[lprim[2][1];ldual[2][2:end-1];lprim[2][end]])
#     # )
#     # @test voxel_edges(g2, withpml=true) ≈ ((ldual[1],lprim[1]), (ldual[2],lprim[2]))
# end  # @testset "Grid{2}"

@testset "isproper_blochphase" begin
    @test isproper_blochphase(1im, true)
    @test isproper_blochphase(cos(π/0.9)+sin(π/0.9)im, true)
    @test !isproper_blochphase(2., true)
    @test !isproper_blochphase(2., false)
end

end  # @testset "grid"


text:
module Axis

using Printf
using Markdown
using Cairo
using ComplexVisual
@ComplexVisual.import_huge

import Main.DocGenerator: DocSource, DocCreationEnvironment, DocContext,
        Document, substitute_marker_in_markdown, create_doc_icon, append_md

"""
# [![./Axis_docicon.png]({image_from_canvas: get_doc_icon()}) Axes](./Axis.md)

## Quicklinks

`cv_format_ticks`   `CV_TickLabelAppearance`   `CV_Ruler`
`↦:Tuple{CV_TickLabelAppearance, Tuple{Vararg{CV_TickLabel{LocT}, N}} where {N, LocT}}`
`cv_ticks_labels`
`⇒:Tuple{AbstractString, Real}`(label, location)
`⇒:Tuple{AbstractString, AbstractVector{T} where T}`(format, locations)

## Names

Let's start with some vocabulary

* ![./Axis_ticks.png]({image_from_canvas: vocab_tick()})
  tick: A marker for a mathematical coordinate
* ![./Axis_ticklabels.png]({image_from_canvas: vocab_ticklabel()})
  tick label: A text below a tick (may be the empty string, i.e. no text)
* ![./Axis_ruler.png]({image_from_canvas: vocab_ruler()})
  ruler: several ticks (with optional tick labels) and information about
  styling
* ![./Axis_axis.png]({image_from_canvas: vocab_axis()})
  axis: several rulers
"""
axes_intro() = nothing

"""
## Tick labels

### `doc: CV_TickLabel`

### `doc: ⇒:Tuple{AbstractString, Real}`

### `doc: cv_format_ticks`

### `doc: ⇒:Tuple{AbstractString, AbstractVector{T} where T}`

### `doc: CV_TickLabelAppearance`
"""
help_ticklabel() = nothing

"""
### `doc: CV_Ruler`

### `doc: ↦:Tuple{CV_TickLabelAppearance, Tuple{Vararg{CV_TickLabel{LocT}, N}} where {N, LocT}}`

### `doc: cv_create_2daxis_canvas`

### `doc: cv_ticks_labels`
"""
help_rulers() = nothing

function vocab_basic(format_string="", app=CV_TickLabelAppearance())
    layout = CV_2DLayout()

    math_canvas = CV_Math2DCanvas(0.0 +1.0im, 8.0 + 0.0im, 30)

    ticks = cv_format_ticks(format_string, 0.0:1.0:8.0...)
    rulers = (CV_Ruler(ticks, app), )

    axis_canvas = cv_create_2daxis_canvas(math_canvas, cv_south, rulers)
    axis_canvas_l = cv_add_canvas!(layout, axis_canvas, (0,0), (0,0))

    cv_add_padding!(layout, 5)
    can_layout = cv_canvas_for_layout(layout)
    cv_create_context(can_layout) do con_layout
        axis_canvas_l(con_layout)
    end

    return can_layout
end

vocab_tick() = vocab_basic("", CV_TickLabelAppearance())
vocab_ticklabel() = vocab_basic("%.0f", CV_TickLabelAppearance())
vocab_ruler() = vocab_basic("%.0f", CV_TickLabelAppearance(;
    tick_style=cv_color(0,0,.5) → cv_linewidth(2),
    label_style=cv_color(0,.5,0) → cv_fontface("serif") → cv_fontsize(15)))

function vocab_axis()
    layout = CV_2DLayout()

    math_canvas = CV_Math2DCanvas(0.0 +1.0im, 8.0 + 0.0im, 30)

    loc1 = 0.0:2.0:8.0
    ticks1 = cv_format_ticks("%0.f", loc1...)
    app1 = CV_TickLabelAppearance(;
        tick_style=cv_color(0,0,.5) → cv_linewidth(2),
        label_style=cv_color(0,.5,0) → cv_fontface("serif") → cv_fontsize(15))

    loc2 = 1.0:2.0:7.0
    ticks2 = cv_format_ticks("", loc2...)
    app2 = CV_TickLabelAppearance(; tick_length=7)

    loc3 = setdiff(setdiff(0.0:0.25:8.0, loc1), loc2)
    ticks3 = cv_format_ticks("", loc3...)
    app3 = CV_TickLabelAppearance(; tick_length=4)

    rulers = (CV_Ruler(ticks1, app1), CV_Ruler(ticks2, app2),
              CV_Ruler(ticks3, app3), )

    axis_canvas = cv_create_2daxis_canvas(math_canvas, cv_south, rulers)
    axis_canvas_l = cv_add_canvas!(layout, axis_canvas, (0,0), (0,0))

    cv_add_padding!(layout, 5)
    can_layout = cv_canvas_for_layout(layout)
    cv_create_context(can_layout) do con_layout
        axis_canvas_l(con_layout)
    end

    return can_layout
end

"""
## Advanced Example

Here is an example showing the construction of rulers.

![./Axis_cmin_example.png]({image_from_canvas: example_inches_cm()})

```julia
{func: example_inches_cm}
```
"""
function example_inches_cm()
    layout = CV_2DLayout()
    math_canvas = CV_Math2DCanvas(0.0 +0.0im, 1.0 - 8.0im, 60)

    in_unit = 2.54
    col_cm, col_in = cv_color(0, 0, 0.5), cv_color(0, 0.5, 0)
    font_large = cv_fontface("serif") → cv_fontsize(15)
    font_small = cv_fontface("serif") → cv_fontsize(10)

    app1(col) = CV_TickLabelAppearance(; tick_length=30, gap=7,
        tick_style=col → cv_linewidth(2), label_style=col → font_large)
    app2(col) = CV_TickLabelAppearance(; tick_length=15,
        tick_style=col → cv_linewidth(2), label_style=col → font_large)
    app3(col) = CV_TickLabelAppearance(; tick_length=9,
        tick_style=col → cv_linewidth(2), label_style=col → font_small)
    app4(col) = CV_TickLabelAppearance(; tick_length=0, gap=5,
        tick_style=col → cv_linewidth(2), label_style=col → font_large)

    loc1 = 0.0:1.0:8.0;  ticks1 = cv_format_ticks("%0.f", loc1...)
    loc2 = 0.5:1.0:7.5;  ticks2 = cv_format_ticks("", loc2...)
    loc3 = setdiff(setdiff(0.0:0.1:8.0, loc1), loc2)

    ticks3 = cv_format_ticks("", loc3...)

    rulers1 = (
        CV_Ruler(ticks1, app1(col_cm)), CV_Ruler(ticks2, app2(col_cm)),
        CV_Ruler(ticks3, app3(col_cm)),)
    axis1_canvas = cv_create_2daxis_canvas(math_canvas, cv_west, rulers1)

    loc1 = 0.0:1.0:3.0
    ticks1 = tuple(
        map( v -> CV_TickLabel(v*in_unit, @sprintf("%0.f", v)), loc1)...)
    loc2 = setdiff(0.0:0.5:2.5, loc1)
    ticks2 = tuple(map(v -> CV_TickLabel(v*in_unit, "½"), loc2)...)
    loc3 = setdiff(setdiff(0.0:0.25:3.0, loc1), loc2)
    ticks3 = tuple(map(e -> CV_TickLabel(e[2]*in_unit,
            isodd(e[1]) ? "¼" : "¾"), enumerate(loc3))...)

    rulers2 = (
        CV_Ruler(ticks1, app1(col_in)), CV_Ruler(ticks2, app2(col_in)),
        CV_Ruler(ticks3, app3(col_in)),)
    axis2_canvas = cv_create_2daxis_canvas(math_canvas, cv_east, rulers2)

    axis1_canvas_l = cv_add_canvas!(layout, axis1_canvas,
        cv_anchor(axis1_canvas, :default), (0,0))
    axis2_canvas_l = cv_add_canvas!(layout, axis2_canvas,
        cv_anchor(axis2_canvas, :default), (0,0))

    cm_text = cv_text("cm", col_cm → font_large)
    cm_l = cv_add_canvas!(layout, cm_text,
        cv_anchor(cm_text, :baseline_center),
        cv_translate(cv_anchor(axis1_canvas_l, :north), 0, -20))

    in_text = cv_text("inch", col_in → font_large)
    in_l = cv_add_canvas!(layout, in_text,
        cv_anchor(in_text, :baseline_center),
        cv_translate(
            cv_anchor(axis2_canvas_l, :north, axis1_canvas_l, :north), 0, -20))

    cv_add_padding!(layout, 10)

    can_layout = cv_canvas_for_layout(layout)
    cv_create_context(can_layout) do con_layout
        axis1_canvas_l(con_layout); axis2_canvas_l(con_layout)
        cm_l(con_layout); in_l(con_layout)
    end

    return can_layout
end

function get_doc_icon()
    src_canvas = example_inches_cm()
    icon = create_doc_icon(src_canvas, cv_rect_blwh(Int32, 133, -50, 200, 200))
    return icon
end

function create_document(doc_env::DocCreationEnvironment)
    doc_source = DocSource("Axis", @__MODULE__)
    context = DocContext(doc_env, doc_source)

    md = Markdown.MD()
    for part in (axes_intro, help_ticklabel, help_rulers, example_inches_cm)
        part_md = Base.Docs.doc(part)
        substitute_marker_in_markdown(context, part_md)

        append_md(md, part_md)
    end

    doc = Document(doc_source, md)
    return doc
end

end

# vim:syn=julia:cc=79:fdm=marker:sw=4:ts=4:

ext: jl


text:
type CurrentPlot
    nullableplot::Nullable{AbstractPlot}
end
const CURRENT_PLOT = CurrentPlot(Nullable{AbstractPlot}())

isplotnull() = isnull(CURRENT_PLOT.nullableplot)

function current()
    if isplotnull()
        error("No current plot/subplot")
    end
    get(CURRENT_PLOT.nullableplot)
end
current(plot::AbstractPlot) = (CURRENT_PLOT.nullableplot = Nullable(plot))

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

Base.string(plt::Plot) = "Plot{$(plt.backend) n=$(plt.n)}"
Base.print(io::IO, plt::Plot) = print(io, string(plt))
Base.show(io::IO, plt::Plot) = print(io, string(plt))

getplot(plt::Plot) = plt
getattr(plt::Plot, idx::Int = 1) = plt.attr
convertSeriesIndex(plt::Plot, n::Int) = n

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

"""
The main plot command.  Use `plot` to create a new plot object, and `plot!` to add to an existing one:

```
    plot(args...; kw...)                  # creates a new plot window, and sets it to be the current
    plot!(args...; kw...)                 # adds to the `current`
    plot!(plotobj, args...; kw...)        # adds to the plot `plotobj`
```

There are lots of ways to pass in data, and lots of keyword arguments... just try it and it will likely work as expected.
When you pass in matrices, it splits by columns.  See the documentation for more info.
"""

# this creates a new plot with args/kw and sets it to be the current plot
function plot(args...; kw...)
    d = KW(kw)
    preprocessArgs!(d)

    # create an empty Plot then process
    plt = Plot()
    # plt.user_attr = d
    _plot!(plt, d, args...)
end

# build a new plot from existing plots
# note: we split into plt1 and plts_tail so we can dispatch correctly
function plot(plt1::Plot, plts_tail::Plot...; kw...)
    d = KW(kw)
    preprocessArgs!(d)

    # build our plot vector from the args
    n = length(plts_tail) + 1
    plts = Array(Plot, n)
    plts[1] = plt1
    for (i,plt) in enumerate(plts_tail)
        plts[i+1] = plt
    end

    # compute the layout
    layout = layout_args(d, n)[1]
    num_sp = sum([length(p.subplots) for p in plts])

    # create a new plot object, with subplot list/map made of existing subplots.
    # note: we create a new backend figure for this new plot object
    # note: all subplots and series "belong" to this new plot...
    plt = Plot()

    # TODO: build the user_attr dict by creating "Any matrices" for the args of each subplot

    # TODO: replace this with proper processing from a merged user_attr KW
    # update plot args, first with existing plots, then override with d
    for p in plts
        _update_plot_args(plt, p.attr)
        plt.n += p.n
    end
    _update_plot_args(plt, d)

    # pass new plot to the backend
    plt.o = _create_backend_figure(plt)
    plt.init = true

    # create the layout and initialize the subplots
    plt.layout, plt.subplots, plt.spmap = build_layout(layout, num_sp, copy(plts))
    for (idx, sp) in enumerate(plt.subplots)
        _initialize_subplot(plt, sp)
        serieslist = series_list(sp)
        if sp in sp.plt.inset_subplots
            push!(plt.inset_subplots, sp)
        end
        sp.plt = plt
        sp.attr[:subplot_index] = idx
        for series in serieslist
            push!(plt.series_list, series)
            _series_added(plt, series)
        end
    end

    # # just in case the backend needs to set up the plot (make it current or something)
    # _prepare_plot_object(plt)

    # first apply any args for the subplots
    for (idx,sp) in enumerate(plt.subplots)
        _update_subplot_args(plt, sp, d, idx, remove_pair = false)
    end

    # # now we can get rid of the axis keys without a letter
    # for k in keys(_axis_defaults)
    #     delete!(d, k)
    #     for letter in (:x,:y,:z)
    #         delete!(d, Symbol(letter,k))
    #     end
    # end

    # do we need to link any axes together?
    link_axes!(plt.layout, plt[:link])

    # finish up
    current(plt)
    if get(d, :show, default(:show))
        gui()
    end
    plt
end

# this adds to the current plot, or creates a new plot if none are current
function  plot!(args...; kw...)
    local plt
    try
        plt = current()
    catch
        return plot(args...; kw...)
    end
    plot!(current(), args...; kw...)
end

# this adds to a specific plot... most plot commands will flow through here
function plot!(plt::Plot, args...; kw...)
    d = KW(kw)
    preprocessArgs!(d)
    # merge!(plt.user_attr, d)
    _plot!(plt, d, args...)
end

function strip_first_letter(s::Symbol)
    str = string(s)
    str[1:1], Symbol(str[2:end])
end

# this method recursively applies series recipes when the seriestype is not supported
# natively by the backend
function _apply_series_recipe(plt::Plot, d::KW)
    # replace seriestype aliases
    st = d[:seriestype]
    st = d[:seriestype] = get(_typeAliases, st, st)

    # if it's natively supported, finalize processing and pass along to the backend, otherwise recurse
    if st in supported_types()

        # getting ready to add the series... last update to subplot from anything
        # that might have been added during series recipes
        sp = d[:subplot]
        sp_idx = get_subplot_index(plt, sp)
        _update_subplot_args(plt, sp, d, sp_idx)

        # do we want to override the series type?
        if !is3d(st) && d[:z] != nothing && (size(d[:x]) == size(d[:y]) == size(d[:z]))
            st = d[:seriestype] = (st == :scatter ? :scatter3d : :path3d)
        end

        # change to a 3d projection for this subplot?
        if is3d(st)
            sp.attr[:projection] = "3d"
        end

        # initialize now that we know the first series type
        if !haskey(sp.attr, :init)
            _initialize_subplot(plt, sp)
            sp.attr[:init] = true
        end

        # strip out series annotations (those which are based on series x/y coords)
        # and add them to the subplot attr
        sp_anns = annotations(sp[:annotations])
        anns = annotations(pop!(d, :series_annotations, []))
        if length(anns) > 0
            x, y = d[:x], d[:y]
            nx, ny, na = map(length, (x,y,anns))
            n = max(nx, ny, na)
            anns = [(x[mod1(i,nx)], y[mod1(i,ny)], text(anns[mod1(i,na)])) for i=1:n]
        end
        sp.attr[:annotations] = vcat(sp_anns, anns)

        # adjust extrema and discrete info
        if st == :image
            w, h = size(d[:z])
            expand_extrema!(sp[:xaxis], (0,w))
            expand_extrema!(sp[:yaxis], (0,h))
            sp[:yaxis].d[:flip] = true
        elseif !(st in (:pie, :histogram, :histogram2d))
            expand_extrema!(sp, d)
        end

        # add the series!
        warnOnUnsupported_args(plt.backend, d)
        warnOnUnsupported(plt.backend, d)
        series = Series(d)
        push!(plt.series_list, series)
        # @show series

        _series_added(plt, series)

    else
        # get a sub list of series for this seriestype
        datalist = RecipesBase.apply_recipe(d, Val{st}, d[:x], d[:y], d[:z])
        # datalist = try
        #     RecipesBase.apply_recipe(d, Val{st}, d[:x], d[:y], d[:z])
        # catch
        #     warn("Exception during apply_recipe(Val{$st}, ...) with types ($(typeof(d[:x])), $(typeof(d[:y])), $(typeof(d[:z])))")
        #     rethrow()
        # end

        # assuming there was no error, recursively apply the series recipes
        for data in datalist
            if isa(data, RecipeData)
                _apply_series_recipe(plt, data.d)
            else
                warn("Unhandled recipe: $(data)")
                break
            end
        end
    end
end

function command_idx(kw_list::AVec{KW}, kw::KW)
    kw[:series_plotindex] - kw_list[1][:series_plotindex] + 1
end

# this is the core plotting function.  recursively apply recipes to build
# a list of series KW dicts.
# note: at entry, we only have those preprocessed args which were passed in... no default values yet
function _plot!(plt::Plot, d::KW, args...)
    d[:plot_object] = plt

    # the grouping mechanism is a recipe on a GroupBy object
    # we simply add the GroupBy object to the front of the args list to allow
    # the recipe to be applied
    if haskey(d, :group)
        args = (extractGroupArgs(d[:group], args...), args...)
    end

    # if we were passed a vector/matrix of seriestypes and there's more than one row,
    # we want to duplicate the inputs, once for each seriestype row.
    kw_list = KW[]
    still_to_process = if isempty(args)
        []
    else
        sts = get(d, :seriestype, :path)
        if typeof(sts) <: AbstractArray
            [begin
                dc = copy(d)
                dc[:seriestype] = sts[r,:]
                RecipeData(dc, args)
            end for r=1:size(sts,1)]
        else
            [RecipeData(copy(d), args)]
        end
    end

    # remove subplot and axis args from d... they will be passed through in the kw_list
    if !isempty(args)
        for (k,v) in d
            for defdict in (_subplot_defaults,
                            _axis_defaults,
                            _axis_defaults_byletter)
                if haskey(defdict, k)
                    delete!(d, k)
                end
            end
        end
    end

    # --------------------------------
    # "USER RECIPES"
    # --------------------------------

    # for plotting recipes, swap out the args and update the parameter dictionary
    # we are keeping a queue of series that still need to be processed.
    # each pass through the loop, we pop one off and apply the recipe.
    # the recipe will return a list a Series objects... the ones that are
    # finished (no more args) get added to the kw_list, and the rest go into the queue
    # for processing.
    while !isempty(still_to_process)

        # grab the first in line to be processed and pass it through apply_recipe
        # to generate a list of RecipeData objects (data + attributes)
        next_series = shift!(still_to_process)
        for recipedata in RecipesBase.apply_recipe(next_series.d, next_series.args...)

            # recipedata should be of type RecipeData.  if it's not then the inputs must not have been fully processed by recipes
            if !(typeof(recipedata) <: RecipeData)
                error("Inputs couldn't be processed... expected RecipeData but got: $recipedata")
            end

            if isempty(recipedata.args)
                # when the arg tuple is empty, that means there's nothing left to recursively
                # process... finish up and add to the kw_list
                kw = recipedata.d
                _add_markershape(kw)

                # if there was a grouping, filter the data here
                _filter_input_data!(kw)

                # map marker_z if it's a Function
                if isa(get(kw, :marker_z, nothing), Function)
                    # TODO: should this take y and/or z as arguments?
                    kw[:marker_z] = map(kw[:marker_z], kw[:x], kw[:y], kw[:z])
                end

                # map line_z if it's a Function
                if isa(get(kw, :line_z, nothing), Function)
                    kw[:line_z] = map(kw[:line_z], kw[:x], kw[:y], kw[:z])
                end

                # convert a ribbon into a fillrange
                if get(kw, :ribbon, nothing) != nothing
                    make_fillrange_from_ribbon(kw)
                end

                # add the plot index
                plt.n += 1
                kw[:series_plotindex] = plt.n

                # check that the backend will support the command and add it to the list
                warnOnUnsupported_scales(plt.backend, kw)
                push!(kw_list, kw)

                # handle error bars by creating new recipedata data... these will have
                # the same recipedata index as the recipedata they are copied from
                for esym in (:xerror, :yerror)
                    if get(d, esym, nothing) != nothing
                        # we make a copy of the KW and apply an errorbar recipe
                        errkw = copy(kw)
                        errkw[:seriestype] = esym
                        errkw[:label] = ""
                        errkw[:primary] = false
                        push!(kw_list, errkw)
                    end
                end

                # handle smoothing by adding a new series
                if get(d, :smooth, false)
                    x, y = kw[:x], kw[:y]
                    β, α = convert(Matrix{Float64}, [x ones(length(x))]) \ convert(Vector{Float64}, y)
                    sx = [minimum(x), maximum(x)]
                    sy = β * sx + α
                    push!(kw_list, merge(copy(kw), KW(
                        :seriestype => :path,
                        :x => sx,
                        :y => sy,
                        :fillrange => nothing,
                        :label => "",
                        :primary => false,
                    )))
                end

            else
                # args are non-empty, so there's still processing to do... add it back to the queue
                push!(still_to_process, recipedata)
            end
        end
    end

    # don't allow something else to handle it
    d[:smooth] = false

    # --------------------------------
    # "PLOT RECIPES"
    # --------------------------------

    # "plot recipe", which acts like a series type, and is processed before
    # the plot layout is created, which allows for setting layouts and other plot-wide attributes.
    # we get inputs which have been fully processed by "user recipes" and "type recipes",
    # so we can expect standard vectors, surfaces, etc.  No defaults have been set yet.
    still_to_process = kw_list
    kw_list = KW[]
    while !isempty(still_to_process)
        # Grab the first in line to be processed and pass it through apply_recipe
        # to generate a list of RecipeData objects (data + attributes).
        # If we applied a "plot recipe" without error, then add the returned datalist's KWs,
        # otherwise we just add the original KW.
        next_kw = shift!(still_to_process)
        if !isa(get(next_kw, :seriestype, nothing), Symbol)
            # seriestype was never set, or it's not a Symbol, so it can't be a plot recipe
            push!(kw_list, next_kw)
            continue
        end
        try
            st = next_kw[:seriestype]
            st = next_kw[:seriestype] = get(_typeAliases, st, st)
            datalist = RecipesBase.apply_recipe(next_kw, Val{st}, plt)
            for data in datalist
                if data.d[:seriestype] == st
                    error("Plot recipe $st returned the same seriestype: $(data.d)")
                end
                push!(still_to_process, data.d)
            end
        catch err
            if isa(err, MethodError)
                push!(kw_list, next_kw)
            else
                rethrow()
            end
        end
    end

    # --------------------------------
    # Plot/Subplot/Layout setup
    # --------------------------------

    # merge in anything meant for the Plot
    for kw in kw_list, (k,v) in kw
        haskey(_plot_defaults, k) && (d[k] = pop!(kw, k))
    end

    # TODO: init subplots here
    _update_plot_args(plt, d)
    if !plt.init
        plt.o = _create_backend_figure(plt)

        # create the layout and subplots from the inputs
        plt.layout, plt.subplots, plt.spmap = build_layout(plt.attr)
        for (idx,sp) in enumerate(plt.subplots)
            sp.plt = plt
            sp.attr[:subplot_index] = idx
        end

        plt.init = true
    end

    # handle inset subplots
    insets = plt[:inset_subplots]
    if insets != nothing
        if !(typeof(insets) <: AVec)
            insets = [insets]
        end
        for inset in insets
            parent, bb = is_2tuple(inset) ? inset : (nothing, inset)
            P = typeof(parent)
            if P <: Integer
                parent = plt.subplots[parent]
            elseif P == Symbol
                parent = plt.spmap[parent]
            else
                parent = plt.layout
            end
            sp = Subplot(backend(), parent=parent)
            sp.plt = plt
            sp.attr[:relative_bbox] = bb
            sp.attr[:subplot_index] = length(plt.subplots)
            push!(plt.subplots, sp)
            push!(plt.inset_subplots, sp)
        end
    end

    # we'll keep a map of subplot to an attribute override dict.
    # Subplot/Axis attributes set by a user/series recipe apply only to the
    # Subplot object which they belong to.
    # TODO: allow matrices to still apply to all subplots
    sp_attrs = Dict{Subplot,Any}()
    for kw in kw_list
        # get the Subplot object to which the series belongs.
        sps = get(kw, :subplot, :auto)
        sp = get_subplot(plt, cycle(sps == :auto ? plt.subplots : plt.subplots[sps], command_idx(kw_list,kw)))
        kw[:subplot] = sp

        attr = KW()
        for (k,v) in kw
            for defdict in (_subplot_defaults,
                            _axis_defaults,
                            _axis_defaults_byletter)
                if haskey(defdict, k)
                    attr[k] = pop!(kw, k)
                end
            end
        end
        sp_attrs[sp] = attr
    end

    # override subplot/axis args.  `sp_attrs` take precendence
    for (idx,sp) in enumerate(plt.subplots)
        attr = merge(d, get(sp_attrs, sp, KW()))
        _update_subplot_args(plt, sp, attr, idx, remove_pair = false)
    end

    # do we need to link any axes together?
    link_axes!(plt.layout, plt[:link])

    # !!! note: At this point, kw_list is fully decomposed into individual series... one KW per series.          !!!
    # !!!       The next step is to recursively apply series recipes until the backend supports that series type !!!

    # --------------------------------
    # "SERIES RECIPES"
    # --------------------------------

    for kw in kw_list
        sp = kw[:subplot]
        idx = get_subplot_index(plt, sp)

        # # we update subplot args in case something like the color palatte is part of the recipe
        # _update_subplot_args(plt, sp, kw, idx)

        # set default values, select from attribute cycles, and generally set the final attributes
        _add_defaults!(kw, plt, sp, command_idx(kw_list,kw))

        # now we have a fully specified series, with colors chosen.   we must recursively handle
        # series recipes, which dispatch on seriestype.  If a backend does not natively support a seriestype,
        # we check for a recipe that will convert that series type into one made up of lower-level components.
        # For example, a histogram is just a bar plot with binned data, a bar plot is really a filled step plot,
        # and a step plot is really just a path.  So any backend that supports drawing a path will implicitly
        # be able to support step, bar, and histogram plots (and any recipes that use those components).
        _apply_series_recipe(plt, kw)
    end

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

    current(plt)

    # note: lets ignore the show param and effectively use the semicolon at the end of the REPL statement
    # # do we want to show it?
    # if haskey(d, :show) && d[:show]
    if get(d, :show, default(:show))
        gui()
    end

    plt
end

function _replace_linewidth(d::KW)
    # get a good default linewidth... 0 for surface and heatmaps
    if get(d, :linewidth, :auto) == :auto
        d[:linewidth] = (get(d, :seriestype, :path) in (:surface,:heatmap,:image) ? 0 : 1)
    end
end

# we're getting ready to display/output.  prep for layout calcs, then update
# the plot object after
function prepare_output(plt::Plot)
    _before_layout_calcs(plt)

    w, h = plt.attr[:size]
    plt.layout.bbox = BoundingBox(0mm, 0mm, w*px, h*px)

    # One pass down and back up the tree to compute the minimum padding
    # of the children on the perimeter.  This is an backend callback.
    _update_min_padding!(plt.layout)
    for sp in plt.inset_subplots
        _update_min_padding!(sp)
    end

    # now another pass down, to update the bounding boxes
    update_child_bboxes!(plt.layout)

    # update those bounding boxes of inset subplots
    update_inset_bboxes!(plt)

    # the backend callback, to reposition subplots, etc
    _update_plot_object(plt)
end

function prepared_object(plt::Plot)
    prepare_output(plt)
    plt.o
end

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

# function get_indices(orig, labels)
#     Int[findnext(labels, l, 1) for l in orig]
# end

# # TODO: remove?? this is the old way of handling discrete data... should be
# # replaced by the Axis type and logic
# function setTicksFromStringVector(plt::Plot, d::KW, di::KW, letter)
#     sym = Symbol(letter)
#     ticksym = Symbol(letter * "ticks")
#     pargs = plt.attr
#     v = di[sym]
#
#     # do we really want to do this?
#     typeof(v) <: AbstractArray || return
#     isempty(v) && return
#     trueOrAllTrue(_ -> typeof(_) <: AbstractString, v) || return
#
#     # compute the ticks and labels
#     ticks, labels = if ticksType(pargs[ticksym]) == :ticks_and_labels
#         # extend the existing ticks and labels. only add to labels if they're new!
#         ticks, labels = pargs[ticksym]
#         newlabels = filter(_ -> !(_ in labels), unique(v))
#         newticks = if isempty(ticks)
#             collect(1:length(newlabels))
#         else
#             maximum(ticks) + collect(1:length(newlabels))
#         end
#         ticks = vcat(ticks, newticks)
#         labels = vcat(labels, newlabels)
#         ticks, labels
#     else
#         # create new ticks and labels
#         newlabels = unique(v)
#         collect(1:length(newlabels)), newlabels
#     end
#
#     d[ticksym] = ticks, labels
#     plt.attr[ticksym] = ticks, labels
#
#     # add an origsym field so that later on we can re-compute the x vector if ticks change
#     origsym = Symbol(letter * "orig")
#     di[origsym] = v
#     di[sym] = get_indices(v, labels)
#
#     # loop through existing plt.seriesargs and adjust indices if there is an origsym key
#     for sargs in plt.seriesargs
#         if haskey(sargs, origsym)
#             # TODO: might need to call the setindex function instead to trigger a plot update for some backends??
#             sargs[sym] = get_indices(sargs[origsym], labels)
#         end
#     end
# end

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

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

# function Base.copy(plt::Plot)
#     backend(plt.backend)
#     plt2 = plot(; plt.attr...)
#     for sargs in plt.seriesargs
#         sargs = filter((k,v) -> haskey(_series_defaults,k), sargs)
#         plot!(plt2; sargs...)
#     end
#     plt2
# end

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


text:
#=  Here I compare a wide variety of tiebreaking rules, some of which (STB, MTB, and
    HTB) are described by Ashlagi and Nikzad (2020), and others (XTB, WTB) of my own
    creation. Testing the tiebreaking rules over a hybrid market in which the students
    uniformly prefer a subset of the schools allows us to examine the properties of
    these mechanisms in over- and underdemanded markets at a glance. See the readme
    for further discussion.                     =#

using Permutations
using Plots
using DeferredAcceptance

"""
Plots outcomes for various tiebreaking rules for a market with
the given characteristics.
"""
function plotter_more(n, m_pop, m_unp, samp)

    m = m_pop + m_unp

    descr = "$n students, $m schools ($m_pop popular, $m_unp unpopular), $samp samples"
    println(descr)

    # All schools have unit capacities
    capacities = ones(Int64, m)

    # Blank rank dists
    cdf_STB = zeros(Float64, m)
    cdf_MTB = zeros(Float64, m)
    cdf_HTB = zeros(Float64, m)
    cdf_XTB = zeros(Float64, m)
    cdf_WHTB = zeros(Float64, m)
    cdf_WXTB = zeros(Float64, m)

    # For HTB and WTB
    blend = ones(Float64, 1, m) # Use MTB in all schools except
    blend[1:m_pop] .= 0            # the popular schools, which use STB

    for i in 1:samp
        # Students unilaterally prefer popular schools to unpopular ones
        # Random ranking otherwise
        students = vcat(hcat((randperm(m_pop) for i in 1:n)...),
                   m_pop .+ hcat((randperm(m_unp) for i in 1:n)...))

        # Schools place all students in single priority category
        schools = ones(Int64, n, m)

        # Break ties
        schools_STB = singletiebreaking(schools)
        schools_MTB = multipletiebreaking(schools)
        schools_HTB = hybridtiebreaking(schools, blend)
        schools_XTB = hybridtiebreaking(schools, 0.5)
        schools_WHTB = welfaretiebreaking(students, schools, blend)
        schools_WXTB = welfaretiebreaking(students, schools, 0.5)

        # Update rank risks
        # println("Starting STB $i")
        cdf_STB += DA_rank_dist(students, schools_STB, capacities)
        # println("Starting MTB $i")
        cdf_MTB += DA_rank_dist(students, schools_MTB, capacities)
        # println("Starting HTB $i")
        cdf_HTB += DA_rank_dist(students, schools_HTB, capacities)
        # println("Starting XTB $i")
        cdf_XTB += DA_rank_dist(students, schools_XTB, capacities)
        # println("Starting WHTB $i")
        cdf_WHTB += DA_rank_dist(students, schools_WHTB, capacities)
        # println("Starting WXTB $i")
        cdf_WXTB += DA_rank_dist(students, schools_WXTB, capacities)

        # To display plot as it updates
        # display(plot([cdf_STB, cdf_MTB, cdf_HTB, cdf_XTB, cdf_WHTB, cdf_WXTB],
        #         label = ["DA-STB" "DA-MTB" "DA-HTB" "DA-XTB, λ=0.5" "DA-WHTB" "DA-WXTB, λ=0.5"],
        #         lc = [:dimgray :dimgray :dimgray :crimson :dodgerblue :olivedrab],
        #         ls = [:dot :dash :dashdot :dot :dash :dashdot],
        #         legend = :bottomright,
        #         title = descr, titlefontsize=11,
        #         xlabel = "rank", ylabel= "average number of students"))

    end

    # Norm rank dists against sample size
    for i in (cdf_STB, cdf_MTB, cdf_HTB, cdf_XTB, cdf_WHTB, cdf_WXTB)
        i ./= samp
    end

    return plot([cdf_STB, cdf_MTB, cdf_HTB, cdf_XTB, cdf_WHTB, cdf_WXTB],
                label = ["DA-STB" "DA-MTB" "DA-HTB" "DA-XTB, λ=0.5" "DA-WHTB" "DA-WXTB, λ=0.5"],
                lc = [:dimgray :dimgray :dimgray :crimson :dodgerblue :olivedrab],
                ls = [:dot :dash :dashdot :dot :dash :dashdot],
                legend = :bottomright,
                title = descr, titlefontsize=11,
                xlabel = "rank", ylabel= "average number of students")
end

(n, m_pop, m_unp, samp) = (50, 17, 33, 200)

p = plotter_more(n, m_pop, m_unp, samp)

# savefig(p, string("examples/hybrid/plot.pdf"))
# savefig(p, string("examples/hybrid/plot.png"))

max_stars_repo_name: maxkapur/DeferredAcceptance

max_issues_repo_name: maxkapur/DeferredAcceptance

max_forks_repo_name: maxkapur/DeferredAcceptance


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "MMG"
version = v"5.6.0"

# Collection of sources required to build MMG
sources = [
    ArchiveSource("https://github.com/MmgTools/mmg/archive/refs/tags/v$(version).tar.gz",
                  "bbf9163d65bc6e0f81dd3acc5a51e4a8c47a7fdae849abc26277e01154fe2437"),
    DirectorySource("./bundled")
]

# Bash recipe for building across all platforms
script = raw"""
# Install genheader for host platform
cp -r ${WORKSPACE}/srcdir/mmg-* ${WORKSPACE}/srcdir/mmg-genheader
cd ${WORKSPACE}/srcdir/mmg-genheader
atomic_patch -p1 "${WORKSPACE}/srcdir/patches/genheader.patch"
mkdir build
cd build
cmake .. \
    -DCMAKE_INSTALL_PREFIX=${host_prefix} \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_HOST_TOOLCHAIN} \
    -DCMAKE_BUILD_TYPE=Release
make -j${nproc}
make install
cd ${WORKSPACE}/srcdir && rm -r ${WORKSPACE}/srcdir/mmg-genheader

# Install MMG
cd ${WORKSPACE}/srcdir/mmg-*
if [[ "${target}" == *mingw* ]]; then
    atomic_patch -p1 "${WORKSPACE}/srcdir/patches/MMG.mingw.patch"
    USE_SCOTCH=OFF
else
    atomic_patch -p1 "${WORKSPACE}/srcdir/patches/MMG.patch"
    USE_SCOTCH=ON
fi
mkdir build
cd build
cmake .. \
    -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_BUILD_TYPE=Release \
    -DBUILD_SHARED_LIBS=ON \
    -DUSE_SCOTCH=${USE_SCOTCH} \
    -DUSE_ELAS=ON \
    -DUSE_VTK=OFF
make -j${nproc}
make install
install_license ../LICENSE
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct("libmmg", :libmmg),
    LibraryProduct("libmmg2d", :libmmg2d),
    LibraryProduct("libmmg3d", :libmmg3d),
    LibraryProduct("libmmgs", :libmmgs),
    ExecutableProduct("mmg2d_O3", :mmg2d_O3),
    ExecutableProduct("mmg3d_O3", :mmg3d_O3),
    ExecutableProduct("mmgs_O3", :mmgs_O3)
]

# SCOTCH is only available on non-Windows platforms
scotch_platforms = filter(!Sys.iswindows, platforms)

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("LinearElasticity_jll"),
    Dependency("SCOTCH_jll", platforms=scotch_platforms)
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6", preferred_gcc_version=v"9")


text:
function derrivatives(state::Array,GM)
    r = sqrt(state[1]^2+state[2]^2)
    dvx = -GM*state[1]/r^3
    dvy = -GM*state[2]/r^3
    drx = state[3]
    dry = state[4]
  return [drx,dry,dvx,dvy]
end

function leapfrog_update_pos(state::Array,dt)
  Rx = state[1] + .5*dt*state[3]
  Ry = state[2] + .5*dt*state[4]
  newstate = [Rx,Ry,state[3],state[4]]
  return newstate
end

function leapfrog_update_both(state::Array,ds,dt)
  vx = state[3]+dt*ds[3]
  vy = state[4]+dt*ds[4]
  [state[1]+.5*dt*vx,state[2]+.5*dt*vy,vx,vy]
end

function integrate_leapfrog_student(state::Array,dt,duration,GM=1)
  r_x,r_y,v_x,v_y = [],[],[],[]
  N = iceil(duration/dt)
  t = 0.
  for i in 0:N
    dt_tmp = (t+dt<=duration) ? dt : duration-t;

    state = leapfrog_update_pos(state,dt_tmp)
    ds = derrivatives(state,GM)
    state = leapfrog_update_both(state,ds,dt_tmp)
    t += dt_tmp

    r_x = vcat(r_x,[state[1]])
    r_y = vcat(r_y,[state[2]])
    v_x = vcat(v_x,[state[3]])
    v_y = vcat(v_y,[state[4]])
  end
  return (r_x,r_y,v_x,v_y)
end

function calc_end_distance_leapfrog_student(dur::Integer, dt::Float64 = 2pi/200.0, state::Vector{Float64} = [1., 0., 0., 1.] )
  #state = [1.,0.,0.,1.];
  dist_init = sqrt(state[1]^2+state[2]^2)
  phase_init = atan2(state[2],state[1])
  datalog = integrate_leapfrog_student(state,dt,dur*2pi)
  # Since student version didn't modify state (no !), we need to update state with the final value
  state = [datalog[1][end], datalog[2][end], datalog[3][end], datalog[4][end]]
  # Calculate three metrics of the accuracy of the integration
  dist = sqrt(state[1]^2+state[2]^2)
  phase = atan2(state[2],state[1])
  offset = sqrt(sum((state[1:2].-[1.0,0.0]).^2))
  return (dist-dist_init,phase-phase_init,offset)
end

using Base.Test
  function test_leapfrog_student(dur::Integer, dt::Float64 = 2pi/200.0)
  err = calc_end_distance_leapfrog_student(dur,dt)
  @test_approx_eq_eps(err[1], 0., 1e-6)
  @test_approx_eq_eps(err[2], 0., 1e-1)
  @test_approx_eq_eps(err[3], 0., 1e-1)
end


text:
#===============================================================================

    AmpModel.jl

    heirarchy of abstract types
    describing available amplification models

    Author: Tom Price
    Date:   June 2019

===============================================================================#

abstract type AmpModel                      end
abstract type DFCModel          <: AmpModel end
abstract type MAK2              <: DFCModel end
abstract type MAK3              <: DFCModel end
abstract type MAKERGAUL3        <: DFCModel end
abstract type MAKERGAUL4        <: DFCModel end
abstract type SFCModel          <: AmpModel end

const AMPMODELS = ## baseline models
    [   subtypes(AmpModel)...,
        subtypes(DFCModel)...   ]

max_stars_repo_name: MakerButt/chaipcr

max_issues_repo_name: MakerButt/chaipcr

max_forks_repo_name: MakerButt/chaipcr


text:
#!/bin/julia

import Printf.@printf
using Combinatorics

const SUITS = "SDCH"
const RANKS = "A23456789TJQK"

@enum Suit begin
    Spades = 1
    Diamonds
    Clubs
    Hearts
end

@enum Rank begin
    Ace = 1
    Deuce
    Three
    Four
    Five
    Six
    Seven
    Eight
    Nine
    Ten
    Jack
    Queen
    King
end

struct Card
    rank::Rank
    suit::Suit
end

function Base.show(io::IO, card::Card)
    print(io, RANKS[Int(card.rank)])
    print(io, SUITS[Int(card.suit)])
end

function order(card::Card)::Int
    return Int(card.rank)
end

Cards = Array{Card, 1}

struct Hand
    cards::Cards
    Hand() = new(Cards())
end

function push(hand::Hand, card::Card)
    @assert !(card in hand.cards)
    push!(hand.cards, card)
end

function pop(hand::Hand)
    pop!(hand.cards)
end

function length(hand::Hand)::Int
    return Base.length(hand.cards)
end

function value(hand::Hand, index::Int)::Int
    return min(10, Int(hand.cards[index].rank))
end

function Base.show(io::IO, hand::Hand)
    sep = false
    for card in hand.cards
        if sep
            print(io, ' ')
        end
        print(io, card)
        sep = true
    end
end

function make_hand(s::String)::Hand
    hand = Hand()
    rank = nothing

    for c in uppercase(s)
        i = findfirst(c, SUITS)
        if !isnothing(i)
            if isnothing(rank)
                error("Malformed hand: $s")
            end
            push(hand, Card(rank, Suit(i)))
            rank = nothing
            continue
        end

        i = findfirst(c, RANKS)
        if !isnothing(i)
            if !isnothing(rank)
                error("Malformed hand: $s")
            end
            rank = Rank(i)
            continue
        end

        if c != ' ' && c != '-'
            error("Malformed hand: $s")
        end
    end

    if !isnothing(rank)
        error("Malformed hand: $s")
    end

    return hand
end

@assert "5H 5C 5S JD 5D" == string(make_hand("5H 5C 5S JD 5D"))
@assert "5H 5C 5S JD 5D" == string(make_hand("5h5c5sjd5d"))
@assert "AH AS JH AC AD" == string(make_hand("ah-as-jh-ac-ad"))

function score_15s(hand::Hand)::Int
    @assert length(hand) == 5

    a = value(hand, 1)
    b = value(hand, 2)
    c = value(hand, 3)
    d = value(hand, 4)
    e = value(hand, 5)

    num_15s = 0

    # five cards - C(5,5)=1
    if a + b + c + d + e == 15
        num_15s += 1
    end

    # four cards - C(5,4)=5
    if a + b + c + d == 15
        num_15s += 1
    end
    if a + b + c + e == 15
        num_15s += 1
    end
    if a + b + d + e == 15
        num_15s += 1
    end
    if a + c + d + e == 15
        num_15s += 1
    end
    if b + c + d + e == 15
        num_15s += 1
    end

    # three cards - C(5,3)=10
    if a + b + c == 15
        num_15s += 1
    end
    if a + b + d == 15
        num_15s += 1
    end
    if a + b + e == 15
        num_15s += 1
    end
    if a + c + d == 15
        num_15s += 1
    end
    if a + c + e == 15
        num_15s += 1
    end
    if a + d + e == 15
        num_15s += 1
    end
    if b + c + d == 15
        num_15s += 1
    end
    if b + c + e == 15
        num_15s += 1
    end
    if b + d + e == 15
        num_15s += 1
    end
    if c + d + e == 15
        num_15s += 1
    end

    # two cards - C(5,2)=10
    if a + b == 15
        num_15s += 1
    end
    if a + c == 15
        num_15s += 1
    end
    if a + d == 15
        num_15s += 1
    end
    if a + e == 15
        num_15s += 1
    end
    if b + c == 15
        num_15s += 1
    end
    if b + d == 15
        num_15s += 1
    end
    if b + e == 15
        num_15s += 1
    end
    if c + d == 15
        num_15s += 1
    end
    if c + e == 15
        num_15s += 1
    end
    if d + e == 15
        num_15s += 1
    end

    return 2 * num_15s
end

@assert  4 == score_15s(make_hand("AH 2H 3H JH QH"))
@assert  8 == score_15s(make_hand("5H 2H 3H JH QH"))
@assert 16 == score_15s(make_hand("5H 5S 5C 5D TH"))
@assert  8 == score_15s(make_hand("6C 6D 4D 4S 5D"))

function score_pairs(hand::Hand)::Int
    num_pairs = 0
    for (ai, a) in enumerate(hand.cards)
        for b in Iterators.drop(hand.cards, ai)
            if a.rank == b.rank
                num_pairs += 1
            end
        end
    end
    return 2 * num_pairs
end

@assert 12 == score_pairs(make_hand("5H 5S 5C 5D TH"))
@assert  8 == score_pairs(make_hand("TS 5S 5C 5D TH"))
@assert  4 == score_pairs(make_hand("6C 6D 4D 4S 5D"))

const X = -1            # match any rank
const RUN_PATTERNS = (
    # (score, (deltas...)),
    (12, (0, 1, 1, 0)), # AA233
    ( 9, (1, 1, 0, 0)), # A2333
    ( 9, (1, 0, 0, 1)), # A2223
    ( 9, (0, 0, 1, 1)), # AAA23
    ( 8, (1, 1, 1, 0)), # A2344
    ( 8, (1, 1, 0, 1)), # A2334
    ( 8, (1, 0, 1, 1)), # A2234
    ( 8, (0, 1, 1, 1)), # AA234
    ( 6, (X, 1, 1, 0)), # xA233
    ( 6, (X, 1, 0, 1)), # xA223
    ( 6, (X, 0, 1, 1)), # xAA23
    ( 6, (1, 1, 0, X)), # A233x
    ( 6, (1, 0, 1, X)), # A223x
    ( 6, (0, 1, 1, X)), # AA23x
    ( 5, (1, 1, 1, 1)), # A2345
    ( 4, (X, 1, 1, 1)), # xA234
    ( 4, (1, 1, 1, X)), # A234x
    ( 3, (X, X, 1, 1)), # xxA23
    ( 3, (X, 1, 1, X)), # xA23x
    ( 3, (1, 1, X, X)), # A23xx
)

function score_runs(hand::Hand)::Int
    @assert length(hand) == 5

    # Make a sorted sequence of the orders of the cards in the hand.  The
    # order of Ace is 1, ..., King is 13
    orders = map(card -> order(card), hand.cards)
    sort!(orders)

    # Compare the sorted hand to the patterns.  Look at the difference between
    # the two cards in each pair of adjacent cards.  Stop at the first match
    for (score, deltas) in RUN_PATTERNS
        previous = orders[1]
        j = 1
        while true
            delta = deltas[j]
            order = orders[j + 1]
            if delta != X && delta != order - previous
                break
            end
            previous = order
            j += 1
            if j == 5
                return score
            end
        end
    end
    return 0
end

@assert  9 == score_runs(make_hand("AH 2H 3H 3D 3C"))
@assert  9 == score_runs(make_hand("KH KD KC JH QH"))  # same pattern A2333
@assert  9 == score_runs(make_hand("AH 2H 2D 2C 3H"))
@assert  9 == score_runs(make_hand("AH AD AC 2H 3H"))
@assert  8 == score_runs(make_hand("AH 2H 3H 4H 4D"))
@assert  8 == score_runs(make_hand("AH 2H 3H 3D 4H"))
@assert  8 == score_runs(make_hand("AH 2H 2C 3H 4H"))
@assert  8 == score_runs(make_hand("AS AH 2H 3H 4H"))
@assert  6 == score_runs(make_hand("JH AH 2H 3D 3H"))
@assert  6 == score_runs(make_hand("JH AH 2S 2H 3H"))
@assert  6 == score_runs(make_hand("JH AH AS 2H 3H"))
@assert  6 == score_runs(make_hand("AH 2H 3S 3H JH"))
@assert  6 == score_runs(make_hand("AH 2H 2S 3H JH"))
@assert  6 == score_runs(make_hand("AH AS 2H 3H JH"))
@assert  5 == score_runs(make_hand("AH 2H 3H 4H 5H"))
@assert  4 == score_runs(make_hand("JH AH 2H 3H 4H"))
@assert  4 == score_runs(make_hand("AH 2H 3H 4H JH"))
@assert  3 == score_runs(make_hand("JH QH AH 2H 3H"))
@assert  3 == score_runs(make_hand("JH AH 2H 3H TH"))
@assert  3 == score_runs(make_hand("AH 2H 3H JH TH"))
@assert  0 == score_runs(make_hand("AH 8H 3H JH TH"))
@assert 12 == score_runs(make_hand("6C 6D 4D 4S 5D"))

function score_flush(hand::Hand, is_crib::Bool)::Int
    @assert length(hand) == 5
    suit = hand.cards[1].suit
    for card in Iterators.take(Iterators.drop(hand.cards, 1), 3)
        if suit != card.suit
            return 0
        end
    end
    # First 4 are the same suit, check the cut card
    if suit == hand.cards[5].suit
        return 5
    end
    # In the crib, a flush counts only if all five cards are the same suit
    if is_crib
        return 0
    end
    return 4
end

@assert 5 == score_flush(make_hand("5H 6H 7H 8H 9H"), false)
@assert 4 == score_flush(make_hand("5H 6H 7H 8H 9D"), false)
@assert 0 == score_flush(make_hand("5H 6H 7H 8H 9D"), true)
@assert 0 == score_flush(make_hand("5H 6H 7H 8D 9D"), false)

function score_nobs(hand::Hand)::Int
    # nobs: one point for the Jack of the same suit as the cut card
    @assert length(hand) == 5
    cut_suit = hand.cards[5].suit
    for card in Iterators.take(hand.cards, 4)
        if card.rank == Jack && card.suit == cut_suit
            return 1
        end
    end
    return 0
end

@assert 1 == score_nobs(make_hand("JH 2C 3C 4C 5H"))
@assert 0 == score_nobs(make_hand("JH 2C 3C 4C 5C"))

function score_hand(hand::Hand, is_crib::Bool)::Int
    return score_15s(hand) +
           score_pairs(hand) +
           score_runs(hand) +
           score_flush(hand, is_crib) +
           score_nobs(hand)
end

function score_hand(hand::String, is_crib::Bool)::Int
    return score_hand(make_hand(hand), is_crib)
end

@assert 12 == score_hand("AH AS JH AC AD", false) # 4oak ("of a kind")
@assert 13 == score_hand("AH AS JD AC AD", false) # ...plus right jack
@assert  5 == score_hand("AH 3H 7H TH JH", false) # 5 hearts
@assert  5 == score_hand("AH 3H 7H TH JH", true)  # 5 hearts but crib
@assert  4 == score_hand("AH 3H 7H TH JS", false) # 4 hearts
@assert  0 == score_hand("AH 3H 7S TH JH", false) # 4 hearts but with cut
@assert  0 == score_hand("AH 3H 7H TH JS", true)  # 4 hearts but crib
@assert  7 == score_hand("AH 2S 3C 5D JH", false) # 15/4 + run/3
@assert 20 == score_hand("7H 7S 7C 8D 8H", false) # 15/12 + 3oak + 2oak
@assert 15 == score_hand("AH 2H 3H 3S 3D", false) # triple run/3
@assert 15 == score_hand("3H AH 3S 2H 3D", false) # triple run/3
@assert 29 == score_hand("5H 5C 5S JD 5D", false)
@assert 28 == score_hand("5H 5C 5S 5D JD", false)
@assert 24 == score_hand("6C 4D 6D 4S 5D", false)

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

function make_deck(exclude::Hand)::Hand
    # Make an entire deck of cards but leave out any cards in `exclude`
    deck = Hand()
    for suit in instances(Suit)
        for rank in instances(Rank)
            card = Card(rank, suit)
            if !(card in exclude.cards)
                push(deck, card)
            end
        end
    end
    return deck
end

make_deck() = make_deck(Hand())

const MAX_SCORE = 29 + 24  # 29 in hand, 24 in crib (44665)
const MIN_SCORE = -29      # 0 in hand, 29 in opp crib
const NUM_SCORES = MAX_SCORE - MIN_SCORE + 1

struct Tally
    scores::Array{Int32, 1}

    Tally() = new(zeros(Int32, NUM_SCORES))
end

function count(tally::Tally, score::Int64)
    tally.scores[score - MIN_SCORE + 1] += 1
end

struct Statistics
    mean::Float32
    stdev::Float32
    min_score::Int
    max_score::Int
end

function Base.show(io::IO, stats::Statistics)
    @printf(io, "%.1f %.1f %d..%d", stats.mean, stats.stdev, stats.min_score, stats.max_score)
end

function make_statistics(tally::Tally, num_hands::Int64)::Statistics
    # Convert `tally` to a Statistics object.  `tally` is the number of times
    # each score MIN_SCORE..MAX_SCORE was achieved over `num_hands` hands.
    mins = 0
    for score in MIN_SCORE:MAX_SCORE
        if tally.scores[score - MIN_SCORE + 1] != 0
            mins = score
            break
        end
    end

    maxs = 0
    for score in MAX_SCORE:-1:MIN_SCORE
        if tally.scores[score - MIN_SCORE + 1] != 0
            maxs = score
            break
        end
    end

    sums = 0.0
    for score in mins:maxs
        sums += score * tally.scores[score - MIN_SCORE + 1]
    end
    mean = sums / num_hands

    sumdev = 0.0
    for score in mins:maxs
        d = score - mean
        sumdev += d * d
    end
    stdev = sqrt(sumdev / num_hands)

    return Statistics(mean, stdev, mins, maxs)
end

function analyze_hand(hand::Hand)
    # Find all possible pairs of cards to discard to the crib.
    # There are C(6,2)=15 possible discards in a cribbage hand.

    for discarding in combinations(hand.cards, 2)
        hold = Hand()
        for card in hand.cards
            if !(card in discarding)
                push(hold, card)
            end
        end

        deck = make_deck(hand)
        @assert length(deck) == 46

        mine_tally = Tally()   # scores when the crib is mine
        theirs_tally = Tally() # scores then the crib is theirs
        num_hands = 0
        for chosen in combinations(deck.cards, 2)
            card1 = chosen[1]
            card2 = chosen[2]

            crib = Hand()
            push(crib, discarding[1])
            push(crib, discarding[2])
            push(crib, card1)
            push(crib, card2)
            @assert length(crib) == 4

            for cut in deck.cards
                if cut == card1 || cut == card2
                    continue
                end

                push(hold, cut)
                hold_score = score_hand(hold, false)
                pop(hold)

                push(crib, cut)
                crib_score = score_hand(crib, true)
                pop(crib)

                mine_score = hold_score + crib_score
                theirs_score = hold_score - crib_score

                num_hands += 1

                count(mine_tally, mine_score)
                count(theirs_tally, theirs_score)
            end
        end

        # deck size: 46, C(46,2)=1035
        # remaining_deck size: 44
        @assert num_hands == 1035 * 44

        if_mine = make_statistics(mine_tally, num_hands)
        if_theirs = make_statistics(theirs_tally, num_hands)

        println("$(discarding[1]) $(discarding[2]) [$if_mine] [$if_theirs]")
    end
end

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

function main()
    for arg in ARGS
        hand = make_hand(arg)
        if length(hand) != 6
            error("Wrong number of cards in hand: $hand")
        end

        println("[ $hand ]")
        analyze_hand(hand)
        println()
    end
end

main()

max_stars_repo_name: MichaelCook/cribbage

max_issues_repo_name: MichaelCook/cribbage

max_forks_repo_name: MichaelCook/cribbage


text:
"""
    Gradient(x, y, z)

Create a `Gradient` object representing x, y, and z B0 gradients. Units are
G/cm.

# Properties
- `x::Real`: x gradient
- `y::Real`: y gradient
- `z::Real`: z gradient
"""
struct Gradient{T<:Real}
    x::T
    y::T
    z::T
end

Gradient(x, y, z) = Gradient(promote(x, y, z)...)

Base.show(io::IO, grad::Gradient) = print(io, "Gradient(", grad.x, ", ", grad.y, ", ", grad.z, ")")
Base.show(io::IO, ::MIME"text/plain", grad::Gradient{T}) where {T} =
    print(io, "Gradient{$T}:\n Gx = ", grad.x, " G/cm\n Gy = ", grad.y, " G/cm\n Gz = ", grad.z, " G/cm")

"""
    gradient_frequency(grad, pos)

Compute the off-resonance frequency in Hz induced by the given B0 gradient
`grad` at position `pos`.
"""
function gradient_frequency(grad::Gradient, pos::Position)

    return GAMBAR * (grad.x * pos.x + grad.y * pos.y + grad.z * pos.z) # Hz

end

"""
    AbstractSpoiling

Abstract type for representing spoiling.
"""
abstract type AbstractSpoiling end

"""
    IdealSpoiling() <: AbstractSpoiling

Represents ideal spoiling, i.e., setting the transverse (x and y) components of
a spin's magnetization to 0.
"""
struct IdealSpoiling <: AbstractSpoiling end

"""
    GradientSpoiling(grad, Tg) <: AbstractSpoiling
    GradientSpoiling(gx, gy, gz, Tg)

Represents gradient spoiling, i.e., applying a gradient
`grad = Gradient(gx, gy, gz)` for time `Tg` (ms).
"""
struct GradientSpoiling{T<:Real} <: AbstractSpoiling
    gradient::Gradient{T}
    Tg::Float64

    function GradientSpoiling(gradient::Gradient{T}, Tg::Real) where {T}

        new{T}(gradient, Tg)

    end
end

GradientSpoiling(gx, gy, gz, Tg) = GradientSpoiling(Gradient(gx, gy, gz), Tg)

Base.show(io::IO, s::GradientSpoiling) = print(io, "GradientSpoiling(", s.gradient, ", ", s.Tg, ")")

function Base.show(io::IO, ::MIME"text/plain", s::GradientSpoiling{T}) where {T}

    print(io, "GradientSpoiling{$T}:")
    print(io, "\n gradient = ", s.gradient, " G/cm")
    print(io, "\n Tg = ", s.Tg, " ms")

end

# I almost removed RFSpoiling but decided against it because it can be used to
# implement, e.g., phase cycling for bSSFP
"""
    RFSpoiling(Δθ = 117°) <: AbstractSpoiling

Represents RF spoiling, i.e., quadratically incrementing the phase of the RF
pulses from TR to TR.
"""
struct RFSpoiling{T<:Real} <: AbstractSpoiling
    Δθ::T
end

RFSpoiling() = RFSpoiling(deg2rad(117))

Base.show(io::IO, s::RFSpoiling) = print(io, "RFSpoiling(", s.Δθ, ")")
Base.show(io::IO, ::MIME"text/plain", s::RFSpoiling{T}) where {T} =
    print(io, "RFSpoiling{$T}:\n Δθ = ", s.Δθ, " rad")

"""
    RFandGradientSpoiling(grad_spoiling, rf_spoiling) <: AbstractSpoiling

Represents both RF and gradient spoiling.
"""
struct RFandGradientSpoiling{T1<:Real,T2<:Real} <: AbstractSpoiling
    gradient::GradientSpoiling{T1}
    rf::RFSpoiling{T2}
end

RFandGradientSpoiling(grad::Gradient, Tg::Real, rf::RFSpoiling) = RFandGradientSpoiling(GradientSpoiling(grad, Tg), rf)
RFandGradientSpoiling(grad::NTuple{3,Real}, Tg::Real, rf::RFSpoiling) = RFandGradientSpoiling(GradientSpoiling(grad..., Tg), rf)
RFandGradientSpoiling(gx, gy, gz, Tg, rf::RFSpoiling) = RFandGradientSpoiling(GradientSpoiling(gx, gy, gz, Tg), rf)
RFandGradientSpoiling(grad::GradientSpoiling, Δθ) = RFandGradientSpoiling(grad, RFSpoiling(Δθ))
RFandGradientSpoiling(grad::Union{<:Gradient,<:NTuple{3,Real}}, Tg, Δθ) = RFandGradientSpoiling(grad, Tg, RFSpoiling(Δθ))
RFandGradientSpoiling(grad::GradientSpoiling) = RFandGradientSpoiling(grad, RFSpoiling())
RFandGradientSpoiling(grad::Union{<:Gradient,<:NTuple{3,Real}}, Tg) = RFandGradientSpoiling(grad, Tg, RFSpoiling())
RFandGradientSpoiling(rf::Union{<:RFSpoiling,<:Real}, grad::GradientSpoiling) = RFandGradientSpoiling(grad, rf)
RFandGradientSpoiling(rf::Union{<:RFSpoiling,<:Real}, grad::Union{<:Gradient,<:NTuple{3,Real}}, Tg) = RFandGradientSpoiling(grad, Tg, rf)
RFandGradientSpoiling(rf::RFSpoiling, gx, gy, gz, Tg) = RFandGradientSpoiling(gx, gy, gz, Tg, rf)

Base.show(io::IO, s::RFandGradientSpoiling) = print(io, "RFandGradientSpoiling(", s.gradient, ", ", s.rf, ")")

function Base.show(io::IO, ::MIME"text/plain", s::RFandGradientSpoiling{T1,T2}) where {T1,T2}

    print(io, "RFandGradientSpoiling{$T1,$T2}:")
    print(io, "\n gradient = GradientSpoiling{$T1}:")
    print(io, "\n  gradient = ", s.gradient.gradient, " G/cm")
    print(io, "\n  Tg = ", s.gradient.Tg, " ms")
    print(io, "\n rf = RFSpoiling{$T2}:")
    print(io, "\n  Δθ = ", s.rf.Δθ, " rad")

end

"""
    spoiler_gradient(spoiling)

Get the `Gradient` object used for gradient spoiling.
"""
spoiler_gradient(s::GradientSpoiling) = s.gradient
spoiler_gradient(s::RFandGradientSpoiling) = spoiler_gradient(s.gradient)

"""
    spoiler_gradient_duration(spoiling)

Return the duration of the spoiler gradient (ms).
"""
spoiler_gradient_duration(::AbstractSpoiling) = 0
spoiler_gradient_duration(s::GradientSpoiling) = s.Tg
spoiler_gradient_duration(s::RFandGradientSpoiling) = spoiler_gradient_duration(s.gradient)

"""
    rfspoiling_increment(spoiling)

Return the quadratic phase increment used for RF spoiling.
"""
rfspoiling_increment(::AbstractSpoiling) = 0
rfspoiling_increment(s::RFSpoiling) = s.Δθ
rfspoiling_increment(s::RFandGradientSpoiling) = rfspoiling_increment(s.rf)

"""
    spoil(spin, spoiling, [nothing])
    spoil(spinmc, spoiling, [workspace])

Simulate gradient or ideal spoiling for the given spin. Returns `(A, B)`, such
that `A * M + B` applies spoiling to the magnetization `M`. If `B` is `nothing`
(as is the case for `IdealSpoiling`), then `A * M` applies spoiling, and if both
`A` and `B` are `nothing` (as is the case for `RFSpoiling`) then there is no
spoiling.

For `SpinMC` objects and for `GradientSpoiling` and `RFandGradientSpoiling`,
`workspace isa BlochMcConnellWorkspace`. Pass in `nothing` instead to use an
approximate matrix exponential to solve the Bloch-McConnell equation.

## Note
This function only simulates gradient or ideal spoiling, *not* RF spoiling. RF
spoiling must be implemented by updating the phase of the RF pulse(s) in your
sequence from TR to TR.

For an in-place version, see [`spoil!`](@ref).
```
"""
spoil(::AbstractSpin, ::IdealSpoiling = IdealSpoiling(), ::Any = nothing) = (idealspoiling, nothing)
spoil(::AbstractSpin, ::RFSpoiling, ::Any = nothing) = (nothing, nothing)

function spoil(
    spin::AbstractSpin,
    spoiling::AbstractSpoiling,
    workspace = spin isa Spin ? nothing : BlochMcConnellWorkspace(spin)
)

    return freeprecess(spin, spoiler_gradient_duration(spoiling), spoiler_gradient(spoiling), workspace)

end

"""
    spoil!(spin)

Apply ideal spoiling to the given spin.
"""
function spoil!(spin::Spin{T}) where {T}

    spin.M.x = zero(T)
    spin.M.y = zero(T)
    return nothing

end

function spoil!(spin::SpinMC{T}) where {T}

    for Mc in spin.M
        Mc.x = zero(T)
        Mc.y = zero(T)
    end

end

applydynamics!(spin::AbstractSpin, ::IdealSpoilingMatrix) = spoil!(spin)

"""
    spoil!(A, B, spin, spoiling, [nothing])
    spoil!(A, B, spinmc, spoiling, [workspace])

Simulate gradient or ideal spoiling, overwriting `A` and `B` (in-place version
of [`spoil`](@ref)).
"""
spoil!(A::IdealSpoilingMatrix, ::Nothing, ::AbstractSpin, ::IdealSpoiling, ::Any = nothing) = nothing
spoil!(::Nothing, ::Nothing, ::AbstractSpin, ::RFSpoiling, ::Any = nothing) = nothing

function spoil!(
    A,
    B,
    spin::AbstractSpin,
    spoiling::AbstractSpoiling,
    workspace = spin isa Spin ? nothing : BlochMcConnellWorkspace(spin)
)

    freeprecess!(A, B, spin, spoiler_gradient_duration(spoiling), spoiler_gradient(spoiling), workspace)

end


text:
include("basisfunctions.jl")
include("fit.jl")
include("designmatrix.jl")


text:
# one parameter predicates

"""Returns nonzero iff the midpoint and radius of x are both finite floating-point numbers, i.e. not infinities or NaN."""
function isfinite(x::ArbFloat{P}) where {P}
    return 0 != ccall(@libarb(arb_is_finite), Int, (Ref{ArbFloat{P}},), x)
end

"""isnan or isinf"""
function notfinite(x::ArbFloat{P}) where {P}
    return 0 == ccall(@libarb(arb_is_finite), Int, (Ref{ArbFloat{P}},), x)
end

function isnan(x::ArbFloat{P}) where {P}
    return isnan(convert(ArfFloat{P},x))
end
function notnan(x::ArbFloat{P}) where {P}
    return notnan(convert(ArfFloat{P},x))
end

function isinf(x::ArbFloat{P}) where {P}
    return isinf(convert(ArfFloat{P},x))
end
function notinf(x::ArbFloat{P}) where {P}
    return notinf(convert(ArfFloat{P},x))
end

function isposinf(x::ArbFloat{P}) where {P}
    return isposinf(convert(ArfFloat{P},x))
end
function notposinf(x::ArbFloat{P}) where {P}
    return notposinf(convert(ArfFloat{P},x))
end

function isneginf(x::ArbFloat{P}) where {P}
    return isneginf(convert(ArfFloat{P},x))
end
function notneginf(x::ArbFloat{P}) where {P}
    return notneginf(convert(ArfFloat{P},x))
end

"""midpoint(x) and radius(x) are zero"""
function iszero(x::ArbFloat{P}) where {P}
    return 0 != ccall(@libarb(arb_is_zero), Int, (Ref{ArbFloat{P}},), x)
end
"""true iff midpoint(x) or radius(x) are not zero"""
function notzero(x::ArbFloat{P}) where {P}
    return 0 == ccall(@libarb(arb_is_zero), Int, (Ref{ArbFloat{P}},), x)
end
notzero(x::T) where {T <: Real} = (x != zero(T))

"""true iff zero is not within [upperbound(x), lowerbound(x)]"""
function nonzero(x::ArbFloat{P}) where {P}
    return 0 != ccall(@libarb(arb_is_nonzero), Int, (Ref{ArbFloat{P}},), x)
end
nonzero(x::T) where {T <: Real} = (x != zero(T))

"""true iff midpoint(x) is one and radius(x) is zero"""
function isone(x::ArbFloat{P}) where {P}
    return 0 != ccall(@libarb(arb_is_one), Int, (Ref{ArbFloat{P}},), x)
end
isone(x::T) where {T <: Real} = (x == one(T))

"""true iff midpoint(x) is not one or midpoint(x) is one and radius(x) is nonzero"""
function notone(x::ArbFloat{P}) where {P}
    return 0 == ccall(@libarb(arb_is_one), Int, (Ref{ArbFloat{P}},), x)
end
notone(x::T) where {T <: Real} = (x != one(T))

"""true iff radius is zero"""
function isexact(x::ArbFloat{P}) where {P}
    return 0 != ccall(@libarb(arb_is_exact), Int, (Ref{ArbFloat{P}},), x)
end
"""true iff radius is nonzero"""
function notexact(x::ArbFloat{P}) where {P}
    return 0 == ccall(@libarb(arb_is_exact), Int, (Ref{ArbFloat{P}},), x)
end

isexact(x::T) where {T <: Integer} = true
notexact(x::T) where {T <: Integer} = false

"""true iff midpoint(x) is an integer and radius(x) is zero"""
function isinteger(x::ArbFloat{P}) where {P}
    return 0 != ccall(@libarb(arb_is_int), Int, (Ref{ArbFloat{P}},), x)
end
"""true iff midpoint(x) is not an integer or radius(x) is nonzero"""
function notinteger(x::ArbFloat{P}) where {P}
    return 0 == ccall(@libarb(arb_is_int), Int, (Ref{ArbFloat{P}},), x)
end

notinteger(x::T) where {T <: Integer} = false

"""true iff lowerbound(x) is positive"""
function ispositive(x::ArbFloat{P}) where {P}
    return 0 != ccall(@libarb(arb_is_positive), Int, (Ref{ArbFloat{P}},), x)
end

"""true iff upperbound(x) is negative"""
function isnegative(x::ArbFloat{P}) where {P}
    return  0 != ccall(@libarb(arb_is_negative), Int, (Ref{ArbFloat{P}},), x)
end

"""true iff upperbound(x) is zero or negative"""
function notpositive(x::ArbFloat{P}) where {P}
    return 0 != ccall(@libarb(arb_is_nonpositive), Int, (Ref{ArbFloat{P}},), x)
end

"""true iff lowerbound(x) is zero or positive"""
function notnegative(x::ArbFloat{P}) where {P}
    return 0 != ccall(@libarb(arb_is_nonnegative), Int, (Ref{ArbFloat{P}},), x)
end

# two parameter predicates

"""true iff midpoint(x)==midpoint(y) and radius(x)==radius(y)"""
function areequal(x::ArbFloat{P}, y::ArbFloat{P}) where {P}
    return 0 != ccall(@libarb(arb_equal), Int, (Ref{ArbFloat{P}}, Ref{ArbFloat{P}}), x, y)
end

"""true iff midpoint(x)!=midpoint(y) or radius(x)!=radius(y)"""
function notequal(x::ArbFloat{P}, y::ArbFloat{P}) where {P}
    return 0 == ccall(@libarb(arb_equal), Int, (Ref{ArbFloat{P}}, Ref{ArbFloat{P}}), x, y)
end

"""true iff x and y have a common point"""
function overlap(x::ArbFloat{P}, y::ArbFloat{P}) where {P}
    return 0 != ccall(@libarb(arb_overlaps), Int, (Ref{ArbFloat{P}}, Ref{ArbFloat{P}}), x, y)
end

"""true iff x and y have no common point"""
function donotoverlap(x::ArbFloat{P}, y::ArbFloat{P}) where {P}
    return 0 == ccall(@libarb(arb_overlaps), Int, (Ref{ArbFloat{P}}, Ref{ArbFloat{P}}), x, y)
end

"""true iff x spans (covers) all of y"""
function contains(x::ArbFloat{P}, y::ArbFloat{P}) where {P}
    return 0 != ccall(@libarb(arb_contains), Int, (Ref{ArbFloat{P}}, Ref{ArbFloat{P}}), x, y)
end

"""true iff y spans (covers) all of x"""
function iscontainedby(x::ArbFloat{P}, y::ArbFloat{P}) where {P}
    return 0 != ccall(@libarb(arb_contains), Int, (Ref{ArbFloat{P}}, Ref{ArbFloat{P}}), y, x)
end

"""true iff x does not span (cover) all of y"""
function doesnotcontain(x::ArbFloat{P}, y::ArbFloat{P}) where {P}
    return 0 != ccall(@libarb(arb_contains), Int, (Ref{ArbFloat{P}}, Ref{ArbFloat{P}}), x, y)
end

"""true iff y does not span (cover) all of x"""
function isnotcontainedby(x::ArbFloat{P}, y::ArbFloat{P}) where {P}
    return 0 == ccall(@libarb(arb_contains), Int, (Ref{ArbFloat{P}}, Ref{ArbFloat{P}}), y, x)
end

"""true if it is quite likely that the arguments indicate the same value"""
function approxeq(x::ArbFloat{P}, y::ArbFloat{P}) where {P}
    z =
        if contains(x,y)
           true
        elseif contains(y,x)
           true
        else
           x = tidy(x)
           y = tidy(y)
           contains(x,y) || contains(y,x)
        end
    return z
end

(≊)(x::ArbFloat{P}, y::ArbFloat{P}) where {P} = approxeq(x,y)
for F in (:overlap, :donotoverlap, :contains, :doesnotcontain, :iscontainedby, :isnotcontainedby, :approxeq)
  @eval ($F)(x::ArbFloat{P}, y::ArbFloat{Q}) where {P,Q} = ($F)(promote(x,y)...)
end

"""true iff x contains an integer"""
function includes_integer(x::ArbFloat{P}) where {P}
    return 0 != ccall(@libarb(arb_contains_int), Int, (Ref{ArbFloat{P}},), x)
end

"""true iff x does not contain an integer"""
function excludes_integer(x::ArbFloat{P}) where {P}
    return 0 == ccall(@libarb(arb_contains_int), Int, (Ref{ArbFloat{P}},), x)
end

"""true iff x contains zero"""
function includes_zero(x::ArbFloat{P}) where {P}
    return 0 != ccall(@libarb(arb_contains_int), Int, (Ref{ArbFloat{P}},), x)
end

"""true iff x does not contain zero"""
function excludes_zero(x::ArbFloat{P}) where {P}
    return 0 == ccall(@libarb(arb_contains_int), Int, (Ref{ArbFloat{P}},), x)
end

"""true iff x contains a positive value"""
function includes_positive(x::ArbFloat{P}) where {P}
    return 0 != ccall(@libarb(arb_contains_positive), Int, (Ref{ArbFloat{P}},), x)
end

"""true iff x contains a negative value"""
function includes_negative(x::ArbFloat{P}) where {P}
    return 0 != ccall(@libarb(arb_contains_negative), Int, (Ref{ArbFloat{P}},), x)
end

"""true iff x contains a nonpositive value"""
function includes_nonpositive(x::ArbFloat{P}) where {P}
    return 0 != ccall(@libarb(arb_contains_nonpositive), Int, (Ref{ArbFloat{P}},), x)
end

"""true iff x contains a nonnegative value"""
function includes_nonnegative(x::ArbFloat{P}) where {P}
    return 0 != ccall(@libarb(arb_contains_nonnegative), Int, (Ref{ArbFloat{ArbFloat{P}}},), x)
end


text:
# This file is a part of ParallelProcessingTools.jl, licensed under the MIT License (MIT).

using ParallelProcessingTools
using Test

using Base.Threads

@testset "threadsafe" begin

    @testset "ThreadSafeReentrantLock" begin
        tsReLock = @inferred ThreadSafeReentrantLock()
        @static if VERSION < v"1.2-DEV.28"
            lock(tsReLock)
            @test islocked(tsReLock.thread_lock)
            @test islocked(tsReLock.task_lock)
            unlock(tsReLock)
            @test !islocked(tsReLock.thread_lock)
            @test !islocked(tsReLock.task_lock)
        end
    end

    @testset "LockableValue" begin
        lv = @inferred LockableValue{Array{Float64,1}}([2.5, 3.4, 7.0])
        @test typeof(lv) <: LockableValue{Array{Float64, 1}}

        f = x::Array{Float64,1} -> x .+ 1.0
        res = @inferred broadcast(f, lv)
        @test res ≈ [3.5, 4.4, 8.0]
        res = @inferred map(f, lv)
        @test res ≈ [3.5, 4.4, 8.0]
    end

    @testset "LockableIO" begin
        lv = @inferred LockableIO(IOBuffer())
        @test typeof(lv) <: LockableIO{Base.Base.GenericIOBuffer{Array{UInt8,1}}}

        f = s -> write(s, 10)
        broadcast(f, lv)
        broadcast(seekstart, lv)
        @test read(lv, Int) == 10

        map(f, lv)
        map(seekstart, lv)
        @test read(lv, Int) == 10
        map(seekstart, lv)
        @test read!(lv, Int) == 10

        map(seekstart, lv)
        write(lv, 11)
        map(seekstart, lv)
        @test read(lv, Int) == 11
    end

    @testset "critical" begin
        @threads for i in 1:nthreads()
            @critical @info Base.Threads.threadid()
        end
    end
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    Sphere(center, radius)

A sphere with `center` and `radius`.
"""
struct Sphere{Dim,T} <: Primitive{Dim,T}
  center::Point{Dim,T}
  radius::T
end

Sphere(center::Tuple, radius) = Sphere(Point(center), radius)

paramdim(::Type{<:Sphere{Dim}}) where {Dim} = Dim - 1

center(s::Sphere) = s.center
radius(s::Sphere) = s.radius

# https://en.wikipedia.org/wiki/N-sphere#Volume_and_surface_area
function measure(s::Sphere{Dim}) where {Dim}
  r, n = s.radius, Dim
  2π^(n/2)*r^(n-1) / gamma(n/2)
end

function Base.in(p::Point, s::Sphere)
  x = coordinates(p)
  c = coordinates(s.center)
  r = s.radius
  sum(abs2, x - c) ≈ r^2
end


text:
function make_ner(;
    n_frames,
    height,
    width,
    color="black",
    words_in,
    words_out,
    margin=width/10
)

marg_width = width ÷ 2 - margin
marg_height = height ÷ 2 - margin
time_steps = floor.(Int, range(1, n_frames, length = 5))

Object(
    1:time_steps[end],
    (args...) -> JUtils.draw_text_doc(
        pos=O,
        sizex=100,
        sizey=100,
        color=color,
        n_lines=10
    )
)
n_in = length(words_in)
start_points = [(min(marg_width, marg_height) ÷ 1.5) * Point(cos(θ), -sin(θ)) for θ in (n_in > 1 ? range(π/2, -π/2, length=n_in) : [π - π/3])]
inwords = [Object(
    1:time_steps[3],
    @JShape begin
        sethue(color)
        fontsize(30)
        text(
            inword,
            start_point,
            halign=:center,
            valign=:middle
        )
    end
) for (start_point, inword) in zip(start_points, words_in)]
for (start_point, inword) in zip(start_points, inwords)
    act!(inword, Action(time_steps[2]:time_steps[3], anim_scale(0.0)))
    act!(
        inword,
        Action(time_steps[2]:time_steps[3], anim_translate(O - start_point))
    )
end

n_out = length(words_out)
outwords = [
    Object(
        time_steps[3]:time_steps[5],
        @JShape begin
            sethue(color)
            fontsize(30)
            text(
                outword,
                O,
                halign=:center,
                valign=:middle
            )
    end
) for outword in words_out
]
end_points = [
    (min(marg_width, marg_height) ÷ 1.5) * Point(cos(θ), -sin(θ))
    for θ in range(-π/2, π/2, length = n_out)
    ]
for (end_point, outword) in zip(end_points, outwords)
    act!(outword, Action(1:1, anim_scale(0.1)))
    act!(outword, Action(time_steps[1]:time_steps[2], anim_scale(1.0)))
    act!(outword, Action(time_steps[1]:time_steps[2], anim_translate(end_point)))
end
end


text:
"""
This file defined functions using Optim.jl to automatically discover circular orbits with the `BoyerLindquist`
metric parameters. This method is able to achieve approx 10^-9 numerical accuracy on the true E and Lz of a
circular orbit.
"""

# using Revise

using GeodesicTracer
using ComputedGeodesicEquations
using GeodesicBase
using CarterBoyerLindquist

using StaticArrays
using Optim

using Plots

# quality of stability function which has a minimum at stable orbits
# this is just a sum of the normalised residuals
Qs(rs) = sqrt(sum((rs ./ rs[1] .- 1.0) .^ 2) / length(rs))

# utility function for tracing a single orbit given some `vϕ`
function trace_single_geodesic(m, u, vϕ)
    v = @SVector [0.0, 0.0, 0.0, vϕ]
    tracegeodesics(
        m,
        u,
        v,
        # we pick a duration that is quite long
        # but you can toy with this
        # the method is pretty good for `r > r_isco` at even just
        # (0.0, 30.0) for the affine time
        (0.0, 300.0),
        μ = 1.0,
        # amazingly, we don't even need _that_ good tolerances
        abstol = 1e-9,
        reltol = 1e-9,
        # important, else it will spam your terminal
        verbose = false,
    )
end

# a little utility method to we don't always have to specify `u` everywhere
function geodesic_for(m, r_init, vϕ)
    u = @SVector [0.0, r_init, deg2rad(90.0), 0.0]
    trace_single_geodesic(m, u, vϕ)
end

# utility function which returns to stability after integrating a geodesic
# for some trial `vϕ`
function estimate_stability(m, u, vϕ)
    sol = trace_single_geodesic(m, u, vϕ)
    rs = selectdim(sol, 1, 6)
    Qs(rs)
end

function find_vϕ_for_orbit(m, r_init; lower_bound = 0.0, upper_bound = 0.1)
    u = @SVector [0.0, r_init, deg2rad(90.0), 0.0]
    # here we use the Optim.jl library
    res = optimize(
        # use an anonymous function to wrap our stability estimator into a univariate function
        vϕ -> estimate_stability(m, u, vϕ),
        lower_bound,
        upper_bound,
        # there are to univariate solvers, the other being `Brent`, but
        # personally I've found more success with `GoldenSection` (https://en.wikipedia.org/wiki/Golden-section_search)
        # even though it is slower
        GoldenSection(),
    )
    # uncomment this to see information about each solve
    # @show res

    # return the `vϕ` which minimized our function
    Optim.minimizer(res)
end

# this uses a sliding window for the lower and upper bounds on the `vϕ` estimate
# basically, if we overestimate the boundaries, the solution takes too many iterations and
# doesn't converge, and if we miss the estimation, we risk excluding the true solution
#
# the sliding window uses the a priori knowledge that radii closer to the origin will have
# greater `vϕ`, and so adjusts the lower and upper bounds accordingly (though still permits
# slightly lower `vϕ` to avoid possible numerical errors)
#
# it seems to work pretty well, but you may need to adjust the upper bound a little
# the tolerance is good to about a factor of 1000, so if your true velocity is 1 and you set
# the bounds to 1000, you will still find it. setting it higher than this, e.g. 10_000 will
# not converge in the set number of iterations
function find_vϕ_for_orbit_range(m, rs; lower = 0.0, upper = 0.1)
    lower = lower
    upper = upper

    # we solve backwards, from the furthest radii inwards
    # so that we can use a asymptotically flat space guess for the window (which should be metric
    # independent), and roll our way in from there
    map(reverse(rs)) do r
        vϕ = find_vϕ_for_orbit(m, r; lower_bound = lower, upper_bound = upper)
        # continuously adjust based on heuristic
        lower = vϕ * 0.99
        upper = vϕ * 2.0
        vϕ
        # reverse again to make sure the order is correct
    end |> reverse
end

# range is set from approx schwarzschild radius to whatever you like it to be
# maybe it's best just to do from r_isco outwards?
function find_orbit_range(; r_range = 2.0:0.1:10.0, a = -0.4, upper = 0.1)
    m = BoyerLindquist(M = 1.0, a = a)
    vϕs = @time find_vϕ_for_orbit_range(m, r_range; upper = upper)

    # calculate the paths
    geodesics = map(i -> geodesic_for(m, r_range[i], vϕs[i]), eachindex(r_range))

    # assemble a matrix of energies and angular momenta
    # such that `e_lz[i, 1]` gives you the energy of geodesic `i``
    e_lz = [
        f(m, geo.u[1].x[2], geo.u[1].x[1]) for geo in geodesics,
        f in (GeodesicBase.E, GeodesicBase.Lz)
    ]

    # get the time velocity
    vts = [geo.u[1].x[1][1] for geo in geodesics]

    # join time velocities in the last column
    # and angular velocities as the second last column
    hcat(e_lz, vϕs, vts)
    return geodesics
end

# test function for a single radius
# incase you want to trial indivial things
function test_single(; r_init = 10.0, a = 0.0, upper = 0.1)

    m = BoyerLindquist(M = 1.0, a = a)
    vϕ = find_vϕ_for_orbit(m, r_init; upper_bound = upper)
    sol = geodesic_for(m, r_init, vϕ)

    plt = plot(sol, vars = (8, 6), projection = :polar, range = (0.0, r_init), legend = false)
    # display(plt)
end

using Plots
gr()

pl = plot()

m = BoyerLindquist(M=1.0, a=0.998)
CarterBoyerLindquist.rms(m.M, m.a)
# m = JohannsenAD(M=1.0, a=a)
using AccretionFormulae
R_isco = AccretionFormulae.r_isco(m.a, m.M)
r_range = R_isco:3.0:30
# vϕs = @time find_vϕ_for_orbit_range(m, r_range; upper=0.1)
# geodesics = map(i -> geodesic_for(m, r_range[i], vϕs[i]), eachindex(r_range))
geodesics = find_orbit_range()

for sol in geodesics
    plot!(
        sol,
        vars=(8, 6),
        projection=:polar,
        range=(0.0, last(r_range)),
        legend=false
    )
end

pl

# pl = plot()

# m = BoyerLindquist(M=1.0, a=-0.5)
# r_range = 3.0:1.0:13.0
# vϕs = @time find_vϕ_for_orbit_range(m, r_range; upper=0.2)
# geodesics = map(i -> geodesic_for(m, r_range[i], vϕs[i]), eachindex(r_range))

# for sol in geodesics
#     plot!(
#         sol,
#         vars=(8, 6),
#         projection=:polar,
#         range=(0.0, last(r_range)),
#         legend=false
#     )
# end

# plot!(
#     _ -> GeodesicBase.inner_radius(m),
#     0.0:0.01:2π,
#     c=:black,
#     lw=5
# )
# plot!(
#     _ -> CarterBoyerLindquist.rms(m.M, m.a),
#     0.0:0.01:2π,
#     c=:black,
#     lw=2,
#     ls=:dot
# )

# pl

# savefig("example_circs.svg")


text:
restr_path = Base.source_dir() * "/SampleFiles/Restricted/"
if isdir(restr_path) == false
  for i in ["/data/", "/data2"]
    restr_dir = i * "SeisIO-TestFiles/SampleFiles/Restricted"
    if isdir(restr_dir)
      try
        run(`cp -r $restr_dir SampleFiles/`)
        println("copied SampleFiles/Restricted/ from ", i)
      catch err
        @warn(err)
      end
    end
  end
else
  println("nothing to copy, SampleFiles/Restricted/ exists")
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    DensityRatioWeighting(tdata, [vars]; [options])

Density ratio weights based on empirical distribution of
variables in target data `tdata`. Default to all variables.

## Optional parameters

* `estimator` - Density ratio estimator (default to `LSIF()`)
* `optlib`    - Optimization library (default to `default_optlib(estimator)`)

### Notes

Estimators from `DensityRatioEstimation.jl` are supported.
"""
struct DensityRatioWeighting <: WeightingMethod
  tdata
  vars
  dre
  optlib
end

function DensityRatioWeighting(tdata, vars=nothing; estimator=LSIF(),
                               optlib=default_optlib(estimator))
  validvars = collect(name.(variables(tdata)))
  wvars = isnothing(vars) ? validvars : vars
  @assert wvars ⊆ validvars "invalid variables ($wvars) for spatial data"
  DensityRatioWeighting(tdata, wvars, estimator, optlib)
end

function weight(sdata, method::DensityRatioWeighting)
  # retrieve method parameters
  tdata  = method.tdata
  vars   = method.vars
  dre    = method.dre
  optlib = method.optlib

  @assert vars ⊆ name.(variables(sdata)) "invalid variables ($vars) for spatial data"

  # numerator and denominator samples
  Ωnu = view(tdata, vars)
  Ωde = view(sdata, vars)
  xnu = collect(Tables.rows(values(Ωnu)))
  xde = collect(Tables.rows(values(Ωde)))

  # perform denstiy ratio estimation
  ratios = densratio(xnu, xde, dre, optlib=optlib)

  SpatialWeights(domain(sdata), ratios)
end


text:
using FFTW, LinearAlgebra, Plots
gr()

include("src/dataset.jl")
include("src/micmac.jl")

const epsilons = [10.0^(-i) for i in 0:6]
const xmin     = 0
const xmax     = 2π
const T        = 2π
const size_x   = [64]
const size_tau = [64]
const Tfinal   = 0.25
const nb_dt    = 5 # different values of dt

tabdt  = zeros(Float64, nb_dt)
taberr = zeros(Float64, (length(epsilons), nb_dt))

etime = @elapsed for nx in size_x, ntau in size_tau

    println(" nx                 : $nx ")
    println(" ntau               : $ntau ")

    for (kk, epsilon) in enumerate(epsilons)

        print(" $epsilon: ")

        data = CosSin(xmin, xmax, nx, epsilon, T, Tfinal)

        for hh in 0:nb_dt-1

            println("dt = $(2.0^(-hh)) ")

            dtmicmac = 2.0^(-hh) * data.Tfinal / 16

            solver = MicMac(data, ntau)

            u, v = solve(solver, dtmicmac)

            tabdt[hh+1] = dtmicmac

            err = compute_error(u,v,data)
            println(" error = $err ")
            taberr[kk,hh+1] = err

        end

        println()

    end

end

println("Elapsed time :", etime)

p = plot(layout=(1,2))
xlabel!(p[1,1],"dt")
ylabel!(p[1,1],"error")

for j in 1:size(taberr)[1]
    plot!(p[1,1],
          tabdt, taberr[j, :],
          markershape=:circle,
          xaxis=:log,
          yaxis=:log,
          label="epsilon=$(epsilons[j])")
end

for j in 1:size(taberr)[2]
    plot!(p[1,2],
          epsilons, taberr[:, j],
          markershape=:circle,
          xaxis=:log,
          yaxis=:log,
          label="dt=$(tabdt[j])")
end

xlabel!(p[1,2],"epsilon")

p


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

mutable struct ContainerServiceProperties <: SwaggerModel
    provisioningState::Any # spec type: Union{ Nothing, String } # spec name: provisioningState
    orchestratorProfile::Any # spec type: Union{ Nothing, ContainerServiceOrchestratorProfile } # spec name: orchestratorProfile
    customProfile::Any # spec type: Union{ Nothing, ContainerServiceCustomProfile } # spec name: customProfile
    servicePrincipalProfile::Any # spec type: Union{ Nothing, ContainerServiceServicePrincipalProfile } # spec name: servicePrincipalProfile
    masterProfile::Any # spec type: Union{ Nothing, ContainerServiceMasterProfile } # spec name: masterProfile
    agentPoolProfiles::Any # spec type: Union{ Nothing, Vector{ContainerServiceAgentPoolProfile} } # spec name: agentPoolProfiles
    windowsProfile::Any # spec type: Union{ Nothing, ContainerServiceWindowsProfile } # spec name: windowsProfile
    linuxProfile::Any # spec type: Union{ Nothing, ContainerServiceLinuxProfile } # spec name: linuxProfile
    diagnosticsProfile::Any # spec type: Union{ Nothing, ContainerServiceDiagnosticsProfile } # spec name: diagnosticsProfile

    function ContainerServiceProperties(;provisioningState=nothing, orchestratorProfile=nothing, customProfile=nothing, servicePrincipalProfile=nothing, masterProfile=nothing, agentPoolProfiles=nothing, windowsProfile=nothing, linuxProfile=nothing, diagnosticsProfile=nothing)
        o = new()
        validate_property(ContainerServiceProperties, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        validate_property(ContainerServiceProperties, Symbol("orchestratorProfile"), orchestratorProfile)
        setfield!(o, Symbol("orchestratorProfile"), orchestratorProfile)
        validate_property(ContainerServiceProperties, Symbol("customProfile"), customProfile)
        setfield!(o, Symbol("customProfile"), customProfile)
        validate_property(ContainerServiceProperties, Symbol("servicePrincipalProfile"), servicePrincipalProfile)
        setfield!(o, Symbol("servicePrincipalProfile"), servicePrincipalProfile)
        validate_property(ContainerServiceProperties, Symbol("masterProfile"), masterProfile)
        setfield!(o, Symbol("masterProfile"), masterProfile)
        validate_property(ContainerServiceProperties, Symbol("agentPoolProfiles"), agentPoolProfiles)
        setfield!(o, Symbol("agentPoolProfiles"), agentPoolProfiles)
        validate_property(ContainerServiceProperties, Symbol("windowsProfile"), windowsProfile)
        setfield!(o, Symbol("windowsProfile"), windowsProfile)
        validate_property(ContainerServiceProperties, Symbol("linuxProfile"), linuxProfile)
        setfield!(o, Symbol("linuxProfile"), linuxProfile)
        validate_property(ContainerServiceProperties, Symbol("diagnosticsProfile"), diagnosticsProfile)
        setfield!(o, Symbol("diagnosticsProfile"), diagnosticsProfile)
        o
    end
end # type ContainerServiceProperties

const _property_map_ContainerServiceProperties = Dict{Symbol,Symbol}(Symbol("provisioningState")=>Symbol("provisioningState"), Symbol("orchestratorProfile")=>Symbol("orchestratorProfile"), Symbol("customProfile")=>Symbol("customProfile"), Symbol("servicePrincipalProfile")=>Symbol("servicePrincipalProfile"), Symbol("masterProfile")=>Symbol("masterProfile"), Symbol("agentPoolProfiles")=>Symbol("agentPoolProfiles"), Symbol("windowsProfile")=>Symbol("windowsProfile"), Symbol("linuxProfile")=>Symbol("linuxProfile"), Symbol("diagnosticsProfile")=>Symbol("diagnosticsProfile"))
const _property_types_ContainerServiceProperties = Dict{Symbol,String}(Symbol("provisioningState")=>"String", Symbol("orchestratorProfile")=>"ContainerServiceOrchestratorProfile", Symbol("customProfile")=>"ContainerServiceCustomProfile", Symbol("servicePrincipalProfile")=>"ContainerServiceServicePrincipalProfile", Symbol("masterProfile")=>"ContainerServiceMasterProfile", Symbol("agentPoolProfiles")=>"Vector{ContainerServiceAgentPoolProfile}", Symbol("windowsProfile")=>"ContainerServiceWindowsProfile", Symbol("linuxProfile")=>"ContainerServiceLinuxProfile", Symbol("diagnosticsProfile")=>"ContainerServiceDiagnosticsProfile")
Base.propertynames(::Type{ ContainerServiceProperties }) = collect(keys(_property_map_ContainerServiceProperties))
Swagger.property_type(::Type{ ContainerServiceProperties }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ContainerServiceProperties[name]))}
Swagger.field_name(::Type{ ContainerServiceProperties }, property_name::Symbol) =  _property_map_ContainerServiceProperties[property_name]

function check_required(o::ContainerServiceProperties)
    (getproperty(o, Symbol("masterProfile")) === nothing) && (return false)
    (getproperty(o, Symbol("agentPoolProfiles")) === nothing) && (return false)
    (getproperty(o, Symbol("linuxProfile")) === nothing) && (return false)
    true
end

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


text:
using RDatasets, DataArrays, Gadfly

plot(dataset("plm", "Cigar"), x="Sales", y="Year", Geom.bar(orientation=:horizontal))


text:
"""
GenX: An Configurable Capacity Expansion Model
Copyright (C) 2021,  Massachusetts Institute of Technology
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
A complete copy of the GNU General Public License v2 (GPLv2) is available
in LICENSE.txt.  Users uncompressing this from an archive may not have
received this license file.  If not, see <http://www.gnu.org/licenses/>.
"""

function write_esr_prices(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)
	dfESR = DataFrame(ESR_Price = convert(Array{Union{Missing, Float64}}, dual.(EP[:cESRShare])))
	if setup["ParameterScale"] == 1
		dfESR[!,:ESR_Price] = dfESR[!,:ESR_Price] * ModelScalingFactor # Converting MillionUS$/GWh to US$/MWh
	end
	CSV.write(joinpath(path, "ESR_prices.csv"), dfESR)
	return dfESR
end

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM


text:
export AbstractLearner, AbstractModelBasedLearner, AbstractModelFreeLearner,
       update!

"""
    AbstractLearner

Supertype for learners. A learner takes in a state and returns an action.
"""
abstract type AbstractLearner end
function update! end

abstract type AbstractModelFreeLearner <: AbstractLearner end
abstract type AbstractMonteCarloLearner <: AbstractModelFreeLearner end

abstract type AbstractModelBasedLearner <: AbstractLearner end
function update! end