text:
#
# Date created: 2022-03-14
# Author: aradclif
#
#
############################################################################################
function sizeblock(N::Int)
    block = Expr(:block)
    for d = 1:N
        ex = Expr(:(=), Symbol(:D_, d), Expr(:call, :size, :A, d))
        push!(block.args, ex)
    end
    block
end

function sizeproductsblock(N::Int)
    block = Expr(:block)
    for k = 3:N
        ex = Expr(:(=), Symbol(:D_, ntuple(identity, k - 1)...),
                  Expr(:call, :*, ntuple(d -> Symbol(:D_, d), k - 1)...))
        push!(block.args, ex)
    end
    block
end

function sumprodliteral(N::Int)
    e = Expr(:call, :+, Symbol(:i_, 1))
    for k = 2:N
        if k == 2
            ex = Expr(:call, :-, Expr(:call, :*, Symbol(:D_, 1), Symbol(:i_, 2)), Symbol(:D_, 1))
            push!(e.args, ex)
        else
            ex = Expr(:call, :-,
                      Expr(:call, :*, ntuple(d -> Symbol(:D_, d), k - 1)..., Symbol(:i_, k)),
                      Expr(:call, :*, ntuple(d -> Symbol(:D_, d), k - 1)...))
            push!(e.args, ex)
        end
    end
    e
end

function sumprodprecomputed(N::Int)
    e = Expr(:call, :+, Symbol(:i_, 1))
    for k = 2:N
        if k == 2
            ex = Expr(:call, :-, Expr(:call, :*, Symbol(:D_, 1), Symbol(:i_, 2)), Symbol(:D_, 1))
            push!(e.args, ex)
        else
            ex = Expr(:call, :-,
                      Expr(:call, :*, Symbol(:D_, ntuple(identity, k - 1)...), Symbol(:i_, k)),
                      Symbol(:D_, ntuple(identity, k - 1)...))
            push!(e.args, ex)
        end
    end
    e
end

function sumprodprecomputed2(N::Int)
    e = Expr(:call, :+, Symbol(:i_, 1))
    for k = 2:N
        if k == 2
            ex = Expr(:call, :*, Symbol(:D_, 1), Symbol(:i_, 2))
            push!(e.args, ex)
        else
            ex = Expr(:call, :*, Symbol(:D_, ntuple(identity, k - 1)...), Symbol(:i_, k))
            push!(e.args, ex)
        end
    end
    e
end

function sumprodconstant(N::Int)
    Expr(:(=), :D_sp,
         Expr(:call, :+, ntuple(d -> Expr(:call, :*, ntuple(i -> Symbol(:D_, i), d)...), N - 1)...))
    # # # Equivalent, but perhaps more clear
    # e = Expr(:(=), :D_sp)
    # r = Expr(:call, :+)
    # for k = 1:(N - 1)
    #     ex = Expr(:call, :*, ntuple(d -> Symbol(:D_, d), k)...)
    #     push!(r.args, ex)
    # end
    # push!(e.args, r)
    # e
end

function sumprodconstant2(N::Int)
    Expr(:(=), :D_sp, Expr(:call, :+, ntuple(d -> Symbol(:D_, ntuple(identity, d)...), N - 1)...))
end

prog = "a == b ? 5 : 10"
ex = Meta.parse(prog)
Meta.show_sexpr(ex)
# Expr(:?, Expr(:call, :(==), :a, :b), 5, 10)
prog2 = "if a == b
    5
    else
    10
    end"
ex2 = Meta.parse(prog2)
Meta.show_sexpr(ex2)
prog3 = "ifelse(a == b, 5, 10)"
ex3 = Meta.parse(prog3)
Meta.show_sexpr(ex3)
prog4 = "a = -a"
ex4 = Meta.parse(prog4)
Meta.show_sexpr(ex4)

function compareblock(N::Int, D)
    block = Expr(:block)
    params = D.parameters
    a = Expr(:ref, :A, ntuple(d -> Symbol(:i_, d), N)...)
    b = Expr(:ref, :B, ntuple(d -> params[d] == Static.One ? 1 : Symbol(:i_, d), N)...)
    c = Expr(:ref, :C, ntuple(d -> params[d] == Static.One ? 1 : Symbol(:i_, d), N)...)
    # y = Expr(:(=), :yes, Expr(:call, :(==), a, b))
    # e = Expr(:(=), c, Expr(:if, y, sumprodliteral(N), c))
    # push!(block.args, y)
    e = Expr(:(=), c, Expr(:if, Expr(:call, :(==), a, b), sumprodliteral(N), c))
    push!(block.args, e)
    block
end

function compareblock2(N::Int, D)
    block = Expr(:block)
    params = D.parameters
    a = Expr(:ref, :A, ntuple(d -> Symbol(:i_, d), N)...)
    b = Expr(:ref, :B, ntuple(d -> params[d] == Static.One ? 1 : Symbol(:i_, d), N)...)
    c = Expr(:ref, :C, ntuple(d -> params[d] == Static.One ? 1 : Symbol(:i_, d), N)...)
    d = sumprodprecomputed2(N)
    push!(d.args, :D_sp)
    e = Expr(:(=), c, Expr(:if, Expr(:call, :(==), a, b), d, c))
    push!(block.args, e)
    block
end

# function compareblock3(N::Int, D)
#     block = Expr(:block)
#     params = D.parameters
#     a = Expr(:ref, :A, ntuple(d -> Symbol(:i_, d), N)...)
#     b = Expr(:ref, :B, ntuple(d -> params[d] == Static.One ? 1 : Symbol(:i_, d), N)...)
#     c = Expr(:ref, :C, ntuple(d -> params[d] == Static.One ? 1 : Symbol(:i_, d), N)...)
#     d = sumprodprecomputed2(N)
#     push!(d.args, :D_sp)
#     # e = Expr(:if, Expr(:call, :(==), a, b), Expr(:(=), c, d))
#     j = Expr(:(=), :j, c)
#     e = Expr(:(=), c, Expr(:if, Expr(:call, :(==), a, b), d, :j))
#     push!(block.args, j)
#     push!(block.args, e)
#     block
# end

function compareblock4(N::Int, D)
    block = Expr(:block)
    params = D.parameters
    a = Expr(:ref, :A, ntuple(d -> Symbol(:i_, d), N)...)
    b = Expr(:ref, :B, ntuple(d -> params[d] == Static.One ? 1 : Symbol(:i_, d), N)...)
    c = Expr(:ref, :C, ntuple(d -> Symbol(:i_, d), N)...)
    e = Expr(:(=), c, Expr(:call, :(==), a, b))
    push!(block.args, e)
    block
end

function compareblock6(N::Int, D)
    block = Expr(:block)
    params = D.parameters
    a = Expr(:ref, :A, ntuple(d -> Symbol(:i_, d), N)...)
    b = Expr(:ref, :B, ntuple(d -> params[d] == Static.One ? 1 : Symbol(:i_, d), N)...)
    c = Expr(:ref, :C, ntuple(d -> params[d] == Static.One ? 1 : Symbol(:i_, d), N)...)
    d = sumprodprecomputed2(N)
    push!(d.args, :D_sp)
    e = Expr(:(=), c, Expr(:call, :ifelse, Expr(:call, :(==), a, b), d, c))
    push!(block.args, e)
    block
end

# p. 26
function outerloopgen(N::Int, D)
    loops = Expr(:for)
    block = Expr(:block)
    params = D.parameters
    for d = N:-1:1
        if params[d] != Static.One
            ex = Expr(:(=), Symbol(:i_, d), Expr(:call, :axes, :A, d))
            push!(block.args, ex)
        end
    end
    push!(loops.args, block)
    loops
end
function innerloopgen(N::Int, D)
    loops = Expr(:for)
    block = Expr(:block)
    params = D.parameters
    for d = N:-1:1
        if params[d] == Static.One
            ex = Expr(:(=), Symbol(:i_, d), Expr(:call, :axes, :A, d))
            push!(block.args, ex)
        end
    end
    push!(loops.args, block)
    loops
end

function innerpost(N::Int, D)
    params = D.parameters
    block = Expr(:block)
    b = Expr(:ref, :B, ntuple(d -> params[d] == Static.One ? 1 : Symbol(:i_, d), N)...)
    c = Expr(:ref, :C, ntuple(d -> params[d] == Static.One ? 1 : Symbol(:i_, d), N)...)
    e1 = Expr(:(=), b, :m)
    e2 = Expr(:(=), c, :j)
    push!(block.args, e1)
    push!(block.args, e2)
    block
end

function compareblock5(N::Int, D)
    block = Expr(:block)
    params = D.parameters
    a = Expr(:ref, :A, ntuple(d -> Symbol(:i_, d), N)...)
    d = sumprodprecomputed2(N)
    push!(d.args, :D_sp)
    yₑ = Expr(:(=), :y, Expr(:call, :(>), a, :m))
    mₑ = Expr(:(=), :m, Expr(:if, :y, a, :m))
    jₑ = Expr(:(=), :j, Expr(:if, :y, d, :j))
    push!(block.args, yₑ)
    push!(block.args, mₑ)
    push!(block.args, jₑ)
    block
end

function findmax5_quote(N::Int, D)
    block1 = sizeblock(N)
    block2 = sizeproductsblock(N)
    block3 = Expr(:block, sumprodconstant(N), Expr(:(=), :D_sp, Expr(:call, :-, :D_sp)))
    outerloops = outerloopgen(N, D)
    block4 = Expr(:block, Expr(:(=), :j, 1), Expr(:(=), :m, Expr(:call, :typemin, :T)))
    # push!(outerloops.args, block4)
    innerloops = innerloopgen(N, D)
    block5 = compareblock5(N, D)
    push!(innerloops.args, block5)
    push!(block4.args, innerloops)
    # push!(outerloops.args, block4)
    # push!(outerloops.args, innerloops)
    block6 = innerpost(N, D)
    # push!(outerloops.args, block6)
    push!(block4.args, block6.args...)
    push!(outerloops.args, block4)
    return quote
        $block1
        $block2
        $block3
        $outerloops
    end
end

function findequal_quote(N::Int, D)
    loops = loopgen(N)
    block1 = sizeblock(N)
    block2 = compareblock(N, D)
    push!(loops.args, block2)
    return quote
        $block1
        $loops
    end
end

function findequal_quote2(N::Int, D)
    loops = loopgen(N)
    block1 = sizeblock(N)
    block2 = sizeproductsblock(N)
    block3 = Expr(:block, sumprodconstant(N), Expr(:(=), :D_sp, Expr(:call, :-, :D_sp)))
    block4 = compareblock2(N, D)
    push!(loops.args, block4)
    return quote
        $block1
        $block2
        $block3
        @turbo $loops
    end
end
# function findequal_quote3(N::Int, D)
#     loops = loopgen(N)
#     block1 = sizeblock(N)
#     block2 = sizeproductsblock(N)
#     block3 = Expr(:block, sumprodconstant(N), Expr(:(=), :D_sp, Expr(:call, :-, :D_sp)))
#     block4 = compareblock3(N, D)
#     push!(loops.args, block4)
#     return quote
#         $block1
#         $block2
#         $block3
#         @turbo $loops
#     end
# end

function findequal_quote4(N::Int, D)
    loops = loopgen(N)
    block = compareblock4(N, D)
    push!(loops.args, block)
    return quote
        @turbo $loops
    end
end

function findequal_quote6(N::Int, D)
    loops = loopgen(N)
    block1 = sizeblock(N)
    block2 = sizeproductsblock(N)
    block3 = Expr(:block, sumprodconstant(N), Expr(:(=), :D_sp, Expr(:call, :-, :D_sp)))
    block4 = compareblock6(N, D)
    push!(loops.args, block4)
    return quote
        $block1
        $block2
        $block3
        @turbo $loops
    end
end

@generated function findequal!(C::AbstractArray{Tₒ, N}, A::AbstractArray{T, N},
                               B::AbstractArray{T, N}, dims::D) where {Tₒ, T, N, D}
    findequal_quote2(N, D)
end
# @generated function findequal3!(C::AbstractArray{Tₒ, N}, A::AbstractArray{T, N},
#                                B::AbstractArray{T, N}, dims::D) where {Tₒ, T, N, D}
#     findequal_quote3(N, D)
# end
@generated function findequal4!(C::AbstractArray{Tₒ, N}, A::AbstractArray{T, N},
                                B::AbstractArray{T, N}, dims::D) where {Tₒ, T, N, D}
    findequal_quote4(N, D)
end

@generated function findmax5!(C::AbstractArray{Tₒ, N}, A::AbstractArray{T, N},
                              B::AbstractArray{T, N}, dims::D) where {Tₒ, T, N, D}
    findmax5_quote(N, D)
end

@generated function findequal6!(C::AbstractArray{Tₒ, N}, A::AbstractArray{T, N},
                                B::AbstractArray{T, N}, dims::D) where {Tₒ, T, N, D}
    findequal_quote6(N, D)
end

A = reshape([1:(4*3*5);], 4, 3, 5);
A = rand(1:10, 4, 3, 5);
dims = (2,);
Dᴮ′ = ntuple(d -> d ∈ dims ? StaticInt(1) : size(A, d), ndims(A));
findequal_quote(ndims(A), typeof(Dᴮ′))
findequal_quote2(ndims(A), typeof(Dᴮ′))
findequal_quote3(ndims(A), typeof(Dᴮ′))
findequal_quote6(ndims(A), typeof(Dᴮ′))
findmax5_quote(ndims(A), typeof(Dᴮ′))

# C0 = deepcopy(C);
B = maximum(A, dims=dims);
C = ones(Int, size(B));

CartesianIndices(A)[C] == argmax(A, dims=dims)

compareblock(ndims(A), typeof(Dᴮ′))
compareblock2(ndims(A), typeof(Dᴮ′))
compareblock3(ndims(A), typeof(Dᴮ′))

function lvfindmax(A::AbstractArray{T, N}, dims::NTuple{M, Int}) where {T, N, M}
    B = lvmaximum(A, dims=dims)
    Dᴮ′ = ntuple(d -> d ∈ dims ? StaticInt(1) : size(A, d), N)
    C = ones(Int, size(B))
    findequal!(C, A, B, Dᴮ′)
    C
end
C2 = lvfindmax(A, dims)

CartesianIndices(A)[C2] == argmax(A, dims=dims)
LinearIndices(A)[argmax(A, dims=dims)]

# function lvfindmax3(A::AbstractArray{T, N}, dims::NTuple{M, Int}) where {T, N, M}
#     B = lvmaximum(A, dims=dims)
#     Dᴮ′ = ntuple(d -> d ∈ dims ? StaticInt(1) : size(A, d), N)
#     C = ones(Int, size(B))
#     findequal3!(C, A, B, Dᴮ′)
#     C
# end
# C3 = lvfindmax3(A, dims)

# function lvfindmax4(A::AbstractArray{T, N}, dims::NTuple{M, Int}) where {T, N, M}
#     B = lvmaximum(A, dims=dims)
#     Dᴮ′ = ntuple(d -> d ∈ dims ? StaticInt(1) : size(A, d), N)
#     C = similar(A, Bool)
#     findequal4!(C, A, B, Dᴮ′)
#     C
# end
# C4 = lvfindmax4(A, dims)
# reshape(findall(C4), size(B)) == argmax(A, dims=dims) == reshape(CartesianIndices(A)[C4], size(B))

function lvfindmax6(A::AbstractArray{T, N}, dims::NTuple{M, Int}) where {T, N, M}
    B = maximum(A, dims=dims)
    Dᴮ′ = ntuple(d -> d ∈ dims ? StaticInt(1) : size(A, d), N)
    C = ones(Int, size(B))
    findequal6!(C, A, B, Dᴮ′)
    B, CartesianIndices(A)[C]
end
B6, C6 = lvfindmax6(A, dims)
lvfindmax6(A, dims) == findmax(A, dims=dims)

function lvfindmax5(A::AbstractArray{T, N}, dims::NTuple{M, Int}) where {T, N, M}
    Dᴮ′ = ntuple(d -> d ∈ dims ? StaticInt(1) : size(A, d), N)
    B = similar(A, Dᴮ′)
    C = similar(B, Int)
    findmax5!(C, A, B, Dᴮ′)
    B, CartesianIndices(A)[C]
end
C5 = lvfindmax5(A, dims)
CartesianIndices(A)[C5] == argmax(A, dims=dims)
lvfindmax5(A, dims)
findmax(A, dims=dims)
lvfindmax5(A, dims) == findmax(A, dims=dims)

function lvfindequal(A::AbstractArray{T, 3}, B::AbstractArray{T, 3}) where {T}
    C = similar(B, Int)
    D_1 = size(A, 1)
    D_2 = size(A, 2)
    D_3 = size(A, 3)
    D_12 = D_1 * D_2
    D_sp = (*)(D_1) + D_1 * D_2
    D_sp = -D_sp
    @turbo for i_3 = axes(A, 3), i_1 = axes(A, 1)
        j = 1
        m = B[i_1, 1, i_3]
        for i_2 = axes(A, 2)
            # j = if A[i_1, i_2, i_3] == m
            #     i_1 + D_1 * i_2 + D_12 * i_3 + D_sp
            # else
            #     j
            # end
            yes = A[i_1, i_2, i_3] == m
            j = yes ? i_1 + D_1 * i_2 + D_12 * i_3 + D_sp : j
        end
        C[i_1, 1, i_3] = j
    end
    C
end
C5 = lvfindequal(A, B)
CartesianIndices(A)[C5] == argmax(A, dims=dims)

# an experiment -- seems to work
function lvfindmax5(A::AbstractArray{T, 3}, dims::NTuple{M, Int}) where {T, M}
    Dᴮ′ = ntuple(d -> d ∈ dims ? StaticInt(1) : size(A, d), 3)
    B = similar(A, Dᴮ′)
    C = similar(B, Int)
    D_1 = size(A, 1)
    D_2 = size(A, 2)
    D_3 = size(A, 3)
    D_12 = D_1 * D_2
    D_sp = (*)(D_1) + D_1 * D_2
    D_sp = -D_sp
    @turbo for i_3 = axes(A, 3), i_1 = axes(A, 1)
        j = 1
        m = typemin(T)
        for i_2 = axes(A, 2)
            # j = if A[i_1, i_2, i_3] == m
            #     i_1 + D_1 * i_2 + D_12 * i_3 + D_sp
            # else
            #     j
            # end
            yes = A[i_1, i_2, i_3] > m
            m = yes ? A[i_1, i_2, i_3] : m
            j = yes ? i_1 + D_1 * i_2 + D_12 * i_3 + D_sp : j
        end
        B[i_1, 1, i_3] = m
        C[i_1, 1, i_3] = j
    end
    B, C
end

B5, C5 = lvfindmax5(A, (1,))
A[C5] == B5
b5, c5 = findmax(A, dims=1)
b5 == B5
c5 == CartesianIndices(A)[C5]


text:
# using Base.Test
# using MIToS.MSA
# using MIToS.Information

print("""

Estimation of measures on a MSA
===============================
""")

print("""
This is the example of MI(X, Y)/H(X, Y) from:
Gao, H., Dou, Y., Yang, J., & Wang, J. (2011). New methods to measure residues coevolution in proteins. BMC bioinformatics, 12(1), 206.
""")
let aln = read(joinpath(pwd(), "data", "Gaoetal2011.fasta"), FASTA), result = Float64[ 0     0     0     0     0     0
                                                                     0     0     0     0     0     0
                                                                     0     0     0     1     1     0.296
                                                                     0     0     1     0     1     0.296
                                                                     0     0     1     1     0     0.296
                                                                     0     0     0.296 0.296 0.296 0 ]

  @test_approx_eq_eps triu(estimateincolumns(aln, ResidueCount{Int, 2, false}, MutualInformationOverEntropy{Float64}()), 1) triu(result, 1) 0.0001
  @test_approx_eq_eps triu(estimateincolumns(aln, Int, ResidueProbability{Float64, 2, false}, MutualInformationOverEntropy{Float64}()), 1) triu(result, 1) 0.0001
  @test_approx_eq_eps triu(estimateincolumns(aln, Float64, ResidueProbability{Float64, 2, false}, MutualInformationOverEntropy{Float64}()), 1) triu(result, 1) 0.0001
  @test_approx_eq_eps triu(estimateincolumns(aln, BigFloat, ResidueProbability{BigFloat, 2, false}, MutualInformationOverEntropy{BigFloat}()), 1) triu(result, 1) 0.0001

  @test_approx_eq_eps estimateincolumns(aln, ResidueCount{Int, 2, false}, MutualInformationOverEntropy{Float64}(),
                                        zero(AdditiveSmoothing{Int}), NoClustering(), false, zero(Float64)) result 0.0001
  @test_approx_eq_eps estimateincolumns(aln, Int, ResidueProbability{Float64, 2, false}, MutualInformationOverEntropy{Float64}(),
                                        zero(AdditiveSmoothing{Int}), NoClustering(), false, zero(Float64)) result 0.0001
  @test_approx_eq_eps estimateincolumns(aln, Float64, ResidueProbability{Float64, 2, false}, MutualInformationOverEntropy{Float64}(),
                                        zero(AdditiveSmoothing{Float64}), NoClustering(), false, zero(Float64)) result 0.0001
  @test_approx_eq_eps estimateincolumns(aln, BigFloat, ResidueProbability{BigFloat, 2, false}, MutualInformationOverEntropy{BigFloat}(),
                                        zero(AdditiveSmoothing{BigFloat}), NoClustering(), false, zero(BigFloat)) result 0.0001
end

ext: jl


text:
#   Unit "plots" of the FourierAnalysis Package for julia language
#   v 0.0.1 - last update 7th of August 2019
#
#   MIT License
#   Copyright (c) 2019, Marco Congedo, CNRS, Grenobe, France:
#   https://sites.google.com/site/marcocongedo/home

# ? CONTENTS :
#   This unit implements utilities fr plotting data produced by
#   Fourier Analysis

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

"""
```
plot( h      :: Taper,
    showkind :: Bool = false;
   args...)
```

Plot the data of a [Taper](@ref) object, that is, the tapering
window(s) hold in its `.y` field.
This is a single tapering window for all kind of tapering windows,
but for Slepian's discrete prolate spheroid sequences,
which are multi-tapering.

If optional keyword argument `showkind` is true, the output of the
[`taperinfo`](@ref) function is printed as title of the plot.
By defaut `showkind` is false.

The function allows to pass any other optional keyword argument
of the standard Julia plot function.

**Examples**:
```
using FourierAnalysis, Plots
sr, t = 128, 128
H=taper(hamming, t)

# plot tapering window
plot(H; title="my own title") # title is a Plot.jl argument
plot(H; showkind=true)

# plot Slepian's dpss
plot(slepians(sr, t*4, 2); showkind=true)
plot(slepians(sr, t*16, 1); title="my own title")
```
"""
plot( h      :: Taper;
    showkind :: Bool = false,
   args...)  =
  plot(h.y;
    title=showkind ? taperinfo(h) : "", gridcolor=:white, labels="", args...)

"""
```
plot(S      :: Spectra;
    maxf    :: Union{Real, Int} = 0,
    space   :: Int              = 4,
    ylabel  :: String = "Power (\\muV²)",
   args...)
```
Produce the line plot of the data in a [Spectra](@ref) object,
that is, plot all its spectra with separate lines.

Optional keyword argument `maxf` limits the plot on the right
to the given frequency (in Hz).

Optional keyword argument `space` (in Hz) determines the spacing between
frequency ticks and labels (x-axis).

The function allows to pass any other optional keyword argument
of the standard Julia plot function.

**See also**: [`bar`](@ref)

**Examples**:
```
using FourierAnalysis, Plots
sr, t, f, a = 128, 128, 16, 0.5
# create a sinusoidal at 16Hz superimposed to white noise
v=sinusoidal(a, f, sr, t*16, 0) + randn(t*16)
# create a data matrix
X=broadcast(+, v, randn(t*16, 3)) * randn(3, 3)
# compute spectra using hamming tapering window
H=taper(hamming, t)
S=spectra(X, sr, t; tapering=H)

# plot spectra
plot(S)
plot(S; maxf=32)
plot(S; maxf=32, space=2)

# compute and plot amplitude spectra
S=spectra(X, sr, t; tapering=H, func=sqrt)
plot(S; maxf=32, space=2, ylabel="Amplitude (\\muV)")

```
"""
plot(S      :: Spectra;
    maxf    :: Union{Real, Int} = 0,
    space   :: Int              = 4,
    ylabel  :: String = "Power (\\muV²)",
   args...) =
  plot(S.y[1:(maxf==0 ? f2b(S.sr/2, S.sr, S.wl; DC=S.DC) : f2b(maxf, S.sr, S.wl; DC=S.DC)), :];
     _spectraplotargs(S;
                  maxf   = maxf,
                  space  = space,
                  ylabel = ylabel)...,
     args...)

 """
 ```
 bar(S      :: Spectra;
     maxf    :: Union{Real, Int} = 0,
     space   :: Int              = 4,
     ylabel  :: String = "Power (\\muV²)",
    args...)
```

 Produce the bar plot of the data in a [Spectra](@ref) object,
 that is, plot all its spectra with separate bars.
 Use the same syntax as the [`plot`](@ref) fuction,
 which produces a line plot instead.
 Bar plots for spectra should be used insted of the line plots
 only for univariate spectra.

 **See**: [`plot`](@ref)
 """
bar( S      :: Spectra;
    maxf    :: Union{Real, Int} = 0,
    space   :: Int              = 4,
    ylabel  :: String = "Power (\\muV²)",
   args...) =
  bar(S.y[1:(maxf==0 ? f2b(S.sr/2, S.sr, S.wl; DC=S.DC) : maxf), :];
      _spectraplotargs(S;
                   maxf   = maxf,
                   space  = space,
                   ylabel = ylabel)...,
      args...)

# internal function:
# generate appropriate arguments for plotting spectra using the PLOTS package,
# given a Spectra object and optional keyword argumants max frequency `maxf`,
# `space` and `ylabel`. This function is called by the plot and bar
# functions declared in this unit.
# It needs unpacking (...), for example:
# plot(S, _spectraplotargs(128, 512, 32; space=4)...)
function _spectraplotargs( S  :: Spectra;
                  maxf   :: Union{Real, Int} = 0,
                  space  :: Int              = 4,
                  ylabel :: String = "Power (\\muV²)")
  if maxf==0 maxf=S.sr/2 end
  b=f2b(maxf, S.sr, S.wl)
  xticstring=["$(space*i)" for i=1:maxf÷space]
  pushfirst!(xticstring, "$(fres(S.sr, S.wl))")
  xticrange=[1:S.wl/S.sr*space:b+1;]
  for i=2:length(xticrange) xticrange[i]=xticrange[i]-1 end
  return (labels        = "",
          dpi           = 300, #linecolor=:grey,
          xaxis         = ("Frequency (Hz)", font(11, "Courier New")),
          xticks        = (xticrange, xticstring),
          yaxis         = (ylabel, font(11, "Courier New")),
          left_margin   = 2mm,
          bottom_margin = 2mm)
end

"""
```
function heatmap(Y    :: TFAnalyticSignal,
                 func :: Function;
            c :: Symbol = :default,
            args...)
```

Generate heatmaps from a [TFAnalyticSignal](@ref) object.
Since the analytic signal (AS) is complex, a function `func`
must be specified and must transform the AS data
matrix into a real matrix. Typical functions are:

- `real`: heatmap of the real part of the AS
- `imag`: heatmap of the imaginary part of the AS
- `amplitude`: heatmap of the AS [`amplitude`](@ref)
- `phase`: heatmap of the AS [`phase`](@ref)

Note:

Appropriate colors for plotting analytic signal and phase are `:pu_or`, `:bluesreds`.
Appropriate colors for plotting amplitude are `:amp`, `:fire`, `:dimgray`, `:gwv`.

**Examples**:
```
using FourierAnalysis, Plots
sr, t, f, a = 128, 128, 8, 2
# create a sinusoidal at 8Hz
wave=sinusoidal(a, f, sr, t*4, 0)
# create a modulating sinusoidal at 0.25Hz
modu=sinusoidal(a, 0.25, sr, t*4, 0)
# create modulated wave and superimpose white noise
v=wave .* modu + randn(t*4)
plot(v)

# compute analytic signal (AS) of vector v
Y=TFanalyticsignal(v, sr, t; fmax=32)

# plot the real part of the AS
heatmap(Y, real)

# ...the imaginary part of the AS
heatmap(Y, imag)

# ...the amplitude of the AS
heatmap(Y, amplitude)

# ...the amplitude of the AS smoothed in the freq. dim.
heatmap(smooth(hannSmoother, noSmoother, Y), amplitude)

# ...the amplitude of the AS smoothed in freq. and time
heatmap(smooth(hannSmoother, hannSmoother, Y), amplitude)

# ...the phase of the AS and use custom colors
heatmap(Y, phase; c=:pu_or)
```
"""
function heatmap(Y    :: TFAnalyticSignal,
                 func :: Function;
            c :: Symbol = :amp,
            args...)
   if c==:default
      (func==real || func==imag || func==phase) ? c=:bluesreds :
      c=:amp
   end
   heatmap(_tfaxis(Y)..., func(Y.y), yflip = false, c=c, args...)
end

"""
```
function heatmap(Y :: TFAmplitude;
        c :: Symbol = :default,
        args...) =
```

Generate heatmaps of the data in a [TFAmplitude](@ref) object.

Note: appropriate colors for plotting amplitude are `:amp`, `:fire`, `:dimgray`, `:gwv`.

**Examples**:
```
using FourierAnalysis, Plots
sr, t, f, a = 128, 128, 8, 2
# create a sinusoidal at 8Hz
wave=sinusoidal(a, f, sr, t*4, 0)
# create a modulating sinusoidal at 0.25Hz
modu=sinusoidal(a, 0.25, sr, t*4, 0)
# create modulated wave and superimpose white noise
v=wave .* modu + randn(t*4)
plot(v)

# compute amplitude of vector v
A=TFamplitude(v, sr, t; fmax=32)

heatmap(A)
# or, quicly: heatmap(TFamplitude(v, sr, t; fmax=32))

heatmap(A; c=:fire)
```
"""
heatmap(Y :: TFAmplitude;
  c :: Symbol = :amp,
  args...) =
heatmap(_tfaxis(Y)..., Y.y, yflip = false, c=c, args...)

"""
```
function heatmap(Y :: TFPhase;
              c :: Symbol = :default,
              args...)
```
Generate heatmaps of the data in a [TFPhase](@ref) object.
Default colors for the heatmap are `c=:amp` if the phase object is unwrapped,
`c=:bluesreds` otherwise.

Note:

Appropriate colors for plotting phase: `:pu_or`, `:bluesreds`
Appropriate colors for plotting unwrapped phase are `:amp`, `:fire`, `:dimgray`, `:gwv`.

**See also**: [`isUnwrapped`](@ref)

**Examples**:
```
using FourierAnalysis, Plots
sr, t, f, a = 128, 128, 8, 2
# create a sinusoidal at 8Hz
wave=sinusoidal(a, f, sr, t*4, 0)
# create a modulating sinusoidal at 0.25Hz
modu=sinusoidal(a, 0.25, sr, t*4, 0)
# create modulated wave and superimpose white noise
v=wave .* modu + randn(t*4)
plot(v)

# heatmap of the phase of vector v
heatmap(TFphase(v, sr, t; fmax=32))

# heatmap of the unwrapped phase of vector v
heatmap(unwrapPhase(TFphase(v, sr, t; fmax=32)); c=:fire)
```
"""
function heatmap(Y :: TFPhase;
              c :: Symbol = :default,
              args...)
  c==:default ? (Y.unwrapped ? c=:amp : c=:bluesreds) : nothing
  heatmap(_tfaxis(Y)..., Y.y, yflip = false, c=c, args...)
end

#Generate labels for the axes to be used in heatmap plots of [TFobjects](@ref).
#On the y-axis the labels correspon to the center frequencies of the band-pass
#filter bank used to obtain the time-frequency representation, that is,
#to the .flabels field of the object.

## Examples:
#    using FourierAnalysis, Plots
#    sr, t, f, a = 128, 128, 8, 2
    # create a sinusoidal at 8Hz
#    wave=sinusoidal(a, f, sr, t*4, 0)
    # create a modulating sinusoidal at 0.25Hz
#    modu=sinusoidal(a, 0.25, sr, t*4, 0)
    # create modulated wave and superimpose white noise
#    v=wave .* modu + randn(t*4)

    # compute analytic signal of vector v
#    Y=TFanalyticsignal(v, sr, t; fmax=32)
    # plot amplitude
#    heatmap(_tfaxis(Y)..., amplitude(Y.y), yflip=false, c=:amp)
    # or create amplitude object
#    A=TFamplitude(Y)
    # and plot it
#    heatmap(_tfaxis(A)..., A.y, c=::pu_or)

    # plot phase
#    heatmap(_tfaxis(Y)..., phase(Y.y), c=:bluesreds)
    # plot phase in [0, 2π]
#    heatmap(_tfaxis(Y)..., phase(Y.y, func=x->x+π), c=:amp)
    # unwrapped phase
#    heatmap(_tfaxis(Y)..., phase(Y; unwrapped=true), c=:amp)
    # or create phase object
#    ϴ=TFphase(Y)
    # and plot it
#    heatmap(_tfaxis(ϴ)..., ϴ.y, c=:bluesreds)

# Be careful with phase! Smoothing phase is unappropriate
# since the phase is a discontinous function, however it makes sense to smooth
# phase if the phase is unwrapped.
# Now, this will 'illegally' plot the smoothed phase nonetheless:
## heatmap(tfaxis(Z)..., phase(Z.y), c=:bkr) # bluesreds
# The line here below would not create anything and return an error
# ϴ=smooth(blackmanSmoother, noSmoother, TFphase(Y))
# The line here below instead will smooth the unwrapped phase (in frequencies)
## ϴ=smooth(blackmanSmoother, noSmoother, TFphase(Y, unwrapped=true))
# plot it
## heatmap(tfaxis(ϴ)..., ϴ.y, c=:amp) # bluesreds

#Note:
#appropriate colors for plotting analytic signal and phase: :pu_or, :bluesreds
#appropriate colors for plotting amplitude: :amp, :fire, :dimgray, ;gwv
function _tfaxis(Y::TFobjects)
    function formatf(r::Real)
        s="$r"
        ss=split(s, ".")
        length(ss)>1 && parse(Float64, ss[2])==0.0 ? (return ss[1]) : (return s)
    end
    xs = [string(i) for i = 1:size(Y.y, 2)]
    ys = Vector{String}([formatf(f) for f in Y.flabels])
    return xs, ys
end


text:
"""
    NoneSelector()

Assign all dates to the validation set, select no holdout dates.
"""
struct NoneSelector <: DateSelector end

Iterators.partition(dates::StepRange{Date, Day}, ::NoneSelector) = _getdatesets(dates, Date[])


text:
@doc raw"""
`HDAE`: Hamiltonian Differential Algebraic Equation *EXPERIMENTAL*

Defines a Hamiltonian differential algebraic initial value problem, that is
a canonical Hamiltonian system of equations subject to Dirac constraints,
```math
\begin{align*}
\dot{q} (t) &= v_1(t, q(t), p(t)) + v_2(t, q(t), p(t), \lambda(t)) + v_3(t, q(t), p(t), \lambda(t), \gamma(t)) , & q(t_{0}) &= q_{0} , \\
\dot{p} (t) &= f_1(t, q(t), p(t)) + f_2(t, q(t), p(t), \lambda(t)) + f_3(t, q(t), p(t), \lambda(t), \gamma(t)) , & p(t_{0}) &= p_{0} , \\
0 &= \phi (t, q(t), p(t)) , \\
0 &= \psi (t, q(t), p(t), \dot{q}(t), \dot{p}(t)) ,
\end{align*}
```
with vector fields ``v_i`` and ``f_i`` for ``i = 1 ... 3``,
primary constraint ``\phi(q,p)=0`` and secondary constraint ``\psi(q,p,\lambda)=0``,
initial conditions ``(q_{0}, p_{0})``, the dynamical variables ``(q,p)``
taking values in ``\mathbb{R}^{d} \times \mathbb{R}^{d}`` and
the algebraic variables ``(\lambda, \gamma)`` taking values in
``\mathbb{R}^{n} \times \mathbb{R}^{d}``.

### Fields

* `d`: dimension of dynamical variables ``q`` and ``p`` as well as the vector fields ``v`` and ``f``
* `m`: dimension of algebraic variables ``\lambda`` and ``\gamma`` and the constraints ``\phi`` and ``\psi``
* `n`: number of initial conditions
* `v`: tuple of functions computing the vector fields ``v_i``, ``i = 1 ... 3``
* `f`: tuple of functions computing the vector fields ``f_i``, ``i = 1 ... 3``
* `ϕ`: primary constraints
* `ψ`: secondary constraints
* `t₀`: initial time
* `q₀`: initial condition for dynamical variable ``q``
* `p₀`: initial condition for dynamical variable ``p``

"""
struct HDAE{dType <: Number, tType <: Number, vType <: Tuple, fType <: Tuple, ϕType <: Function, ψType <: Function, N} <: AbstractEquationPDAE{dType, tType}
    d::Int
    m::Int
    n::Int
    v::vType
    f::fType
    ϕ::ϕType
    ψ::ψType
    t₀::tType
    q₀::Array{dType, N}
    p₀::Array{dType, N}

    function HDAE{dType,tType,vType,fType,ϕType,ψType,N}(d, m, n, v, f, ϕ, ψ, t₀, q₀, p₀) where {dType <: Number, tType <: Number, vType <: Tuple, fType <: Tuple, ϕType <: Function, ψType <: Function, N}
        @assert d == size(q₀,1) == size(p₀,1)
        @assert n == size(q₀,2) == size(p₀,2)
        @assert 2n ≥ m

        @assert dType == eltype(q₀)
        @assert dType == eltype(p₀)

        @assert ndims(q₀) == ndims(p₀) == N ∈ (1,2)

        new(d, m, n, v, f, ϕ, ψ, t₀, q₀, p₀)
    end
end

function HDAE(v::VT, f::FT, ϕ::ΦT, ψ::ΨT, m::Int, t₀::TT, q₀::AbstractArray{DT}, p₀::AbstractArray{DT}) where {DT,TT,VT,FT,ΦT,ΨT}
    @assert size(q₀) == size(p₀)
    HDAE{DT, TT, VT, FT, ΦT, ΨT, ndims(q₀)}(size(q₀, 1), m, size(q₀, 2), v, f, ϕ, ψ, t₀, q₀, p₀)
end

function HDAE(v::Function, f::Function, ϕ::Function, ψ::Function, t₀, q₀, p₀)
    HDAE(v, f, ϕ, ψ, size(q₀,1), t₀, q₀, p₀)
end

function HDAE(v::Function, f::Function, ϕ::Function, ψ::Function, q₀, p₀)
    HDAE(v, f, ϕ, ψ, size(q₀,1), zero(Float64), q₀, p₀)
end

Base.hash(dae::HDAE, h::UInt) = hash(dae.d, hash(dae.m, hash(dae.n, hash(dae.v, hash(dae.f, hash(dae.ϕ, hash(dae.ψ, hash(dae.t₀, hash(dae.q₀, hash(dae.p₀, h))))))))))
Base.:(==)(dae1::HDAE, dae2::HDAE) = (
                                dae1.d == dae2.d
                             && dae1.m == dae2.m
                             && dae1.n == dae2.n
                             && dae1.v == dae2.v
                             && dae1.f == dae2.f
                             && dae1.ϕ == dae2.ϕ
                             && dae1.ψ == dae2.ψ
                             && dae1.t₀ == dae2.t₀
                             && dae1.q₀ == dae2.q₀
                             && dae1.p₀ == dae2.p₀)

@inline Base.ndims(dae::HDAE) = ode.d

@inline periodicity(equation::HDAE) = equation.periodicity


text:
function write_model_reactions_to_file(data_dictionary::Dict{String,Any}, path_to_reaction_file::String)

    reaction_list = data_dictionary["list_of_chemical_reaction_strings"]
    writedlm(path_to_reaction_file,reaction_list)

end


text:
# https://github.com/stevengj/PyPlot.jl

is_marker_supported(::PyPlotBackend, shape::Shape) = true

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

# problem: https://github.com/tbreloff/Plots.jl/issues/308
# solution: hack from @stevengj: https://github.com/stevengj/PyPlot.jl/pull/223#issuecomment-229747768
otherdisplays = splice!(Base.Multimedia.displays, 2:length(Base.Multimedia.displays))
append!(Base.Multimedia.displays, otherdisplays)
pycolors = PyPlot.pyimport("matplotlib.colors")
pypath = PyPlot.pyimport("matplotlib.path")
mplot3d = PyPlot.pyimport("mpl_toolkits.mplot3d")
axes_grid1 = PyPlot.pyimport("mpl_toolkits.axes_grid1")
pypatches = PyPlot.pyimport("matplotlib.patches")
pyfont = PyPlot.pyimport("matplotlib.font_manager")
pyticker = PyPlot.pyimport("matplotlib.ticker")
pycmap = PyPlot.pyimport("matplotlib.cm")
pynp = PyPlot.pyimport("numpy")
pynp."seterr"(invalid="ignore")
pytransforms = PyPlot.pyimport("matplotlib.transforms")
pycollections = PyPlot.pyimport("matplotlib.collections")
pyart3d = PyPlot.art3D
pyrcparams = PyPlot.PyDict(PyPlot.matplotlib."rcParams")

# "support" matplotlib v1.5
set_facecolor_sym = if PyPlot.version < v"2"
    @warn("You are using Matplotlib $(PyPlot.version), which is no longer officialy supported by the Plots community. To ensure smooth Plots.jl integration update your Matplotlib library to a version >= 2.0.0")
    :set_axis_bgcolor
else
    :set_facecolor
end

# PyCall API changes in v1.90.0
if !isdefined(PyPlot.PyCall, :_setproperty!)
    @warn "Plots no longer supports PyCall < 1.90.0 and PyPlot < 2.8.0. Either update PyCall and PyPlot or pin Plots to a version <= 0.23.2."
end

# # convert colorant to 4-tuple RGBA
# py_color(c::Colorant, α=nothing) = map(f->float(f(convertColor(c,α))), (red, green, blue, alpha))
# py_color(cvec::ColorVector, α=nothing) = map(py_color, convertColor(cvec, α).v)
# py_color(grad::ColorGradient, α=nothing) = map(c -> py_color(c, α), grad.colors)
# py_color(scheme::ColorScheme, α=nothing) = py_color(convertColor(getColor(scheme), α))
# py_color(vec::AVec, α=nothing) = map(c->py_color(c,α), vec)
# py_color(c, α=nothing) = py_color(convertColor(c, α))

# function py_colormap(c::ColorGradient, α=nothing)
#     pyvals = [(v, py_color(getColorZ(c, v), α)) for v in c.values]
#     pycolors["LinearSegmentedColormap"][:from_list]("tmp", pyvals)
# end

# # convert vectors and ColorVectors to standard ColorGradients
# # TODO: move this logic to colors.jl and keep a barebones wrapper for pyplot
# py_colormap(cv::ColorVector, α=nothing) = py_colormap(ColorGradient(cv.v), α)
# py_colormap(v::AVec, α=nothing) = py_colormap(ColorGradient(v), α)

# # anything else just gets a bluesred gradient
# py_colormap(c, α=nothing) = py_colormap(default_gradient(), α)

py_color(s) = py_color(parse(Colorant, string(s)))
py_color(c::Colorant) = (red(c), green(c), blue(c), alpha(c))
py_color(cs::AVec) = map(py_color, cs)
py_color(grad::PlotUtils.AbstractColorList) = py_color(color_list(grad))
py_color(c::Colorant, α) = py_color(plot_color(c, α))

function py_colormap(cg::ColorGradient)
    pyvals = collect(zip(cg.values, py_color(PlotUtils.color_list(cg))))
    cm = pycolors."LinearSegmentedColormap"."from_list"("tmp", pyvals)
    cm."set_bad"(color=(0,0,0,0.0), alpha=0.0)
    cm
end
function py_colormap(cg::PlotUtils.CategoricalColorGradient)
    r = range(0, stop = 1, length = 256)
    pyvals = collect(zip(r, py_color(cg[r])))
    cm = pycolors."LinearSegmentedColormap"."from_list"("tmp", pyvals)
    cm."set_bad"(color=(0,0,0,0.0), alpha=0.0)
    cm
end
py_colormap(c) = py_colormap(_as_gradient(c))

function py_shading(c, z)
    cmap = py_colormap(c)
    ls = pycolors."LightSource"(270,45)
    ls."shade"(z, cmap, vert_exag=0.1, blend_mode="soft")
end

# get the style (solid, dashed, etc)
function py_linestyle(seriestype::Symbol, linestyle::Symbol)
    seriestype == :none && return " "
    linestyle == :solid && return "-"
    linestyle == :dash && return "--"
    linestyle == :dot && return ":"
    linestyle == :dashdot && return "-."
    @warn("Unknown linestyle $linestyle")
    return "-"
end

function py_marker(marker::Shape)
    x, y = coords(marker)
    n = length(x)
    mat = zeros(n+1,2)
    for i=1:n
        mat[i,1] = x[i]
        mat[i,2] = y[i]
    end
    mat[n+1,:] = mat[1,:]
    pypath."Path"(mat)
end

# get the marker shape
function py_marker(marker::Symbol)
    marker == :none && return " "
    marker == :circle && return "o"
    marker == :rect && return "s"
    marker == :diamond && return "D"
    marker == :utriangle && return "^"
    marker == :dtriangle && return "v"
    marker == :+ && return "+"
    marker == :x && return "x"
    marker == :star5 && return "*"
    marker == :pentagon && return "p"
    marker == :hexagon && return "h"
    marker == :octagon && return "8"
    marker == :pixel && return ","
    marker == :hline && return "_"
    marker == :vline && return "|"
    haskey(_shapes, marker) && return py_marker(_shapes[marker])

    @warn("Unknown marker $marker")
    return "o"
end

# py_marker(markers::AVec) = map(py_marker, markers)
function py_marker(markers::AVec)
    @warn("Vectors of markers are currently unsupported in PyPlot: $markers")
    py_marker(markers[1])
end

# pass through
function py_marker(marker::AbstractString)
    @assert length(marker) == 1
    marker
end

function py_stepstyle(seriestype::Symbol)
    seriestype == :steppost && return "steps-post"
    seriestype == :steppre && return "steps-pre"
    return "default"
end

function py_fillstepstyle(seriestype::Symbol)
    seriestype == :steppost && return "post"
    seriestype == :steppre && return "pre"
    return nothing
end

# # untested... return a FontProperties object from a Plots.Font
# function py_font(font::Font)
#     pyfont["FontProperties"](
#         family = font.family,
#         size = font.size
#     )
# end

function get_locator_and_formatter(vals::AVec)
    pyticker."FixedLocator"(eachindex(vals)), pyticker."FixedFormatter"(vals)
end

function add_pyfixedformatter(cbar, vals::AVec)
    cbar[:locator], cbar[:formatter] = get_locator_and_formatter(vals)
    cbar[:update_ticks]()
end

function labelfunc(scale::Symbol, backend::PyPlotBackend)
    PyPlot.LaTeXStrings.latexstring ∘ labelfunc_tex(scale)
end

function py_mask_nans(z)
    # pynp["ma"][:masked_invalid](z)))
    PyPlot.PyCall.pycall(pynp."ma"."masked_invalid", Any, z)
    # pynp["ma"][:masked_where](pynp["isnan"](z),z)
end

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

function fix_xy_lengths!(plt::Plot{PyPlotBackend}, series::Series)
    if series[:x] !== nothing
        x, y = series[:x], series[:y]
        nx, ny = length(x), length(y)
        if !isa(get(series.plotattributes, :z, nothing), Surface) && nx != ny
            if nx < ny
                series[:x] = Float64[x[mod1(i,nx)] for i=1:ny]
            else
                series[:y] = Float64[y[mod1(i,ny)] for i=1:nx]
            end
        end
    end
end

py_linecolormap(series::Series)       = py_colormap(series[:linecolor])
py_markercolormap(series::Series)     = py_colormap(series[:markercolor])
py_fillcolormap(series::Series)       = py_colormap(series[:fillcolor])

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

# TODO: these can probably be removed eventually... right now they're just keeping things working before cleanup

# getAxis(sp::Subplot) = sp.o

# function getAxis(plt::Plot{PyPlotBackend}, series::Series)
#     sp = get_subplot(plt, get(series.plotattributes, :subplot, 1))
#     getAxis(sp)
# end

# getfig(o) = o

# ---------------------------------------------------------------------------
# Figure utils -- F*** matplotlib for making me work so hard to figure this crap out

# the drawing surface
py_canvas(fig) = fig."canvas"

# the object controlling draw commands
py_renderer(fig) = py_canvas(fig)."get_renderer"()

# draw commands... paint the screen (probably updating internals too)
py_drawfig(fig) = fig."draw"(py_renderer(fig))
# py_drawax(ax) = ax[:draw](py_renderer(ax[:get_figure]()))

# get a vector [left, right, bottom, top] in PyPlot coords (origin is bottom-left!)
py_extents(obj) = obj."get_window_extent"()."get_points"()

# compute a bounding box (with origin top-left), however pyplot gives coords with origin bottom-left
function py_bbox(obj)
    fl, fr, fb, ft = py_extents(obj."get_figure"())
    l, r, b, t = py_extents(obj)
    BoundingBox(l*px, (ft-t)*px, (r-l)*px, (t-b)*px)
end

py_bbox(::Nothing) = BoundingBox(0mm, 0mm)

# get the bounding box of the union of the objects
function py_bbox(v::AVec)
    bbox_union = defaultbox
    for obj in v
        bbox_union = bbox_union + py_bbox(obj)
    end
    bbox_union
end

# bounding box: union of axis tick labels
function py_bbox_ticks(ax, letter)
    labels = getproperty(ax, Symbol("get_"*letter*"ticklabels"))()
    py_bbox(labels)
end

# bounding box: axis guide
function py_bbox_axislabel(ax, letter)
    pyaxis_label = getproperty(ax, Symbol("get_"*letter*"axis"))().label
    py_bbox(pyaxis_label)
end

# bounding box: union of axis ticks and guide
function py_bbox_axis(ax, letter)
    ticks = py_bbox_ticks(ax, letter)
    labels = py_bbox_axislabel(ax, letter)
    ticks + labels
end

# bounding box: axis title
function py_bbox_title(ax)
    bb = defaultbox
    for s in (:title, :_left_title, :_right_title)
        bb = bb + py_bbox(getproperty(ax, s))
    end
    bb
end

# bounding box: legend
py_bbox_legend(ax) = py_bbox(ax."get_legend"())

function py_thickness_scale(plt::Plot{PyPlotBackend}, ptsz)
    ptsz * plt[:thickness_scaling]
end

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

# Create the window/figure for this backend.
function _create_backend_figure(plt::Plot{PyPlotBackend})
    w,h = map(px2inch, Tuple(s * plt[:dpi] / Plots.DPI for s in plt[:size]))

    # # reuse the current figure?
    fig = if plt[:overwrite_figure]
        PyPlot.gcf()
    else
        fig = PyPlot.figure()
        # finalizer(fig, close)
        fig
    end

    # clear the figure
    # PyPlot.clf()
    fig
end

# Set up the subplot within the backend object.
# function _initialize_subplot(plt::Plot{PyPlotBackend}, sp::Subplot{PyPlotBackend})

function py_init_subplot(plt::Plot{PyPlotBackend}, sp::Subplot{PyPlotBackend})
    fig = plt.o
    proj = sp[:projection]
    proj = (proj in (nothing,:none) ? nothing : string(proj))

    # add a new axis, and force it to create a new one by setting a distinct label
    ax = fig."add_axes"(
        [0,0,1,1],
        label = string(gensym()),
        projection = proj
    )
    sp.o = ax
end

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

# function _series_added(pkg::PyPlotBackend, plt::Plot, plotattributes::KW)
# TODO: change this to accept Subplot??
# function _series_added(plt::Plot{PyPlotBackend}, series::Series)

function py_add_series(plt::Plot{PyPlotBackend}, series::Series)
    # plotattributes = series.plotattributes
    st = series[:seriestype]
    sp = series[:subplot]
    ax = sp.o

    # PyPlot doesn't handle mismatched x/y
    fix_xy_lengths!(plt, series)

    # ax = getAxis(plt, series)
    x, y, z = series[:x], series[:y], series[:z]
    if st == :straightline
        x, y = straightline_data(series)
    elseif st == :shape
        x, y = shape_data(series)
    end

    if ispolar(series)
        # make negative radii positive and flip the angle
        # (PyPlot ignores negative radii)
        for i in eachindex(y)
            if y[i] < 0
                y[i] = -y[i]
                x[i] -= π
            end
        end
    end

    xyargs = (st in _3dTypes ? (x,y,z) : (x,y))

    # handle zcolor and get c/cmap
    needs_colorbar = hascolorbar(sp)
    vmin, vmax = clims = get_clims(sp, series)

    # Dict to store extra kwargs
    if st == :wireframe
        extrakw = KW()          # vmin, vmax cause an error for wireframe plot
    else
        extrakw = KW(:vmin => vmin, :vmax => vmax)
    end

    # holds references to any python object representing the matplotlib series
    handles = []
    discrete_colorbar_values = nothing

    # pass in an integer value as an arg, but a levels list as a keyword arg
    levels = series[:levels]
    levelargs = if isscalar(levels)
        (levels)
    elseif isvector(levels)
        extrakw[:levels] = levels
        ()
    else
        error("Only numbers and vectors are supported with levels keyword")
    end

    # add custom frame shapes to markershape?
    series_annotations_shapes!(series, :xy)

    # for each plotting command, optionally build and add a series handle to the list

    # line plot
    if st in (:path, :path3d, :steppre, :steppost, :straightline)
        if maximum(series[:linewidth]) > 0
            # TODO: check LineCollection alternative for speed
            # if length(segments) > 1 && (any(typeof(series[attr]) <: AbstractVector for attr in (:fillcolor, :fillalpha)) || series[:fill_z] !== nothing) && !(typeof(series[:linestyle]) <: AbstractVector)
            #     # multicolored line segments
            #     n = length(segments)
            #     # segments = Array(Any,n)
            #     segments = []
            #     kw = KW(
            #         :label => series[:label],
            #         :zorder => plt.n,
            #         :cmap => py_linecolormap(series),
            #         :linewidths => py_thickness_scale(plt, get_linewidth.(series, 1:n)),
            #         :linestyle => py_linestyle(st, get_linestyle.(series)),
            #         :norm => pycolors["Normalize"](; extrakw...)
            #     )
            #     lz = _cycle(series[:line_z], 1:n)
            #     handle = if RecipesPipeline.is3d(st)
            #         line_segments = [[(x[j], y[j], z[j]) for j in rng] for rng in segments]
            #         lc = pyart3d["Line3DCollection"](line_segments; kw...)
            #         lc[:set_array](lz)
            #         ax[:add_collection3d](lc, zs=z) #, zdir='y')
            #         lc
            #     else
            #         line_segments = [[(x[j], y[j]) for j in rng] for rng in segments]
            #         lc = pycollections["LineCollection"](line_segments; kw...)
            #         lc[:set_array](lz)
            #         ax[:add_collection](lc)
            #         lc
            #     end
            #     push!(handles, handle)
            # else
            for (i, rng) in enumerate(iter_segments(series, st))
                handle = ax."plot"((arg[rng] for arg in xyargs)...;
                                   label = i == 1 ? series[:label] : "",
                                   zorder = series[:series_plotindex],
                                   color = py_color(single_color(get_linecolor(series, clims, i)), get_linealpha(series, i)),
                                   linewidth = py_thickness_scale(plt, get_linewidth(series, i)),
                                   linestyle = py_linestyle(st, get_linestyle(series, i)),
                                   solid_capstyle = "butt",
                                   dash_capstyle = "butt",
                                   drawstyle = py_stepstyle(st)
                                   )[1]
                push!(handles, handle)
            end
            # end

            a = series[:arrow]
            if a !== nothing && !RecipesPipeline.is3d(st)  # TODO: handle 3d later
                if typeof(a) != Arrow
                    @warn("Unexpected type for arrow: $(typeof(a))")
                else
                    arrowprops = KW(
                        :arrowstyle => "simple,head_length=$(a.headlength),head_width=$(a.headwidth)",
                        :shrinkA => 0,
                        :shrinkB => 0,
                        :edgecolor => py_color(get_linecolor(series)),
                        :facecolor => py_color(get_linecolor(series)),
                        :linewidth => py_thickness_scale(plt, get_linewidth(series)),
                        :linestyle => py_linestyle(st, get_linestyle(series)),
                    )
                    add_arrows(x, y) do xyprev, xy
                        ax."annotate"("",
                            xytext = (0.001xyprev[1] + 0.999xy[1], 0.001xyprev[2] + 0.999xy[2]),
                            xy = xy,
                            arrowprops = arrowprops,
                            zorder = 999
                        )
                    end
                end
            end
        end
    end

    # add markers?
    if series[:markershape] != :none && st in (
        :path, :scatter, :path3d, :scatter3d, :steppre, :steppost, :bar
    )
        for (i, rng) in enumerate(iter_segments(series, :scatter))
            xyargs = if st == :bar && !isvertical(series)
                y[rng], x[rng]
            else
                x[rng], y[rng]
            end

            handle = ax."scatter"(xyargs...;
                label = series[:label],
                zorder = series[:series_plotindex] + 0.5,
                marker = py_marker(_cycle(series[:markershape], i)),
                s =  py_thickness_scale(plt, _cycle(series[:markersize], i)).^ 2,
                facecolors = py_color(get_markercolor(series, i), get_markeralpha(series, i)),
                edgecolors = py_color(get_markerstrokecolor(series, i), get_markerstrokealpha(series, i)),
                linewidths = py_thickness_scale(plt, get_markerstrokewidth(series, i)),
                extrakw...
            )
            push!(handles, handle)
        end
    end

    if st == :hexbin
        handle = ax."hexbin"(x, y,
            label = series[:label],
            C = series[:weights],
            gridsize = series[:bins]==:auto ? 100 : series[:bins],  # 100 is the default value
            linewidths = py_thickness_scale(plt, series[:linewidth]),
            edgecolors = py_color(get_linecolor(series)),
            alpha = series[:fillalpha],
            cmap = py_fillcolormap(series),  # applies to the pcolorfast object
            zorder = series[:series_plotindex],
            # extrakw...  # We are not supporting clims for hexbin as calculation of bins is not trivial
        )
        push!(handles, handle)
    end

    if st in (:contour, :contour3d)
        z = transpose_z(series, z.surf)
	if typeof(x)<:Plots.Surface
            x = Plots.transpose_z(series, x.surf)
        end
        if typeof(y)<:Plots.Surface
            y = Plots.transpose_z(series, y.surf)
        end

        if st == :contour3d
            extrakw[:extend3d] = true
        end

        if typeof(series[:linecolor]) <: AbstractArray
            extrakw[:colors] = py_color.(series[:linecolor])
        else
            extrakw[:cmap] = py_linecolormap(series)
        end

        # contour lines
        handle = ax."contour"(x, y, z, levelargs...;
            label = series[:label],
            zorder = series[:series_plotindex],
            linewidths = py_thickness_scale(plt, series[:linewidth]),
            linestyles = py_linestyle(st, series[:linestyle]),
            extrakw...
        )
        if series[:contour_labels] == true
            ax."clabel"(handle, handle.levels)
        end
        push!(handles, handle)

        # contour fills
        if series[:fillrange] !== nothing
            handle = ax."contourf"(x, y, z, levelargs...;
                label = series[:label],
                zorder = series[:series_plotindex] + 0.5,
                alpha = series[:fillalpha],
                extrakw...
            )
            push!(handles, handle)
        end
    end

    if st in (:surface, :wireframe)
        if typeof(z) <: AbstractMatrix || typeof(z) <: Surface
            x, y, z = map(Array, (x,y,z))
            if !ismatrix(x) || !ismatrix(y)
                x = repeat(x', length(y), 1)
                y = repeat(y, 1, length(series[:x]))
            end
            z = transpose_z(series, z)
            if st == :surface
                if series[:fill_z] !== nothing
                    # the surface colors are different than z-value
                    extrakw[:facecolors] = py_shading(series[:fillcolor], transpose_z(series, series[:fill_z].surf))
                    extrakw[:shade] = false
                else
                    extrakw[:cmap] = py_fillcolormap(series)
                end
            end
            handle = getproperty(ax, st == :surface ? :plot_surface : :plot_wireframe)(x, y, z;
                label = series[:label],
                zorder = series[:series_plotindex],
                rstride = series[:stride][1],
                cstride = series[:stride][2],
                linewidth = py_thickness_scale(plt, series[:linewidth]),
                edgecolor = py_color(get_linecolor(series)),
                extrakw...
            )
            push!(handles, handle)

            # contours on the axis planes
            if series[:contours]
                for (zdir,mat) in (("x",x), ("y",y), ("z",z))
                    offset = (zdir == "y" ? ignorenan_maximum : ignorenan_minimum)(mat)
                    handle = ax."contourf"(x, y, z, levelargs...;
                        zdir = zdir,
                        cmap = py_fillcolormap(series),
                        offset = (zdir == "y" ? ignorenan_maximum : ignorenan_minimum)(mat)  # where to draw the contour plane
                    )
                    push!(handles, handle)
                end
            end

        elseif typeof(z) <: AbstractVector
            # tri-surface plot (http://matplotlib.org/mpl_toolkits/mplot3d/tutorial.html#tri-surface-plots)
            handle = ax."plot_trisurf"(x, y, z;
                label = series[:label],
                zorder = series[:series_plotindex],
                cmap = py_fillcolormap(series),
                linewidth = py_thickness_scale(plt, series[:linewidth]),
                edgecolor = py_color(get_linecolor(series)),
                extrakw...
            )
            push!(handles, handle)
        else
            error("Unsupported z type $(typeof(z)) for seriestype=$st")
        end
    end

    if st == :image
        xmin, xmax = ignorenan_extrema(series[:x])
        ymin, ymax = ignorenan_extrema(series[:y])
        dx = (xmax - xmin) / (length(series[:x]) - 1) / 2
        dy = (ymax - ymin) / (length(series[:y]) - 1) / 2
        img = Array(transpose_z(series, z.surf))
        z = if eltype(img) <: Colors.AbstractGray
            float(img)
        elseif eltype(img) <: Colorant
            map(c -> Float64[red(c),green(c),blue(c),alpha(c)], img)
        else
            z  # hopefully it's in a data format that will "just work" with imshow
        end
        handle = ax."imshow"(z;
            zorder = series[:series_plotindex],
            cmap = py_colormap(cgrad(plot_color([:black, :white]))),
            vmin = 0.0,
            vmax = 1.0,
            extent = (xmin - dx, xmax + dx, ymax + dy, ymin - dy)
        )
        push!(handles, handle)

        # expand extrema... handle is AxesImage object
        xmin, xmax, ymax, ymin = handle."get_extent"()
        expand_extrema!(sp, xmin, xmax, ymin, ymax)
        # sp[:yaxis].series[:flip] = true
    end

    if st == :heatmap
        x, y, z = heatmap_edges(x, sp[:xaxis][:scale]), heatmap_edges(y, sp[:yaxis][:scale]), transpose_z(series, z.surf)

        expand_extrema!(sp[:xaxis], x)
        expand_extrema!(sp[:yaxis], y)
        dvals = sp[:zaxis][:discrete_values]
        if !isempty(dvals)
            discrete_colorbar_values = dvals
        end

        handle = ax."pcolormesh"(x, y, py_mask_nans(z);
            label = series[:label],
            zorder = series[:series_plotindex],
            cmap = py_fillcolormap(series),
            alpha = series[:fillalpha],
            # edgecolors = (series[:linewidth] > 0 ? py_linecolor(series) : "face"),
            extrakw...
        )
        push!(handles, handle)
    end

    if st == :shape
        handle = []
        for (i, rng) in enumerate(iter_segments(series))
            if length(rng) > 1
                path = pypath."Path"(hcat(x[rng], y[rng]))
                patches = pypatches."PathPatch"(
                    path;
                    label = series[:label],
                    zorder = series[:series_plotindex],
                    edgecolor = py_color(get_linecolor(series, clims, i), get_linealpha(series, i)),
                    facecolor = py_color(get_fillcolor(series, clims, i), get_fillalpha(series, i)),
                    linewidth = py_thickness_scale(plt, get_linewidth(series, i)),
                    linestyle = py_linestyle(st, get_linestyle(series, i)),
                    fill = true
                )
                push!(handle, ax."add_patch"(patches))
            end
        end
        push!(handles, handle)
    end

    series[:serieshandle] = handles

    # # smoothing
    # handleSmooth(plt, ax, series, series[:smooth])

    # handle area filling
    fillrange = series[:fillrange]
    if fillrange !== nothing && st != :contour
        for (i, rng) in enumerate(iter_segments(series))
            f, dim1, dim2 = if isvertical(series)
                :fill_between, x[rng], y[rng]
            else
                :fill_betweenx, y[rng], x[rng]
            end
            n = length(dim1)
            args = if typeof(fillrange) <: Union{Real, AVec}
                dim1, _cycle(fillrange, rng), dim2
            elseif is_2tuple(fillrange)
                dim1, _cycle(fillrange[1], rng), _cycle(fillrange[2], rng)
            end

            handle = getproperty(ax, f)(args..., trues(n), false, py_fillstepstyle(st);
                zorder = series[:series_plotindex],
                facecolor = py_color(get_fillcolor(series, clims, i), get_fillalpha(series, i)),
                linewidths = 0
            )
            push!(handles, handle)
        end
    end

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

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

function py_set_lims(ax, sp::Subplot, axis::Axis)
    letter = axis[:letter]
    lfrom, lto = axis_limits(sp, letter)
    getproperty(ax, Symbol("set_", letter, "lim"))(lfrom, lto)
end

function py_surround_latextext(latexstring, env)
    if latexstring[1] == '$' && latexstring[end] == '$'
        unenclosed = latexstring[2:end-1]
    else
        unenclosed = latexstring
    end
    PyPlot.LaTeXStrings.latexstring(env, "{", unenclosed, "}")
end

function py_set_ticks(ax, ticks, letter, env)
    ticks == :auto && return
    axis = getproperty(ax, Symbol(letter,"axis"))
    if ticks == :none || ticks === nothing || ticks == false
        kw = KW()
        for dir in (:top,:bottom,:left,:right)
            kw[dir] = kw[Symbol(:label,dir)] = false
        end
        axis."set_tick_params"(;which="both", kw...)
        return
    end

    ttype = ticksType(ticks)
    if ttype == :ticks
        axis."set_ticks"(ticks)
    elseif ttype == :ticks_and_labels
        axis."set_ticks"(ticks[1])

        if pyrcparams["text.usetex"]
            tick_labels = ticks[2]
        else
            tick_labels = [py_surround_latextext(ticklabel, env) for ticklabel in ticks[2]]
        end

        axis."set_ticklabels"(tick_labels)
    else
        error("Invalid input for $(letter)ticks: $ticks")
    end
end

function py_compute_axis_minval(sp::Subplot, axis::Axis)
    # compute the smallest absolute value for the log scale's linear threshold
    minval = 1.0
    sps = axis.sps
    for sp in sps
        for series in series_list(sp)
            v = series.plotattributes[axis[:letter]]
            if !isempty(v)
                minval = NaNMath.min(minval, ignorenan_minimum(abs.(v)))
            end
        end
    end

    # now if the axis limits go to a smaller abs value, use that instead
    vmin, vmax = axis_limits(sp, axis[:letter])
    minval = NaNMath.min(minval, abs(vmin), abs(vmax))

    minval
end

function py_set_scale(ax, sp::Subplot, axis::Axis)
    scale = axis[:scale]
    letter = axis[:letter]
    scale in supported_scales() || return @warn("Unhandled scale value in pyplot: $scale")
    func = getproperty(ax, Symbol("set_", letter, "scale"))
    kw = KW()
    arg = if scale == :identity
        "linear"
    else
        kw[Symbol(:base,letter)] = if scale == :ln
            ℯ
        elseif scale == :log2
            2
        elseif scale == :log10
            10
        end
        kw[Symbol(:linthresh,letter)] = NaNMath.max(1e-16, py_compute_axis_minval(sp, axis))
        "symlog"
    end
    func(arg; kw...)
end

function py_set_axis_colors(sp, ax, a::Axis)
    for (loc, spine) in ax.spines
        spine."set_color"(py_color(a[:foreground_color_border]))
    end
    axissym = Symbol(a[:letter], :axis)
    if PyPlot.PyCall.hasproperty(ax, axissym)
        tickcolor = sp[:framestyle] in (:zerolines, :grid) ? py_color(plot_color(a[:foreground_color_grid], a[:gridalpha])) : py_color(a[:foreground_color_axis])
        ax."tick_params"(axis=string(a[:letter]), which="both",
                         colors=tickcolor,
                         labelcolor=py_color(a[:tickfontcolor]))
        getproperty(ax, axissym).label.set_color(py_color(a[:guidefontcolor]))
    end
end

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

function _before_layout_calcs(plt::Plot{PyPlotBackend})
    # update the fig
    w, h = plt[:size]
    fig = plt.o
    fig."clear"()
    dpi = plt[:dpi]
    fig."set_size_inches"(w/DPI, h/DPI, forward = true)
    getproperty(fig, set_facecolor_sym)(py_color(plt[:background_color_outside]))
    fig."set_dpi"(plt[:dpi])

    # resize the window
    PyPlot.plt."get_current_fig_manager"().resize(w, h)

    # initialize subplots
    for sp in plt.subplots
        py_init_subplot(plt, sp)
    end

    # add the series
    for series in plt.series_list
        py_add_series(plt, series)
    end

    # update subplots
    for sp in plt.subplots
        ax = sp.o
        if ax === nothing
            continue
        end

        # add the annotations
        for ann in sp[:annotations]
            py_add_annotations(sp, locate_annotation(sp, ann...)...)
        end

        # title
        if sp[:title] != ""
            loc = lowercase(string(sp[:titlelocation]))
            func = if loc == "left"
                :_left_title
            elseif loc == "right"
                :_right_title
            else
                :title
            end
            getproperty(ax, func)."set_text"(sp[:title])
            getproperty(ax, func)."set_fontsize"(py_thickness_scale(plt, sp[:titlefontsize]))
            getproperty(ax, func)."set_family"(sp[:titlefontfamily])
            getproperty(ax, func)."set_color"(py_color(sp[:titlefontcolor]))
            # ax[:set_title](sp[:title], loc = loc)
        end

        # add the colorbar legend
        if hascolorbar(sp)
            # add keyword args for a discrete colorbar
            slist = series_list(sp)
            colorbar_series = slist[findfirst(hascolorbar.(slist))]
            handle = colorbar_series[:serieshandle][end]
            kw = KW()
            if !isempty(sp[:zaxis][:discrete_values]) && colorbar_series[:seriestype] == :heatmap
                locator, formatter = get_locator_and_formatter(sp[:zaxis][:discrete_values])
                # kw[:values] = eachindex(sp[:zaxis][:discrete_values])
                kw[:values] = sp[:zaxis][:continuous_values]
                kw[:ticks] = locator
                kw[:format] = formatter
                kw[:boundaries] = vcat(0, kw[:values] + 0.5)
            elseif any(colorbar_series[attr] !== nothing for attr in (:line_z, :fill_z, :marker_z))
                cmin, cmax = get_clims(sp)
                norm = pycolors."Normalize"(vmin = cmin, vmax = cmax)
                f = if colorbar_series[:line_z] !== nothing
                    py_linecolormap
                elseif colorbar_series[:fill_z] !== nothing
                    py_fillcolormap
                else
                    py_markercolormap
                end
                cmap = pycmap."ScalarMappable"(norm = norm, cmap = f(colorbar_series))
                cmap."set_array"([])
                handle = cmap
            end
            kw[:spacing] = "proportional"

            fig = plt.o

            if RecipesPipeline.is3d(sp)
                cbax = fig."add_axes"([0.9, 0.1, 0.03, 0.8])
                cb = fig."colorbar"(handle; cax=cbax, kw...)
            else
                # divider approach works only with 2d plots
                divider = axes_grid1.make_axes_locatable(ax)
                # width = axes_grid1.axes_size.AxesY(ax, aspect=1.0 / 3.5)
                # pad = axes_grid1.axes_size.Fraction(0.5, width)  # Colorbar is spaced 0.5 of its size away from the ax
                # cbax = divider.append_axes("right", size=width, pad=pad)   # This approach does not work well in subplots
                colorbar_position = "right"
                colorbar_pad = "2.5%"
                colorbar_orientation="vertical"

                if sp[:colorbar] == :left
                    colorbar_position = string(sp[:colorbar])
                    colorbar_pad = "5%"
                elseif sp[:colorbar] == :top
                    colorbar_position = string(sp[:colorbar])
                    colorbar_pad = "2.5%"
                    colorbar_orientation="horizontal"
                elseif sp[:colorbar] == :bottom
                    colorbar_position = string(sp[:colorbar])
                    colorbar_pad = "5%"
                    colorbar_orientation="horizontal"
                end

                cbax = divider.append_axes(colorbar_position, size="5%", pad=colorbar_pad)  # Reasonable value works most of the usecases
                cb = fig."colorbar"(handle; cax=cbax, orientation = colorbar_orientation, kw...)

                if sp[:colorbar] == :left
                    cbax.yaxis.set_ticks_position("left")
                elseif sp[:colorbar] == :top
                    cbax.xaxis.set_ticks_position("top")
                elseif sp[:colorbar] == :bottom
                    cbax.xaxis.set_ticks_position("bottom")
                end

            end

            cb."set_label"(sp[:colorbar_title],size=py_thickness_scale(plt, sp[:yaxis][:guidefontsize]),family=sp[:yaxis][:guidefontfamily], color = py_color(sp[:yaxis][:guidefontcolor]))

            # cb."formatter".set_useOffset(false)   # This for some reason does not work, must be a pyplot bug, instead this is a workaround:
            cb."formatter".set_powerlimits((-Inf, Inf))
            cb."update_ticks"()

            for lab in cb."ax"."yaxis"."get_ticklabels"()
                  lab."set_fontsize"(py_thickness_scale(plt, sp[:yaxis][:tickfontsize]))
                  lab."set_family"(sp[:yaxis][:tickfontfamily])
                  lab."set_color"(py_color(sp[:yaxis][:tickfontcolor]))
            end
            sp.attr[:cbar_handle] = cb
            sp.attr[:cbar_ax] = cbax
        end

        # framestyle
        if !ispolar(sp) && !RecipesPipeline.is3d(sp)
            ax.spines["left"]."set_linewidth"(py_thickness_scale(plt, 1))
            ax.spines["bottom"]."set_linewidth"(py_thickness_scale(plt, 1))
            if sp[:framestyle] == :semi
                intensity = 0.5
                ax.spines["right"]."set_alpha"(intensity)
                ax.spines["top"]."set_alpha"(intensity)
                ax.spines["right"]."set_linewidth"(py_thickness_scale(plt, intensity))
                ax.spines["top"]."set_linewidth"(py_thickness_scale(plt, intensity))
            elseif sp[:framestyle] in (:axes, :origin)
                ax.spines["right"]."set_visible"(false)
                ax.spines["top"]."set_visible"(false)
                if sp[:framestyle] == :origin
                    ax.spines["bottom"]."set_position"("zero")
                    ax.spines["left"]."set_position"("zero")
                end
            elseif sp[:framestyle] in (:grid, :none, :zerolines)
                for (loc, spine) in ax.spines
                    spine."set_visible"(false)
                end
                if sp[:framestyle] == :zerolines
                    ax."axhline"(y = 0, color = py_color(sp[:xaxis][:foreground_color_axis]), lw = py_thickness_scale(plt, 0.75))
                    ax."axvline"(x = 0, color = py_color(sp[:yaxis][:foreground_color_axis]), lw = py_thickness_scale(plt, 0.75))
                end
            end
        end

        # axis attributes
        for letter in (:x, :y, :z)
            axissym = Symbol(letter, :axis)
            PyPlot.PyCall.hasproperty(ax, axissym) || continue
            axis = sp[axissym]
            pyaxis = getproperty(ax, axissym)
            if axis[:mirror] && letter != :z
                pos = letter == :x ? "top" : "right"
                pyaxis."set_label_position"(pos)     # the guides
                pyaxis."set_ticks_position"("both")  # the hash marks
                getproperty(pyaxis, Symbol(:tick_, pos))()        # the tick labels
            end
            if axis[:guide_position] != :auto && letter != :z
                pyaxis."set_label_position"(axis[:guide_position])
            end
            py_set_scale(ax, sp, axis)
            axis[:ticks] == :native ? nothing : py_set_lims(ax, sp, axis)
            if ispolar(sp) && letter == :y
                ax."set_rlabel_position"(90)
            end
            ticks = sp[:framestyle] == :none ? nothing : get_ticks(sp, axis)
            # don't show the 0 tick label for the origin framestyle
            if sp[:framestyle] == :origin && length(ticks) > 1
                ticks[2][ticks[1] .== 0] .= ""
            end

            # Set ticks
            fontProperties = PyPlot.PyCall.PyDict(Dict("family" => axis[:tickfontfamily],
                                                       "size" => py_thickness_scale(plt, axis[:tickfontsize]),
                                                       # "useTex" => false,
                                                       "rotation" => axis[:tickfontrotation]))

            positions = getproperty(ax, Symbol("get_",letter,"ticks"))()
            pyaxis.set_major_locator(pyticker.FixedLocator(positions))
            getproperty(ax, Symbol("set_",letter,"ticklabels"))(positions, fontdict=fontProperties)

            # workaround to set mathtext.fontspec per Text element
            env = "\\mathregular"  # matches the outer fonts https://matplotlib.org/tutorials/text/mathtext.html

            axis[:ticks] == :native ? nothing : py_set_ticks(ax, ticks, letter, env)
            # Tick marks
            intensity = 0.5  # This value corresponds to scaling of other grid elements
            pyaxis."set_tick_params"(
                direction = axis[:tick_direction] == :out ? "out" : "in",
                width=py_thickness_scale(plt, intensity),
                length= 5 * py_thickness_scale(plt, intensity)
                                     )

            getproperty(ax, Symbol("set_", letter, "label"))(axis[:guide])
            if get(axis.plotattributes, :flip, false)
                getproperty(ax, Symbol("invert_", letter, "axis"))()
            end
            pyaxis."label"."set_fontsize"(py_thickness_scale(plt, axis[:guidefontsize]))
            pyaxis."label"."set_family"(axis[:guidefontfamily])

            if (RecipesPipeline.is3d(sp))
                pyaxis."set_rotate_label"(false)
            end

            if (letter == :y && !RecipesPipeline.is3d(sp))
                pyaxis."label"."set_rotation"(axis[:guidefontrotation] + 90)
            else
                pyaxis."label"."set_rotation"(axis[:guidefontrotation])
            end

            if axis[:grid] && !(ticks in (:none, nothing, false))
                fgcolor = py_color(axis[:foreground_color_grid])
                pyaxis."grid"(true,
                    color = fgcolor,
                    linestyle = py_linestyle(:line, axis[:gridstyle]),
                    linewidth = py_thickness_scale(plt, axis[:gridlinewidth]),
                    alpha = axis[:gridalpha])
                ax."set_axisbelow"(true)
            else
                pyaxis."grid"(false)
            end
            #

            if axis[:minorticks] > 1
                pyaxis."set_minor_locator"(PyPlot.matplotlib.ticker.AutoMinorLocator(axis[:minorticks]))
                pyaxis."set_tick_params"(
                    which = "minor",
                    direction = axis[:tick_direction] == :out ? "out" : "in",
                    width=py_thickness_scale(plt, intensity))
            end

            if axis[:minorgrid]
                if !(axis[:minorticks] > 1)  # Check if ticks were already configured
                    ax."minorticks_on"()
                end
                pyaxis."set_tick_params"(
                    which = "minor",
                    direction = axis[:tick_direction] == :out ? "out" : "in",
                    width=py_thickness_scale(plt, intensity))

                pyaxis."grid"(true,
                    which = "minor",
                    color = fgcolor,
                    linestyle = py_linestyle(:line, axis[:minorgridstyle]),
                    linewidth = py_thickness_scale(plt, axis[:minorgridlinewidth]),
                    alpha = axis[:minorgridalpha])
            end

            py_set_axis_colors(sp, ax, axis)
        end

        # showaxis
        if !sp[:xaxis][:showaxis]
            kw = KW()
            if ispolar(sp)
                ax.spines["polar"].set_visible(false)
            end
            for dir in (:top, :bottom)
                if !ispolar(sp)
                    ax.spines[string(dir)].set_visible(false)
                end
                kw[dir] = kw[Symbol(:label,dir)] = false
            end
            ax."xaxis"."set_tick_params"(; which="both", kw...)
        end
        if !sp[:yaxis][:showaxis]
            kw = KW()
            for dir in (:left, :right)
                if !ispolar(sp)
                    ax.spines[string(dir)].set_visible(false)
                end
                kw[dir] = kw[Symbol(:label,dir)] = false
            end
            ax."yaxis"."set_tick_params"(; which="both", kw...)
        end

        # aspect ratio
        aratio = get_aspect_ratio(sp)
        if aratio != :none
            ax."set_aspect"(isa(aratio, Symbol) ? string(aratio) : aratio, anchor = "C")
        end

        #camera/view angle
        if RecipesPipeline.is3d(sp)
            #convert azimuthal to match GR behaviour
            #view_init(elevation, azimuthal) so reverse :camera args
            ax."view_init"((sp[:camera].-(90,0))[end:-1:1]...)
        end

        # legend
        py_add_legend(plt, sp, ax)

        # this sets the bg color inside the grid
        getproperty(ax, set_facecolor_sym)(py_color(sp[:background_color_inside]))

        # link axes
        x_ax_link, y_ax_link = sp[:xaxis].sps[1].o, sp[:yaxis].sps[1].o
        ax != x_ax_link && ax."get_shared_x_axes"()."join"(ax, sp[:xaxis].sps[1].o)
        ax != y_ax_link && ax."get_shared_y_axes"()."join"(ax, sp[:yaxis].sps[1].o)
    end
    py_drawfig(fig)
end

# Set the (left, top, right, bottom) minimum padding around the plot area
# to fit ticks, tick labels, guides, colorbars, etc.
function _update_min_padding!(sp::Subplot{PyPlotBackend})
    ax = sp.o
    ax === nothing && return sp.minpad
    plotbb = py_bbox(ax)

    # TODO: this should initialize to the margin from sp.attr
    # figure out how much the axis components and title "stick out" from the plot area
    # leftpad = toppad = rightpad = bottompad = 1mm
    leftpad   = 0mm
    toppad    = 0mm
    rightpad  = 0mm
    bottompad = 0mm
    for bb in (py_bbox_axis(ax, "x"), py_bbox_axis(ax, "y"), py_bbox_title(ax), py_bbox_legend(ax))
        if ispositive(width(bb)) && ispositive(height(bb))
            leftpad   = max(leftpad,   left(plotbb) - left(bb))
            toppad    = max(toppad,    top(plotbb)  - top(bb))
            rightpad  = max(rightpad,  right(bb)    - right(plotbb))
            bottompad = max(bottompad, bottom(bb)   - bottom(plotbb))
        end
    end

    # optionally add the width of colorbar labels and colorbar to rightpad
    if haskey(sp.attr, :cbar_ax)
        bb = py_bbox(sp.attr[:cbar_handle]."ax"."get_yticklabels"())
        sp.attr[:cbar_width] = width(bb) + (sp[:colorbar_title] == "" ? 0px : 30px)
        rightpad = rightpad + sp.attr[:cbar_width]
    end

    # add in the user-specified margin
    leftpad   += sp[:left_margin]
    toppad    += sp[:top_margin]
    rightpad  += sp[:right_margin]
    bottompad += sp[:bottom_margin]

    dpi_factor = Plots.DPI / sp.plt[:dpi]

    sp.minpad = Tuple(dpi_factor .* [leftpad, toppad, rightpad, bottompad])
end

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

function py_add_annotations(sp::Subplot{PyPlotBackend}, x, y, val)
    ax = sp.o
    ax."annotate"(val, xy = (x,y), zorder = 999, annotation_clip = false)
end

function py_add_annotations(sp::Subplot{PyPlotBackend}, x, y, val::PlotText)
    ax = sp.o
    ax."annotate"(val.str,
        xy = (x,y),
        family = val.font.family,
        color = py_color(val.font.color),
        horizontalalignment = val.font.halign == :hcenter ? "center" : string(val.font.halign),
        verticalalignment = val.font.valign == :vcenter ? "center" : string(val.font.valign),
        rotation = val.font.rotation,
        size = py_thickness_scale(sp.plt, val.font.pointsize),
        zorder = 999,
        annotation_clip = false
    )
end

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

py_legend_pos(pos::Symbol) = get(
    (
        right = "right",
        left = "center left",
        top = "upper center",
        bottom = "lower center",
        bottomleft = "lower left",
        bottomright = "lower right",
        topright = "upper right",
        topleft = "upper left",
        outerright = "center left",
        outerleft = "right",
        outertop = "lower center",
        outerbottom = "upper center",
        outerbottomleft = "lower right",
        outerbottomright = "lower left",
        outertopright = "upper left",
        outertopleft = "upper right",
    ),
    pos,
    "best",
)
py_legend_pos(pos) = "lower left"

py_legend_bbox(pos::Symbol) = get(
    (
        outerright = (1.0, 0.5, 0.0, 0.0),
        outerleft = (-0.15, 0.5, 0.0, 0.0),
        outertop = (0.5, 1.0, 0.0, 0.0),
        outerbottom = (0.5, -0.15, 0.0, 0.0),
        outerbottomleft = (-0.15, 0.0, 0.0, 0.0),
        outerbottomright = (1.0, 0.0, 0.0, 0.0),
        outertopright = (1.0, 1.0, 0.0, 0.0),
        outertopleft = (-0.15, 1.0, 0.0, 0.0),
    ),
    pos,
    (0.0, 0.0, 1.0, 1.0),
)
py_legend_bbox(pos) = pos

function py_add_legend(plt::Plot, sp::Subplot, ax)
    leg = sp[:legend]
    if leg != :none
        # gotta do this to ensure both axes are included
        labels = []
        handles = []
        for series in series_list(sp)
            if should_add_to_legend(series)
                clims = get_clims(sp, series)
                # add a line/marker and a label
                push!(handles,
                    if series[:seriestype] == :shape || series[:fillrange] !== nothing
                        pypatches."Patch"(
                            edgecolor = py_color(single_color(get_linecolor(series, clims)), get_linealpha(series)),
                            facecolor = py_color(single_color(get_fillcolor(series, clims)), get_fillalpha(series)),
                            linewidth = py_thickness_scale(plt, clamp(get_linewidth(series), 0, 5)),
                            linestyle = py_linestyle(series[:seriestype], get_linestyle(series)),
                            capstyle = "butt"
                        )
                    elseif series[:seriestype] in (:path, :straightline, :scatter, :steppre, :steppost)
                        hasline = get_linewidth(series) > 0
                        PyPlot.plt."Line2D"((0, 1),(0,0),
                            color = py_color(single_color(get_linecolor(series, clims)), get_linealpha(series)),
                            linewidth = py_thickness_scale(plt, hasline * sp[:legendfontsize] / 8),
                            linestyle = py_linestyle(:path, get_linestyle(series)),
                            solid_capstyle = "butt", solid_joinstyle = "miter",
                            dash_capstyle = "butt", dash_joinstyle = "miter",
                            marker = py_marker(_cycle(series[:markershape], 1)),
                            markersize = py_thickness_scale(plt, 0.8 * sp[:legendfontsize]),
                            markeredgecolor = py_color(single_color(get_markerstrokecolor(series)), get_markerstrokealpha(series)),
                            markerfacecolor = py_color(single_color(get_markercolor(series, clims)), get_markeralpha(series)),
                            markeredgewidth = py_thickness_scale(plt, 0.8 * get_markerstrokewidth(series) * sp[:legendfontsize] / first(series[:markersize]))   # retain the markersize/markerstroke ratio from the markers on the plot
                        )
                    else
                        series[:serieshandle][1]
                    end
                )
                push!(labels, series[:label])
            end
        end

        # if anything was added, call ax.legend and set the colors
        if !isempty(handles)
            leg = ax."legend"(handles,
                labels,
                loc = py_legend_pos(leg),
                bbox_to_anchor = py_legend_bbox(leg),
                scatterpoints = 1,
                fontsize = py_thickness_scale(plt, sp[:legendfontsize]),
                facecolor = py_color(sp[:background_color_legend]),
                edgecolor = py_color(sp[:foreground_color_legend]),
                framealpha = alpha(plot_color(sp[:background_color_legend])),
                fancybox = false,  # makes the legend box square
                borderpad=0.8      # to match GR legendbox
            )
            frame = leg."get_frame"()
            frame."set_linewidth"(py_thickness_scale(plt, 1))
            leg."set_zorder"(1000)
            if sp[:legendtitle] !== nothing
                leg."set_title"(sp[:legendtitle])
                PyPlot.plt."setp"(leg."get_title"(), color = py_color(sp[:legendtitlefontcolor]), family = sp[:legendtitlefontfamily], fontsize = py_thickness_scale(plt, sp[:legendtitlefontsize]))
            end

            for txt in leg."get_texts"()
                PyPlot.plt."setp"(txt, color = py_color(sp[:legendfontcolor]), family = sp[:legendfontfamily], fontsize = py_thickness_scale(plt, sp[:legendfontsize]))
            end
        end
    end
end

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

# Use the bounding boxes (and methods left/top/right/bottom/width/height) `sp.bbox` and `sp.plotarea` to
# position the subplot in the backend.
function _update_plot_object(plt::Plot{PyPlotBackend})
    for sp in plt.subplots
        ax = sp.o
        ax === nothing && return
        figw, figh = sp.plt[:size]
        figw, figh = figw*px, figh*px
        pcts = bbox_to_pcts(sp.plotarea, figw, figh)
        ax."set_position"(pcts)

        if haskey(sp.attr, :cbar_ax) && RecipesPipeline.is3d(sp)   # 2D plots are completely handled by axis dividers
            cbw = sp.attr[:cbar_width]
            # this is the bounding box of just the colors of the colorbar (not labels)
            cb_bbox = BoundingBox(right(sp.bbox)-cbw - 2mm, top(sp.bbox) + 2mm, _cbar_width-1mm, height(sp.bbox) -  4mm)
            pcts = bbox_to_pcts(cb_bbox, figw, figh)
            sp.attr[:cbar_ax]."set_position"(pcts)
        end
    end
    PyPlot.draw()
end

# -----------------------------------------------------------------
# display/output

_display(plt::Plot{PyPlotBackend}) = plt.o."show"()

for (mime, fmt) in (
    "application/eps"         => "eps",
    "image/eps"               => "eps",
    "application/pdf"         => "pdf",
    "image/png"               => "png",
    "application/postscript"  => "ps",
    "image/svg+xml"           => "svg",
    "application/x-tex"       => "pgf"
)
    @eval function _show(io::IO, ::MIME{Symbol($mime)}, plt::Plot{PyPlotBackend})
        fig = plt.o
        fig."canvas"."print_figure"(
            io,
            format=$fmt,
            # bbox_inches = "tight",
            # figsize = map(px2inch, plt[:size]),
            facecolor = fig."get_facecolor"(),
            edgecolor = "none",
            dpi = plt[:dpi]
        )
    end
end

closeall(::PyPlotBackend) = PyPlot.plt."close"("all")


text:
## Figure 2 - robust ellipsoids

using Plots, Distributions, LaTeXStrings
pyplot()

# creates an ellipse
pts = Plots.partialcircle(0, 2π, 200)
x1, y1 = Plots.unzip(pts)
x2, y2 = Plots.unzip(pts)
x1 = 2x1 .+ 10
y1 += 3y1 .+ 5
x2 = 3x2 .+ 10
y2 += 5.5y2 .+ 5

a1 = rand(Normal(10, 2), 100)
a2 = rand(Normal(5, 3), 100)

scatter(
    a1,
    a2,
    label = "Observations",
    title = L"Distribution of ($a_{1}$, $a_2$)",
    xaxis = (L"a_1"),
    yaxis = (L"a_2"),
    color = 1
)
scatter!((10, 5), shape = :square, label = "Nominal value", color = 2)

savefig("Figure2.pdf")

plot!(Shape(x2, y2), label = L"\Gamma_2=1.5", opacity = .2)
plot!(Shape(x1, y1), label = L"\Gamma_1=1", opacity = .2)

savefig("Figure2_with_ellip.pdf")


text:
module UnfoldMakie

using Makie
using AlgebraOfGraphics
using Unfold

include("plot_design.jl")
include("plot_results.jl")

export plot_results
export plot
# Write your package code here.

end


text:
using Test
using ForestBiometrics

@testset "volume equations" begin
    @test doyle_volume(10.0,16) == 36.0
    @test scribner_volume(10.0,16) == 55.0
    @test international_volume(10.0,16) == 71.44
end


text:
###############################################################################################
#
#  Extract raw Markdown from source files and turn them into docs
#
#  - Use `names(module)` to get a list of exported symbols from the module
#  - Use __META__ to get all documented symbols along with path to the sourcefile
#  - Use `grep` to get line number in file
#    - because MD.meta only has line for Function and not Module or DataType
#    - because even when there is a line number, it is off by one
#  - Use readlines to read in the source file and pull out docs & end line number
#    - because julia compiles docs to incorrect Markdown objects, and discards the raw MD
#  - Re-execute all docstrings because they may include external files
#  - Write .md files and a TOC, add links back to repo with start/end lines for each object
#  - Create mkdocs.yml
#  - Run `mkdocs` to build html docs
#  - Cleanup (optional)
#
#
#              -- Because julia's own documentation libraries cannot handle all of markdown
#
###############################################################################################

###############################################################################################
###
### The following constants may be modified for based on your own config
###

# Module to be documented
const mod = "mPulseAPI"

# Temporary location for generated files, change this if you already have a directory called src/
const doc_src = "src"

# YAML config file for mkdocs. This will be generated.
const mkdocsy = "mkdocs.yml"

# Prefix configuration for mkdocs.  Page names will be appended to this
const mkdocs_config = Dict(
    :site_name => "$(mod).jl Documentation",
    :site_url  => "https://akamai.github.com/$(mod).jl/",
    :repo_url  => "https://github.com/akamai/$(mod).jl/",
    :site_favicon     => "favicon.ico",
    :extra_css => ["css/mkdocs.css"],
    :site_description => "Communicate with the mPulse Query & Repository REST APIs to fetch information about tenants and apps.",
    :copyright => "Akamai, Inc.",
    :docs_dir  => "src",
    :use_directory_urls => false,
    :theme     => "readthedocs",
    :markdown_extensions => [:admonition, :def_list, :attr_list, "toc:\n        permalink: True"],
)

# Don't change this
immutable Page
    name::AbstractString
    title::AbstractString
    pregenerated::Bool

    Page(name::AbstractString, title::AbstractString, pregenerated::Bool=false) = new(name, title, pregenerated)
end

# Pages to build:
# - name
# - title
const pages = Page[
    Page("index", mod),
    Page("apiToken", "How to generate an API Token", true),
    Page("RepositoryAPI", "Repository API"),
    Page("QueryAPI", "Query API"),
    Page("exceptions", "Exceptions"),
    Page("cache_utilities", "Internal Cache Utilities"),
]

###
### End of user configurable section
###
###############################################################################################

using Formatting

eval(parse("using $mod"))

Mod = eval(parse(mod))

function getSymbols(mod::Module; order=[Module, DataType, Function])
    exported = Dict( map( n -> (string(n) => getfield(mod, n)), names(mod) ) )

    declarator = Dict(Function => "(function )?", DataType => "(abstract|immutable|type) ", Module => "module ")

    # Now get all the symbols and mark the exported ones
    function symbol2dict(k)
        k_doc = Docs.doc(k)

        name  = replace(string(k), Regex("^$mod\."), "")

        typ   = typeof(k)

        file  = haskey(k_doc.meta, :path) ?
                    k_doc.meta[:path] :
                length(k_doc.content) == 1 && haskey(k_doc.content[1].meta, :path) ?
                    k_doc.content[1].meta[:path] :
                ""

        api_doc = ""

        if file != ""
            lines = open(readlines, file)
            line = find(x -> ismatch(Regex("^$(declarator[typ])$(name)"), x), lines)
            if length(line) == 0
                println(Regex("^$(declarator[typ]) +$(name)"))
                println(lines)
            end
            line = line[1]

            if typ == Module
                endline = line
                api_doc_start = findnext(lines, "\"\"\"\n", line+1)
                api_doc_end   = findnext(lines, "\"\"\"\n", api_doc_start+1)
                if api_doc_end - api_doc_start > 0
                    api_doc = eval(parse(join(lines[api_doc_start:api_doc_end], "")))
                end
            else
                endline = findnext(lines, "end\n", line)
                api_doc = findprev(lines, "\"\"\"\n", line-2)
                if api_doc > 0
                    api_doc = eval(parse(join(lines[api_doc:line-1], "")))
                end
            end
        else
            line = 0
        end

        api_doc = replace(api_doc, r"^   \*"m, "    *")

        return Dict(
            :name     => name,
            :type     => typ,
            :exported => haskey(exported, replace(string(k), Regex("^$mod\."), "")),
            :file     => replace(file, r"^.*/", ""),
            :line     => line,
            :endline  => endline,
            :doc      => api_doc
        )
    end

    symbols = map( symbol2dict, filter( k -> !isa(k, ObjectIdDict), collect( keys(mod.__META__) ) ) )

    expo_order = Dict(true => "1", false => "2")
    type_order = Dict(zip(order, 1:length(order)))

    sort!(symbols, by = x -> format("{1}.{2}.{3}.{4:04d}.{5}", expo_order[x[:exported]], type_order[x[:type]], x[:file], x[:line], x[:name]))

    return symbols
end

labels = Pair[Module => "", DataType => "Type", Function => "Function"]
exps   = Pair[true => "Exported", false => "Namespaced"]

symbols = getSymbols(Mod)
refids  = Dict(map(s -> (s[:name] => replace(s[:file], r"\.jl$", ".md") * lowercase(string("#", s[:type], "-", s[:name]))), symbols))

function replace_refs(m)
    ref = match(Regex("^\\[(`?(?:$mod\\.)?(\\w+)`?)\\]"), m)

    if ref == nothing
        return m
    end

    txt = ref.captures[1]
    ref = ref.captures[2]

    return "[$txt]($(refids[ref])){: .x-ref}"
end

cd(dirname(@__FILE__)) do

    dir_existed = true

    if !isdir(doc_src)
        mkdir(doc_src)
        dir_existed = false
    end

    open(mkdocsy, "w") do yml

        for (k, v) in mkdocs_config
            print(yml, k, ": ")

            if isa(v, AbstractArray)
                println(yml, mapfoldl(x -> "\n    - $x", *, v))
            else
                println(yml, v)
            end
        end
        println(yml, "pages:")

        for page in pages

            println("INFO    -  Processing $(page.name)")

            # Only generate pages that are not pre-generated
            page.pregenerated || open(joinpath(doc_src, page.name * ".md"), "w") do f
                println(f, """
                    # $(page.title)

                    """)

                local file_symbs = filter(x -> x[:file] == (page.name == "index" ? mod : page.name) * ".jl", symbols)
                if page.name != "index"
                    for s in file_symbs
                        println(f, "* [$(s[:name])]($(refids[s[:name]]))")
                    end
                end

                for (exported, ex_label) in exps
                    for (typ, ty_label) in labels

                        local symbs = filter(x -> x[:exported] == exported && x[:type] == typ, file_symbs)

                        if length(symbs) > 0
                            !isempty(ty_label) && println(f, "## $ex_label $(ty_label)s")

                            exported || println(f, """
                                !!! note
                                    The following methods are not exported by default. You may use them by explicitly
                                    importing them or by prefixing them with the `$(mod).` namespace.

                                """)

                            for s in symbs
                                println(f, """
                                    ##$(typ == Module ? "" : "#") $(lowercase(string(s[:type]))) `$(s[:name])`
                                    [$(s[:file])#$(s[:line])$(s[:endline] != s[:line] ? "-$(s[:endline])" : "")]($(mkdocs_config[:repo_url])tree/master/src/$(s[:file])#L$(s[:line])-L$(s[:endline])){: .source-link}
                                    """)

                                # Replace references with links to actual functions
                                s_doc = replace(s[:doc], Regex("\\[`?(?:$mod\\.)?\\w+`?\\]\\(@ref\\)"), replace_refs)
                                s_doc = replace(s_doc, r"^`(\w+)`"m, m -> (r = replace(m, "`", ""); haskey(refids, r) ? "[$m]($(refids[r]))" : m))

                                # Remove `docs/src/` from any links since we might have that in raw md in our functions
                                s_doc = replace(s_doc, r"/?docs/src/", "")

                                println(f, s_doc, "\n---\n")
                            end
                        end
                    end
                end

                if page.name == "index"
                    println(f, """
                        ## API Reference
                        """)

                    p = Page("", "")
                    for s in symbols
                        (s[:file] == mod * ".jl") && continue

                        if p.name * ".jl" != s[:file]
                            p = filter(p -> p.name * ".jl" == s[:file], pages)

                            if length(p) == 0
                                continue
                            else
                                p = p[1]
                            end

                            println(f, """

                                * [$(p.title)]($(p.name).md)
                                """)
                        end
                        println(f, "    * [$(s[:name])]($(refids[s[:name]]))")

                    end

                end
            end

            println(yml, "    - \"$(page.title)\": \"$(page.name).md\"")

        end

    end

    run(`mkdocs build -c -f $mkdocsy`)

    if any(x -> x=="--delete", ARGS)
        for p in pages
            path = joinpath(doc_src, p.name * ".md")
            println("INFO    -  Removing $path")
            rm(path)
        end

        if !dir_existed
            println("INFO    -  Removing $doc_src/")
            rm(doc_src)
        end

        println("INFO    -  Removing $mkdocsy")
        rm(mkdocsy)
    end

    if any(x -> x=="--add", ARGS)
        info("Adding all documentation changes in $(doc_src) to this commit.")
        success(`git add $(doc_src)`)
    end
end


text:
@doc doc"""
    TranslationGroup{T<:Tuple,F} <: GroupManifold{Euclidean{T,F},AdditionOperation}

Translation group $\mathrm{T}(n)$ represented by translation arrays.

# Constructor
    TranslationGroup(n₁,...,nᵢ; field = 𝔽)

Generate the translation group on
$𝔽^{n₁,\dots,nᵢ}$ = `Euclidean(n₁,...,nᵢ; field = 𝔽)`, which is isomorphic to group itself.
"""
const TranslationGroup{T<:Tuple,F} = GroupManifold{Euclidean{T,F},AdditionOperation}

function TranslationGroup(n::Int...; field::AbstractNumbers = ℝ)
    return TranslationGroup{Tuple{n...},field}(
        Euclidean(n...; field = field),
        AdditionOperation(),
    )
end

function show(io::IO, ::TranslationGroup{N,F}) where {N,F}
    print(io, "TranslationGroup($(join(N.parameters, ", ")); field = $(F))")
end


text:
function luc_seeder{T<:Integer}(n::T)
    a = -ones(BigInt, n)
    a[end] = big(n)
    return a
end

function luc{T<:Integer}(n::T, k::T)
    NFib(n, k, luc_seeder)
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
struct CmdExe <: AbstractShell
    prompt::AbstractString
end

CmdExe() = new("playground>")

function Base.run(shell::CmdExe, env::Environment)
    prompt = getprompt(shell, env)
    run(`cmd /K prompt $prompt`)
end


text:
""
    ispangram(input)

Return `true` if `input` contains every alphabetic character (case insensitive).

"""
ispangram(input) = prod([ch ∈ lowercase(input) for ch in ['a':'z';]])

max_stars_repo_name: Rigante-pl/Exercism-Julia

max_issues_repo_name: Rigante-pl/Exercism-Julia

max_forks_repo_name: Rigante-pl/Exercism-Julia