text:
# module LOLTools

abstract type AbstractDTO end

function (::Type{T})(nt::NT) where {T <: AbstractDTO, NT <: NamedTuple}
    T(get.(Ref(nt), fieldnames(T), nothing)...)
end

using HTTP, JSON2

function lol_api_server(region_or_platform::String)::HTTP.URI
    merge(HTTP.URI("https://"), host=string(region_or_platform, ".api.riotgames.com"), path="/")
end

function http_action(endpoint::HTTP.URI, path::String, headers::Vector{Pair{String,String}}=Vector{Pair{String,String}}(), query::Union{Nothing,Dict{String,String}}=nothing)
    url = string(merge(endpoint, path=path))
    HTTP.get(url, headers; query=query)
end

function nothing_in_event(caller::Function, resp, result::Union{Vector{Any}, Vector{String}, <:AbstractDTO, Array{<:AbstractDTO}, Set{<:AbstractDTO}})
    nothing
end

function call_api(::Type{T}, api_key::String, action::Function, endpoint::HTTP.URI, path::String, event::Function, caller::Function; query::Union{Nothing,Dict{String,String}}=nothing)::T where T
    headers = ["X-Riot-Token" => api_key]
    resp = action(endpoint, path, headers, query)
    data = JSON2.read(IOBuffer(resp.body))
    if T <: Array
        result = first(T.parameters).(data)
    elseif T <: Set
        P = first(T.parameters)
        result = Set{P}(P.(data))
    else
        result = T(data)
    end
    event(caller, resp, result)
    result
end

# module LOLTools


text:
using SafeTestsets
using ThreeBodyDecay

# ThreeBodyDecay structure
@safetestset "Test of the tbs structure" begin
    include("test_threebodymasses.jl")
    include("test_threebodyspins.jl")
    include("test_threebodyparities.jl")
    include("test_invariants.jl")
end

include("wignerrotationdispatch.jl")

@safetestset "four vectors" begin include("test_generation.jl") end

# angular functions
@safetestset "Angular Functions block" begin
    include("circle_relations.jl")
    include("test_representaiton_property.jl")
    include("check_lineshape_functions.jl")
end

@safetestset "Permutations" begin include("test_wigner_permutations.jl") end

@safetestset "Sum rule" begin include("wigner_angle_sumrules.jl") end

# Lineshapes
@safetestset "Lineshapes" begin
    # include("breit-wigner.jl")
    include("scattering_length_X3872.jl")
end

# JPC functions
@safetestset "Couplings" begin
    include("test_couplings.jl")
    include("test_recoupling.jl")
    include("test_ls_amplitude.jl")
    include("test_sum_over_polarization.jl")
end

# Integrals
@safetestset "Integrals" begin
    include("test_three_body_phase_space_integral.jl")
end
# examples
@safetestset "Complete example" begin include("test_Xic_decay.jl") end

@safetestset "square dalitz" begin
    include("test_transformations_of_squaredalitz.jl")
end

include("test_inference.jl")

ext: jl

max_stars_repo_name: mmikhasenko/ThreeBodyDecay

max_issues_repo_name: mmikhasenko/ThreeBodyDecay

max_forks_repo_name: mmikhasenko/ThreeBodyDecay


text:
for fn in (*, /, +, -, ^)
    fns = Symbol(fn)
    @eval $fns(x::TReal, y::TReal) = record!(x.tape, Call, $fn, (x, y))
    # @eval $fns(x::Real, y::TReal) = record!(y.tape, Call, $fn, (x, y))
    # @eval $fns(x::TReal, y::Real) = record!(x.tape, Call, $fn, (x, y))
    @eval $fns(x::TReal, y::Real) = $fn(x, constant(x.tape, y))
    @eval $fns(x::Real, y::TReal) = $fn(constant(y.tape, x), y)
end

^(x::TReal, p::Integer) = record!(x.tape, Call, ^, (x, constant(x.tape, p)))
-(x::TReal) = record!(x.tape, Call, -, (x,))

for fn in (sin, cos, Base.log, Base.exp, abs, abs2, sign, tanh, Base.sqrt)
    fns = Symbol(fn)
    @eval $fns(x::TReal) = record!(x.tape, Call, $fn, (x,))
end

# non-tracked ops

for fn in (>, >=, <, <=, ==, !=)
    fns = Symbol(fn)
    @eval Base.$fns(x::TReal, y::TReal) = $fn(x.val, y.val)
end


text:
# ---
# layout: post
# title: "π day"
# date: 2019-03-13 00:00:00 +0000
# categories: blog
# mathjax: true
# ---

# >In the UK we have started to celebrate π day (the 3rd month's 14th day) every year, even though we don't use the USA's date formatting convention of `monthnumber` followed by `daynumber`. But we can't really celebrate the 31st of April (31/4) or the 3rd of Quatember (?) (3/14), so we'll happily celebrate π day on 14/3 along with everyone else!

# >I set myself a challenge at the beginning of March: make a π-related image using Julia and the Luxor.jl package every day until π day. Some days it worked out well, others didn't, but I've gathered them all here anyway. This post has a fair few images, but not very much code or mathematical content.

# The images here are in low-resolution: they should be available on my [Flickr page](https://www.flickr.com/photos/153311384@N03/) at their full resolution if you want to download or re-use them.

# ### Day 1: Circle packing

# Circle packing may be a well-trodden path, but it always looks neat, and it's a nice easy start. You maintain a list of circles (center point and radius). Then you create a random circle, check it against all the other ones, draw it if it doesn't overlap, or reduce the radius and try again. It's not very efficient but you can set it going and go and make some coffee.

# To make the π shape appear, the code creates a path:

#md fontsize(480)
#md textoutlines("π", O, :path, halign=:center, valign=:middle)
#md πoutline = first(pathtopoly())

# then checks whether each circle's centerpoint is inside or outside the outline of the π shape:

#md isinside(pt, πoutline)

# and colors it accordingly.

# ![image label](IMAGEFOLDER/t-800.png)

# ### Day 2: Dry and wet

# I repeated myself today, thinking I could develop the circles a bit more, and ended up with this glossier wet-look version. The apparently random-looking shapes in the background are Bézier splodges that are supposed to be splashes...

# ![image label](IMAGEFOLDER/pi-reds-balls-wet-800.png)

# ### Day 3: π packing

# This is π packing rather than circle packing, although the code is again quite similar in outline: choose a point at random, find the largest font size at which the π character fits without overlapping others in the list, and then place it and add it to the list. The colors are a bit murky though.

# ![image label](IMAGEFOLDER/pi-swarm-3-800.png)

# ### Day 4: Rainbow

# Combining concentric circles and rainbow colors, this image shows about 350 digits of π.

# ![image label](IMAGEFOLDER/digits-of-pi-avenir-800.png)

# To generate the digits of π, I use this function:

function pidigits(n)
    result = BigInt[]
    k, a, b, a1, b1 = big.([2, 4, 1, 12, 4])
    while n > 0
        p, q, k = k^2, 2k + 1, k + 1
        a, b, a1, b1 = a1, b1, p * a  +  q * a1, p * b  +  q * b1
        d, d1 = a ÷ b, a1 ÷ b1
        while d == d1
            push!(result, d)
            n -= 1
            a, a1 = 10(a % b), 10(a1 % b1)
            d, d1 = a ÷ b, a1 ÷ b1
        end
    end
    return result
end

# It looks like witchcraft to me, but I understand that it's a "spigot" algorithm. I was hoping for a while that it was named after a Professor Spigot, but in fact it's describing the way the digits trickle out one by one like drops of water. It's quick enough for a thousand digits or so, but slows down a lot when you ask for 100_000 or more, probably due to the hard work that the big integer library has to do: even when you're just calculating the first 15 digits of π, the values of `a1` and `b1` are way over the limits of Int64s.

#md julia-1.1> @time pidigits(1000);
#md  0.014522 seconds (44.90 k allocations: 9.425 MiB, 28.97% gc time)

# The image might work better on white:

# ![image label](IMAGEFOLDER/digits-of-pi-avenir-on-white-800.png)

# Sometimes I wanted to check where certain sequences of digits appeared. I couldn't find a built-in function that looked for a sequence of digits in an array, but this worked well enough for my purposes:

function findsubsequence(needle, haystack)
    result = Int64[]
    for k in 1:length(haystack) - length(needle)
        if needle == view(haystack, k:k + length(needle) - 1)
            push!(result, k)
        end
    end
    return result
end

findsubsequence(str::String, digits) =
    findsubsequence(map(x -> parse(Int, x), split(str, "")), digits)

findsubsequence("999999", pidigits(2000)) # => [763]

# ### Day 5: Low-fat

# A chunky typeface like the Avenir Heavy I used yesterday is good for masking and clipping. But I wondered how the narrowest typeface would look. I found Briem Akademi, designed by Gunnlaugur Briem at the Royal Academy of Fine Arts in Copenhagen. Adobe's description says:

# >The most compressed version works best where legibility is less important than dramatic visual effect.

# and I like the abstract look even though it's almost illegible... Would this make nice bathroom tiles?

# ![image label](IMAGEFOLDER/many-digits-of-pi-briem-800.png)

# ### Day 6 Breakfast and Tiffany

# I'm still thinking about using typefaces. I'm a fan of Ed Benguiat's ITC Tiffany font, his nostalgic look back from the 1970s to the age of Edwardian elegance.

# ![image label](IMAGEFOLDER/pi-digits-appearing-800.png)

# It's easy to do this with tables. Like the circle packing, the code checks whether the coordinates of each table cell fall within a given perimeter, and changes the font accordingly.

# ### Day 7 Distinguished

# The excellent Colors.jl package has a function called `distinguishable_colors()` (which fortunately tab-completes). The help text says:

# > This uses a greedy brute-force approach to choose `n` colors that are maximally distinguishable. Given `seed` color(s), and a set of possible hue, chroma, and lightness values (in LCHab space), it repeatedly chooses the next color as the one that maximizes the minimum pairwise distance to any of the colors already in the palette.

# Much to do with color depends on the viewer's perception, but I think it works well here. It starts at the top left, and works from left to right. (That pesky decimal point defaults to using the previous color...) You can spot the Feynman point (`999999`) halfway down on the left (look for the six consecutive sandy brown squares), or the four purple sevens on the bottom row.

# ![image label](IMAGEFOLDER/pi-distinguishable_colors-800.png)

# I remembered to try to choose the color for the small labels (probably unreadable in the low-resolution PNG you see here) so that they're either light on dark, or dark on light.

#md ... r, g, b = color of square
#md gamma = 2.2
#md luminance = 0.2126 * r^gamma + 0.7152 * g^gamma + 0.0722 * b^gamma
#md (luminance > 0.5^gamma) ? sethue("black") : sethue("white")

# ### Day 8 Candy crush edition

# ![image label](IMAGEFOLDER/candy-crush.png)

# I must have seen an advert for CandyCrush yesterday, or perhaps all that talk of gamma and LCHab spaces caused a reaction, but this sugar rush of an image was the result. The SVG version looks tasty but is too big for this web page.

# ### Day 9 Like a circle in a spiral, a wheel within a wheel

# Arranging the sweets in a spiral looks tidy.

# ![image label](IMAGEFOLDER/pi-digits-in-spiral-balls-800.png)

# ### Day 10 π into circumference

# Luxor's `polysample()` function takes a polygon and samples it at regular intervals. This allows the following idea, where each point on a shape (here, the outline of the π character) is slowly moved to a matching location on the circular shape around the outside.

# ![image label](IMAGEFOLDER/pi-to-circle-800.png)

# For a point on the π border `p1`, and a matching point on the circumference polygon `p2`, the intermediate point is given by `between(p1, p2, n)`, where `n` is between 0 and 1.

# I like the almost 3D effect you get from this.

# ### Day 11 Charcoal

# Time for a charcoal sketch:

# ![image label](IMAGEFOLDER/pi-charcoal-1-800.png)

# The crinkly edges of the paper are made by the `polysample()` function on a rectangle then applying simplex-`noise()`-y nudges to the vertices. The paper is textured with `rule()`d lines, and there's some very low values for `setopacity()` smudges. Shifting the Bézier curve handles slightly for each iteration gives a brushy/sketchy feel. (It's fortunate I can copy and paste some of this code from drawings I've made before: I've learnt the hard way that it's better keep things than throw them away...)

# ### Day 12

# I ran out of time on this one, and there are still some problems with the text spacing. The idea is to have the infinite digits of π spiral into some fiery star with some space-y stuff. Probably not the sort of image I should be attempting at all with simple vector-based 2D graphics tools, but it feels like a challenge. Those wispy trails are the same as yesterday's brush strokes, but using custom `setdash()` dashing patterns.

# ![image label](IMAGEFOLDER/pi cosmic spiral-800.png)

# ### Day 13

# The idea here is to show which digit of π is the current leader, in terms of how many times that digit has appeared already. (Yes, a stupid idea, I know!) Then I couldn't decide on how many digits to show, so it's going to be an animated GIF showing the first 1000 digits. At the 200 digit mark poor old "7" is struggling at the back of the field, but the glory days are ahead - after 1000 digits, it's overtaken 0, 4, and 6.

# ![image label](IMAGEFOLDER/200-digits-of-pi-800.png)

# The animation turned into a video rather than a GIF, because I don't like the low resolution of GIFs today.

# And now of course I have to add a suitable audio soundtrack. Luckily I've recently been playing with George Datseris' [MIDI interface for Julia](https://github.com/JuliaMusic), so it was easy enough to make a musical version of the first 1000 digits of π, where the digits from 0 to 9 choose the appropriate note from a reasonably harmonious scale.

using MIDI

function savetrack(track, notes)
    file = MIDIFile()
    addnotes!(track, notes)
    addtrackname!(track, "a track")
    push!(file.tracks, track)
    writeMIDIFile("/tmp/sound-of-pi.mid", file)
end

scales = [46, 48, 51, 53, 55, 57, 58, 60, 62, 65, 67]

function generatetune!(notes)
    pos = 1
    dur = 80
    k = 1
    manypidigits = pidigits(1000)
    for i in manypidigits
        dur = k * 960
        pos += k * 960
        n = scales[i + 1]
        note = Note(n, 76, pos, dur)
        push!(notes, note)
    end
end

notes = Notes()
track = MIDITrack()
generatetune!(notes)
savetrack(track, notes)

# ![image label](IMAGEFOLDER/music-credits.png)

# This "sonification" (or "audification") is just for fun. For a more convincing critique of these sonifications than I can provide, watch the always entertaining [Tantacrul](https://www.youtube.com/watch?v=Ocq3NeudsVk)'s presentation on YouTube.

# And while you're on YouTube, the π video is on [my YouTube channel](https://www.youtube.com/channel/UCfd52kTA5JpzOEItSqXLQxg), and it's my entry for YouTube's Most Boring Video of 2019 competition, but I suspect it won't do very well—competition in this category is fierce, even if sometimes the contestants are unwilling participants.

# Happy π day!

# [2019-03-13]

# ![cormullion signing off](http://steampiano.net/cormullionknot.gif?piday){: .center-image}

using Literate                                                                 #src
# preprocess for notebooks                                                     #src
function setimagefolder(content)                                               #src
    content = replace(content, "IMAGEFOLDER" => "$IMAGEFOLDER")                #src
    return content                                                             #src
end                                                                            #src

# for Jupyter notebook, put images in subfolder                                #src
#IMAGEFOLDER = "images/piday"                                                   #src

#Literate.notebook("source/piday.jl", "notebooks", preprocess = setimagefolder) #src

# for Markdown/Jekyll notebook, put images in "/images"                        #src

IMAGEFOLDER = "/images/piday"                                                  #src

Literate.markdown("source/piday.jl", ".", name="_posts/2019-03-13-piday",      #src
 preprocess = setimagefolder,                                                  #src
 codefence = "{% highlight julia %}" => "{% endhighlight julia %}",            #src
 documenter=false)                                                             #src

#src


text:
using ItemGraphs: ItemGraph, add_edge!, edgeitems

using MuladdMacro
using ..LeapSeconds
using ..Periods
export @transform

abstract type Transformation end

const registry = ItemGraph{TimeScale, Transformation}()

macro transform(from::Symbol, to::Symbol, ep::Symbol, args...)
    trans = Symbol(from, "to", to)
    epoch = Expr(:escape, ep)
    func = quote
        @inline function (::$(esc(trans)))($epoch::Epoch{$from})::Epoch{$to}
            $(args[end])
        end
    end

    arglist = func.args[end].args[end].args[1].args[1].args
    for a in args[1:end-1]
        if a isa Expr && a.head == :(=)
            arg = Expr(:kw, Expr(:escape, a.args[1]), a.args[2:end]...)
        else
            arg = Expr(:escape, a)
        end
        push!(arglist, arg)
    end

    quote
        struct $trans <: Transformation end
        add_edge!(registry, $from, $to, $(esc(trans))())
        $func
        $(esc(Epoch)){$(esc(to))}(ep::Epoch{$from, T}, args...) where {T} = $trans()(ep, args...)
    end
end

"""
   taiutc(jd1, jd2)

Transform a two-part Julian date from `TAI` to `UTC`.

# Example

```jldoctest
julia> tai  = Epoch{TAI}(2.4578265e6, 0.30477440993249416)
2017-03-14T07:18:52.509 TAI
julia> AstroTime.Epochs.taiutc(tai.jd1, tai.jd2)
(2.4578265e6, 0.30434616919175345)
```
"""
@inline function taiutc(jd1, jd2)
    big1 = jd1 >= jd2
    if  big1
        a1 = jd1
        a2 = jd2
    else
        a1 = jd2
        a2 = jd1
    end

    u1 = a1
    u2 = a2
    for i in range(1,3)
        tai1, tai2 = utctai(u1, u2)
        u2 += a1 - tai1
        u2 += a2 - tai2
    end

    if  big1
        date = u1
        date1 = u2
    else
        date = u2
        date1 = u1
    end
    date, date1
end

"""
    tttai(jd1, jd2)

Transform a two-part Julian date from `TT` to `TAI`.

# Example

```jldoctest
julia> tt  = Epoch{TT}(2.4578265e6, 0.30477440993249416)
2017-03-14T07:18:20.325 TT
julia> AstroTime.Epochs.tttai(tt.jd1, tt.jd2)
(2.4578265e6, 0.30440190993249416)
```
"""
@inline function tttai(jd1, jd2)
    dtat = OFFSET_TT_TAI/SECONDS_PER_DAY;
    if jd1 > jd2
        jd1 = jd1
        jd2 -= dtat
    else
        jd1 -= dtat
        jd2 = jd2
    end
    jd1, jd2
end

"""
    taitt(jd1, jd2)

Transform a two-part Julian date from `TAI` to `TT`.

# Example

```jldoctest
julia> tai = Epoch{TAI}(2.4578265e6, 0.30440190993249416)
2017-03-14T07:18:20.325 TAI
julia> AstroTime.Epochs.taitt(tai.jd1, tai.jd2)
(2.4578265e6, 0.30477440993249416)
```
"""
@inline function taitt(jd1, jd2)
    dtat = OFFSET_TT_TAI / SECONDS_PER_DAY
    if jd1 > jd2
        jd1 = jd1
        jd2 = jd2 + dtat
    else
        jd1 = jd1 + dtat
        jd2 = jd2
    end
    jd1, jd2
end

"""
    ut1tai(jd1, jd2, dta)

Transform a two-part Julian date from `UT1` to `TAI`.

# Example

```jldoctest
julia> ut1 = Epoch{UT1}(2.4578265e6, 0.30440190993249416)
2017-03-14T07:18:20.325 UT1
julia> AstroTime.Epochs.ut1tai(ut1.jd1, ut1.jd2, AstroTime.Epochs.dut1(ut1)-AstroTime.Epochs.leapseconds(julian(ut1)))
(2.4578265e6, 0.3048243932182868)
```
"""
@inline function ut1tai(jd1, jd2, dta)
    dtad = dta / SECONDS_PER_DAY
    if jd1 > jd2
        date = jd1
        date1 = jd2 - dtad
    else
        date = jd1 - dtad
        date1 = jd2
    end
    date, date1
end

"""
    taiut1(jd1, jd2, dta)

Transform a two-part Julian date from `TAI` to `UT1`.

# Example

```jldoctest
julia> tai = Epoch{TAI}(2.4578265e6, 0.30440190993249416)
2017-03-14T07:18:20.325 TAI
julia> AstroTime.Epochs.ut1tai(tai.jd1, tai.jd2, AstroTime.Epochs.dut1(tai)-AstroTime.Epochs.leapseconds(julian(tai)))
(2.4578265e6, 0.30477440993249416)
```
"""
@inline function taiut1(jd1, jd2, dta)
    dtad = dta / SECONDS_PER_DAY
    if jd1 > jd2
        date = jd1
        date1 = jd2 + dtad
    else
        date = jd1 + dtad
        date1 = jd2
    end
    date, date1
end

"""
    tcgtt(jd1, jd2)

Transform a two-part Julian date from `TCG` to `TT`.

# Example

```jldoctest
julia> tcg = Epoch{TCG}(2.4578265e6, 0.30440190993249416)
2017-03-14T07:18:20.325 TCG
julia> AstroTime.Epochs.tcgtt(tcg.jd1, tcg.jd2)
(2.4578265e6, 0.30440190993249416)
```
"""
@inline function tcgtt(jd1, jd2)
    t77t = MOD_JD_77 + OFFSET_TT_TAI / SECONDS_PER_DAY
    if jd1 > jd2
        date = jd1
        date1 = jd2 - ((jd1 - MJD) + (jd2 - t77t)) * ELG
    else
        date = jd1 - ((jd2 - MJD) + (jd1 - t77t)) * ELG
        date1 = jd2
    end
    date, date1
end

"""
    tttcg(jd1, jd2)

Transform a two-part Julian date from `TT` to `TCG`.

# Example

```jldoctest
julia> tt = Epoch{TT}(2.4578265e6, 0.30440190993249416)
2017-03-14T07:18:20.325 TT
julia> AstroTime.Epochs.tttgc(tt.jd1, tt.jd2)
(2.4578265e6, 0.30440190993249416)
```
"""
@inline function tttcg(jd1, jd2)
    t77t = MOD_JD_77 + OFFSET_TT_TAI / SECONDS_PER_DAY
    elgg = ELG/(1.0-ELG)
    if jd1 > jd2
        date = jd1
        date1 = jd2 + ((jd1 - MJD) + (jd2 - t77t)) * elgg
    else
        date = jd1 + ((jd2 - MJD) + (jd1 - t77t)) * elgg
        date1 = jd2
    end
    date, date1
end

"""
    tdbtt(jd1, jd2, dtr)

Transform a two-part Julian date from `TDB` to `TT`.

# Example

```jldoctest
julia> tdb = Epoch{TDB}(2.4578265e6, 0.30440190993249416)
2017-03-14T07:18:20.325 TDB
julia> AstroTime.Epochs.tdbtt(tdb.jd1, tdb.jd2, AstroTime.Epochs.deltatr(tdb))
(2.4578265e6, 0.30440190993249416)
```
"""
@inline function tdbtt(jd1, jd2, dtr)
    dtrd = dtr / SECONDS_PER_DAY
    if jd1 > jd2
        date = jd1
        date1 = jd2 - dtrd
    else
        date = jd1 - dtrd
        date1 = jd2
    end
    date, date1
end

"""
    diff_tdb_tt(jd1, jd2)

Computes difference TDB-TT in seconds at time JD (julian days)
The timescale for the input JD can be either TDB or TT.

The accuracy of this routine is approx 40 microseconds in interval 1900-2100 AD.
Note that an accurate transformation betweem TDB and TT depends on the
trajectory of the observer. For two observers fixed on the earth surface
the quantity TDB-TT can differ by as much as about 4 microseconds.

### References ###

1. [https://www.cv.nrao.edu/~rfisher/Ephemerides/times.html#TDB](https://www.cv.nrao.edu/~rfisher/Ephemerides/times.html#TDB)
2. [Issue #26](https://github.com/JuliaAstro/AstroTime.jl/issues/26)

"""
function diff_tdb_tt(jd1, jd2)
    g = 357.53 + 0.9856003((jd1 - J2000) + jd2)
    0.001658sind(g) + 0.000014sind(2g)
end

diff_tdb_tt(ep::Epoch) = diff_tdb_tt(julian1(ep), julian2(ep))

function diff_tdb_tt(jd1, jd2, ut, elong, u, v)
    t = ((jd1 - J2000) + jd2) / DAYS_PER_MILLENNIUM
    # Convert UT to local solar time in radians.
    tsol = mod(ut, 1.0) * 2π  + elong

    # FUNDAMENTAL ARGUMENTS:  Simon et al. 1994.
    # Combine time argument (millennia) with deg/arcsec factor.
    w = t / 3600.0
    # Sun Mean Longitude.
    elsun = deg2rad(mod(280.46645683 + 1296027711.03429 * w, 360.0))
    # Sun Mean Anomaly.
    emsun = deg2rad(mod(357.52910918 + 1295965810.481 * w, 360.0))
    # Mean Elongation of Moon from Sun.
    d = deg2rad(mod(297.85019547 + 16029616012.090 * w, 360.0))
    # Mean Longitude of Jupiter.
    elj = deg2rad(mod(34.35151874 + 109306899.89453 * w, 360.0))
    # Mean Longitude of Saturn.
    els = deg2rad(mod(50.07744430 + 44046398.47038 * w, 360.0))
    # TOPOCENTRIC TERMS:  Moyer 1981 and Murray 1983.
    wt = 0.00029e-10 * u * sin(tsol + elsun - els) +
        0.00100e-10 * u * sin(tsol - 2.0 * emsun) +
        0.00133e-10 * u * sin(tsol - d) +
        0.00133e-10 * u * sin(tsol + elsun - elj) -
        0.00229e-10 * u * sin(tsol + 2.0 * elsun + emsun) -
        0.02200e-10 * v * cos(elsun + emsun) +
        0.05312e-10 * u * sin(tsol - emsun) -
        0.13677e-10 * u * sin(tsol + 2.0 * elsun) -
        1.31840e-10 * v * cos(elsun) +
        3.17679e-10 * u * sin(tsol)

    # =====================
    # Fairhead et al. model
    # =====================

    # T**0
    w0 = 0.0
    for j in eachindex(fairhd0)
        @muladd w0 += fairhd0[j][1] * sin(fairhd0[j][2] * t + fairhd0[j][3])
    end
    # T**1
    w1 = 0.0
    for j in eachindex(fairhd1)
        @muladd w1 += fairhd1[j][1] * sin(fairhd1[j][2] * t + fairhd1[j][3])
    end
    # T**2
    w2 = 0.0
    for j in eachindex(fairhd2)
        @muladd w2 += fairhd2[j][1] * sin(fairhd2[j][2] * t + fairhd2[j][3])
    end
    # T**3
    w3 = 0.0
    for j in eachindex(fairhd3)
        @muladd w3 += fairhd3[j][1] * sin(fairhd3[j][2] * t + fairhd3[j][3])
    end
    # T**4
    w4 = 0.0
    for j in eachindex(fairhd4)
        @muladd w4 += fairhd4[j][1] * sin(fairhd4[j][2] * t + fairhd4[j][3])
    end
    # Multiply by powers of T and combine.
    wf = @evalpoly t w0 w1 w2 w3 w4
    # Adjustments to use JPL planetary masses instead of IAU.
    wj = 0.00065e-6 * sin(6069.776754 * t + 4.021194) +
        0.00033e-6 * sin( 213.299095 * t + 5.543132) +
        (-0.00196e-6 * sin(6208.294251 * t + 5.696701)) +
        (-0.00173e-6 * sin(  74.781599 * t + 2.435900)) +
        0.03638e-6 * t * t
    # ============
    # Final result
    # ============
    # TDB-TT in seconds.
    w = wt + wf + wj
end

@inline function diff_ut1_tt(ep::Epoch)
    leapsec = leapseconds(julian(ep))
    ΔUT1 = dut1(ep)
    OFFSET_TT_TAI + leapsec - ΔUT1
end

"""
    tttdb(jd1, jd2, dtr)

Transform a two-part Julian date from `TT` to `TDB`.

# Example

```jldoctest
julia> tt = Epoch{TT}(2.4578265e6, 0.30440190993249416)
2017-03-14T07:18:20.325 TT
julia> AstroTime.Epochs.tttdb(tt.jd1, tt.jd2, AstroTime.Epochs.deltatr(tdb))
(2.4578265e6, 0.30440190993249416)
```
"""
@inline function tttdb(jd1, jd2, dtr)
    dtrd = dtr / SECONDS_PER_DAY
    if jd1 > jd2
        date = jd1
        date1 = jd2 + dtrd
    else
        date = jd1 + dtrd
        date1 = jd2
    end
    date, date1
end

"""
    ttut1(jd1, jd2, dt)

Transform a two-part Julian date from `TT` to `UT1`.

# Example

```jldoctest
julia> tt = Epoch{tt}(2.4578265e6, 0.30440190993249416)
2017-03-14T07:18:20.325 TT
julia> AstroTime.Epochs.ttut1(tt.jd1, tt.jd2, AstroTime.Epochs.deltat(tt))
(2.4578265e6, 0.30440190993249416)
```
"""
@inline function ttut1(jd1, jd2, dt)
     dtd = dt / SECONDS_PER_DAY
     if jd1 > jd2
        date = jd1
        date1 = jd2 - dtd
     else
        date = jd1 - dtd
        date1 = jd2
    end
    date, date1
end

"""
    ut1tt(jd1, jd2, dt)

Transform a two-part Julian date from `UT1` to `TT`.

# Example

julia> ut1 = Epoch{UT1}(2.4578265e6, 0.30440190993249416)
2017-03-14T07:18:20.325 UT1
julia> AstroTime.Epochs.ut1tt(ut1.jd1, ut1.jd2, AstroTime.Epochs.deltat(ut1))
(2.4578265e6, 0.30440190993249416)
```
"""
@inline function ut1tt(jd1, jd2, dt)
     dtd = dt / SECONDS_PER_DAY
     if jd1 > jd2
        date = jd1
        date1 = jd2 + dtd
     else
        date = jd1 + dtd
        date1 = jd2
     end
     date, date1
end

"""
    tdbtcb(jd1, jd2)

Transform a two-part Julian date from `TDB` to `TCB`.

# Example

julia> tdb = Epoch{TDB}(2.4578265e6, 0.30440190993249416)
2017-03-14T07:18:20.325 TDB
julia> AstroTime.Epochs.tdbtcb(tdb.jd1, tdb.jd2)
(2.4578265e6, 0.30440190993249416)
```
"""
@inline function tdbtcb(jd1, jd2)
    t77td = MJD + MOD_JD_77
    t77tf = OFFSET_TT_TAI/SECONDS_PER_DAY
    jd0 = TDB0/SECONDS_PER_DAY
    elbb = ELB/(1.0-ELB)
    if jd1 > jd2
        d = t77td - jd1
        f  = jd2 - jd0
        date = jd1
        date1 = f - ( d - ( f - t77tf ) ) * elbb
    else
        d = t77td - jd2
        f  = jd1 - jd0
        date = f + ( d - ( f - t77tf ) ) * elbb
        date1 = jd2
    end
    date, date1
end

"""
    tcbtdb(jd1, jd2)

Transform a two-part Julian date from `TCB` to `TDB`.

# Example

julia> ut1 = Epoch{TCB}(2.4578265e6, 0.30440190993249416)
2017-03-14T07:18:20.325 TCB1
julia> AstroTime.Epochs.tcbtdb(tcb.jd1, tcb.jd2)
(2.4578265e6, 0.30440190993249416)
```
"""
@inline function tcbtdb(jd1, jd2)
    t77td = MJD + MOD_JD_77
    t77tf = OFFSET_TT_TAI/SECONDS_PER_DAY
    jd0 = TDB0/SECONDS_PER_DAY
    if jd1 > jd2
        d = jd1 - t77td
        date = jd1
        date1 = jd2 + jd0 - ( d + ( jd2 - t77tf ) ) * ELB
    else
        d = jd2 - t77td;
        date = jd1 + jd0 - ( d + ( jd1 - t77tf ) ) * ELB
        date1 = jd2
    end
    date, date1
end

function jd2cal(jd1, jd2)
    dj = jd1 + jd2
    if dj < JD_MIN || dj > JD_MAX
        throw(ArgumentError("Julian date is outside of the representable range ($JD_MIN, $JD_MAX)."))
    end

    if jd1 >= jd2
        date = jd1
        date1 = jd2
    else
        date = jd2
        date1 = jd1
    end

    date1 -= 0.5

    f1 = mod(date, 1.0)
    f2 = mod(date1, 1.0)
    f = mod(f1 + f2, 1.0)
    if f < 0.0
        f += 1.0
    end
    d = round(date-f1) + round(date1-f2) + round(f1+f2-f)
    jd = Int(round(d) + 1)

    l = jd + 68569
    n = (4 * l) ÷ 146097
    l -= (146097 * n + 3) ÷ 4
    i = (4000 * (l + 1)) ÷ 1461001
    l -= (1461 * i) ÷ 4 - 31
    k = (80 * l) ÷ 2447
    id = Int(floor((l - (2447 * k) ÷ 80)))
    l = k ÷ 11
    im = Int(floor((k + 2 - 12 * l)))
    iy = Int(floor((100 * (n - 49) + i + l)))

    iy, im, id, f

end

function cal2jd(iy, im, id)
    EYEAR_ALLOWED = -4799
    MON_LENGTH = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]

    if iy < EYEAR_ALLOWED
        throw(ArgumentError("Year is outside of the representable range (< $EYEAR_ALLOWED)"))
    end

    if im < 1 || im > 12
        throw(ArgumentError("Month is outside of the range (1,12)"))
    end

    ly = ((im == 2 ) && !Bool(iy % 4!=0) && (Bool(iy % 100!=0) || !Bool(iy % 400!=0)))? 1:0 #check if leap year

    if ((id < 1) || (id > (MON_LENGTH[im] + ly)))
        throw(ArgumentError("Day is outside of permissible range (1, $(MON_LENGTH[im]))"))
    end

    my = (im - 14) ÷ 12
    iypmy = trunc(Int,(iy + my))
    jd = MJD
    jd1 = float((((1461 * (iypmy + 4800)) ÷ 4)
          + (367 * trunc(Int,(im - 2 - 12 * my))) ÷ 12
          - (3 * ((iypmy + 4900) ÷ 100)) ÷ 4
          + trunc(Int,id) - 2432076))

    jd, jd1
end

"""
   utctai(jd1, jd2)

Transform a two-part Julian date from `UTC` to `TAI`.

# Example

```jldoctest
julia> utc  = Epoch{UTC}(2.4578265e6, 0.30477440993249416)
2017-03-14T07:18:52.509 UTC
julia> AstroTime.Epochs.utctai(utc.jd1, utc.jd2)
(2.4578265e6, 0.3052026506732349)
```
"""
function utctai(jd1, jd2)
    big1 = jd1 >= jd2
    if big1
        u1 = jd1
        u2 = jd2
    else
        u1 = jd2
        u2 = jd1
    end

    iy, im, id, fd = jd2cal(u1, u2)
    u2 -= fd
    drift0 = leapseconds(u1 +  u2)
    drift12 = leapseconds(u1 + u2 + 0.5)
    drift24 = leapseconds(u1 + u2 + 1.5)

    dlod = 2.0 * (drift12 - drift0)
    dleap = drift24 - (drift0 + dlod)

    fd *= (SECONDS_PER_DAY + dleap)/SECONDS_PER_DAY
    fd *= (SECONDS_PER_DAY + dlod)/SECONDS_PER_DAY

    z1, z2 = cal2jd(iy, im, id)

    a2 = z1 - u1
    a2 += z2
    a2 += fd + drift0 / SECONDS_PER_DAY
    if big1
        date = u1
        date1 = a2
    else
        date = a2
        date1 = u1
    end
    date, date1
end

"""
    datetime2julian(scale::T, year, month, date, hour, min, sec) where {T <: TimeScale}

Transforms DateTime field to two-part Julian Date, special provision for leapseconds is provided.

#Example
```jldoctest
julia> AstroTime.Epochs.datetime2julian(UTC, 2016, 12, 31, 23, 59, 60)
(2.4577535e6, 0.9999884260598836)
julia> AstroTime.Epochs.datetime2julian(TT, 2017, 2, 1, 23, 59, 59)
(2.4577855e6, 0.999988425925926)
```
"""
function datetime2julian(scale::T, year, month, date, hour, min, sec) where {T <: TimeScale}

    jd = sum(cal2jd(year, month, date))
    seclim = 60.0
    adjusted_seconds_per_day = SECONDS_PER_DAY
    if hour < 0 || hour > 23
        throw(ArgumentError("The input hour value should be between 0 and 23"))
    end
    if min < 0 || min > 59
        throw(ArgumentError("The input minute value should be between 0 and 59"))
    end

    if scale == UTC
        dat0 = leapseconds(jd)
        dat12 = leapseconds(jd + 0.5)
        dat24 = leapseconds(jd + 1.5)

        dleap = dat24 - (2.0 * dat12 - dat0)
        adjusted_seconds_per_day = SECONDS_PER_DAY + dleap
        if hour == 23 && min == 59
            seclim += dleap
        end
    end

    if sec < 0
        throw(ArgumentError("The input second value should be greater than 0"))
    end

    if sec >= seclim
            throw(ArgumentError("Time exceeds the maximum seconds in $(year)-$(month)-$(date)"))
    end

    time  = ( 60.0 * ( 60.0 * hour + min )  + sec ) / adjusted_seconds_per_day
    jd, time
end

function d2tf(ndp, days)
    sign = (days >= 0.0 ) ? '+' : '-'
    a = SECONDS_PER_DAY * abs(days)

    if ndp < 0
        nrs = 1
        for n in range(1,-ndp)
            nrs *= (n == 2 || n == 4) ? 6 : 10
        end
        w = a / nrs
        a = nrs * round(w)
    end

    nrs = 1
    for n in range(1,ndp)
        nrs *= 10
    end
    rm = nrs * 60.0
    rh = rm * 60.0

    a = round(nrs * a)
    ah = floor(a / rh)
    a -= ah * rh
    am = floor(a / rm)
    a -= am * rm
    as = floor(a / nrs)
    af = floor(a - as * nrs)
    sign, Int(ah), Int(am), Int(as), Int(af)
end

"""
   utcut1(jd1, jd2, dut1, dat)

Transform a two-part Julian date from `UTC` to `UT1`.

# Example

```jldoctest
julia> utc  = Epoch{UTC}(2.4578265e6, 0.30477440993249416)
2017-03-14T07:18:52.509 UTC
julia> AstroTime.Epochs.utcut1(utc.jd1, utc.jd2, Epochs.dut1(utc), leapseconds(julian(utc)))
(2.4578265e6, 0.30477440993249416)
```
"""
@inline function utcut1(jd1, jd2, dut1, dat)
    jd = +(jd1, jd2)
    dta = dut1 - dat
    tai1, tai2 = utctai(jd1, jd2)
    date, date1 = taiut1(tai1, tai2, dta)
    date, date1
end

function julian2datetime(scale::T, ndp, jd1, jd2) where T <: TimeScale
    a1 = jd1
    b1 = jd2
    jd = jd1 + jd2

    year1, month1, day1, frac1 = jd2cal(a1, b1)

    leap = 0

    if scale == UTC
        dat0 = leapseconds(jd - frac1)
        dat12 = leapseconds(jd - frac1 +0.5)
        dat24 = leapseconds(jd + 1.5 - frac1)
        dleap = dat24 - (2.0 * dat12 - dat0)
        leap = dleap != 0.0
        if leap
            frac1 += frac1 * dleap / SECONDS_PER_DAY
        end
    end

    sign, hour, min, sec, fracd = d2tf(ndp, frac1)

    if hour > 23
        year2, month2, day2, frac2 = jd2cal(a1+1.5, b1-frac1)
        if !leap
            year1 = year2
            month1 = month2
            day1 = day2
            hour = 0
            min = 0
            sec = 0
        else
            if min > 0
                  year1 = year2
                  month1 = month2
                  day1 = day2
                  hour = 0
                  min = 0
                  sec = 0
            else
                  hour = 23
                  min = 59
                  sec = 60
            end
            if ndp < 0 && min == 60
                  year1 = year2
                  month1 = month2
                  day1 = day2
                  hour = 0
                  min = 0
                  sec = 0
            end
        end
  end

  year1, month1, day1, hour, min, sec, fracd
end

"""
   ut1utc(jd1, jd2, dut1)

Transform a two-part Julian date from `UT1` to `UTC`.

# Example

```jldoctest
julia> ut1  = Epoch{UT1}(2.4578265e6, 0.30477440993249416)
2017-03-14T07:18:52.509 UT1
julia> AstroTime.Epochs.ut1utc(ut1.jd1, ut1.jd2, Epochs.dut1(ut1))
(2.4578265e6, 0.3047686523910154)
```
"""
function ut1utc(jd1, jd2, dut1)
    duts = dut1
    big1 = jd1 >= jd2
    if big1
        u1 = jd1
        u2 = jd2
    else
        u1 = jd2
        u2 = jd1
    end
    d1 = u1
    dats1 = 0
    for i in -1:3
        d2 = u2 + float(i)
        year, month, day, frac = jd2cal(d1, d2)
        dats2 = leapseconds(d1 + d2 - frac)
        if i == - 1
            dats1 = dats2
        end
        ddats = dats2 - dats1
        if abs(ddats) >= 0.5
            if ddats * duts >= 0
                duts -= ddats
            end
            d1, d2 = cal2jd(year, month, day)
            us1 = d1
            us2 = d2 - 1.0 + duts / SECONDS_PER_DAY
            du = u1 - us1
            du += u2 - us2
            if  du > 0
                fd = du * SECONDS_PER_DAY / ( SECONDS_PER_DAY + ddats )
                duts += ddats * ( frac <= 1.0 ? frac : 1.0 )
            end
            break
        end
        dats1 = dats2
    end
    u2 -= duts / SECONDS_PER_DAY
    if big1
        date = u1
        date1 = u2
    else
        date = u2
        date1 = u1
    end
    date, date1
end

# TAI <-> UTC
@transform UTC TAI ep begin
    jd1, jd2 = julian1(ep), julian2(ep)
    date, date1 = utctai(jd1, jd2)
    TAIEpoch(date, date1)
end

@transform TAI UTC ep begin
    jd1, jd2 = julian1(ep), julian2(ep)
    date, date1 = taiutc(jd1, jd2)
    UTCEpoch(date, date1)
end

# UTC <-> UT1
@transform UT1 UTC ep begin
    jd1, jd2 = julian1(ep), julian2(ep)
    date, date1 = ut1utc(jd1, jd2, dut1(ep))
    UTCEpoch(date, date1)
end

@transform UTC UT1 ep begin
    jd1, jd2 = julian1(ep), julian2(ep)
    date, date1 = utcut1(jd1, jd2, dut1(ep), leapseconds(julian(ep)))
    UT1Epoch(date, date1)
end

# TAI <-> UT1
@transform UT1 TAI ep Δt=dut1(ep)-leapseconds(ep) begin
    jd1, jd2 = julian1(ep), julian2(ep)
    date, date1 = ut1tai(jd1, jd2, Δt)
    TAIEpoch(date, date1)
end

@transform TAI UT1 ep Δt=dut1(ep)-leapseconds(ep) begin
    jd1, jd2 = julian1(ep), julian2(ep)
    date, date1 = taiut1(jd1, jd2, Δt)
    UT1Epoch(date, date1)
end

# TT <-> UT1
@transform UT1 TT ep Δt=diff_ut1_tt(ep) begin
    jd1, jd2 = julian1(ep), julian2(ep)
    date, date1 = ut1tt(jd1, jd2, Δt)
    TTEpoch(date, date1)
end

@transform TT UT1 ep Δt=diff_ut1_tt(ep) begin
    jd1, jd2 = julian1(ep), julian2(ep)
    date, date1 = ttut1(jd1, jd2, Δt)
    UT1Epoch(date, date1)
end

# TAI <-> TT
@transform TT TAI ep begin
    jd1, jd2 = julian1(ep), julian2(ep)
    date, date1 = tttai(jd1, jd2)
    TAIEpoch(date, date1)
end

@transform TAI TT ep begin
    jd1, jd2 = julian1(ep), julian2(ep)
    date, date1 = taitt(jd1, jd2)
    TTEpoch(date, date1)
end

# TT <-> TCG
@transform TCG TT ep begin
    jd1, jd2 = julian1(ep), julian2(ep)
    date, date1 = tcgtt(jd1, jd2)
    TTEpoch(date, date1)
end

@transform TT TCG ep begin
    jd1, jd2 = julian1(ep), julian2(ep)
    date, date1 = tttcg(jd1, jd2)
    TCGEpoch(date, date1)
end

# TT <-> TDB
@transform TDB TT ep Δtr=diff_tdb_tt(ep) begin
    jd1, jd2 = julian1(ep), julian2(ep)
    date, date1 = tdbtt(jd1, jd2, Δtr)
    TTEpoch(date, date1)
end

@transform TT TDB ep Δtr=diff_tdb_tt(ep) begin
    jd1, jd2 = julian1(ep), julian2(ep)
    date, date1 = tttdb(jd1, jd2, Δtr)
    TDBEpoch(date, date1)
end

# TDB <-> TCB
@transform TCB TDB ep begin
    jd1, jd2 = julian1(ep), julian2(ep)
    date, date1 = tcbtdb(jd1, jd2)
    TDBEpoch(date, date1)
end

@transform TDB TCB ep begin
    jd1, jd2 = julian1(ep), julian2(ep)
    date, date1 = tdbtcb(jd1, jd2)
    TCBEpoch(date, date1)
end

"""
    Epoch{T}(ep::Epoch{S}) where {T}, S}

Convert an `Epoch` with timescale `S` to an `Epoch` with timescale `T`.

# Example

```jldoctest
julia> Epoch{TT}(Epoch{TAI}(2000, 1, 1))
2000-01-01T00:00:32.184 TT
```
"""
function Epoch{T}(ep::Epoch{S}) where {T,S}
    transformations = edgeitems(registry, S, T)
    rescale(ep, transformations...)
end

Epoch{T}(ep::Epoch{T}) where {T} = ep

@generated function rescale(ep, transformations...)
    ex = :(ep)
    for trans in transformations
        ex = :($trans()($ex))
    end
    ex
end


text:
#   This file is part of Fatou.jl. It is licensed under the MIT license
#   Copyright (C) 2019 Michael Reed

using PyPlot
import PyPlot: plot

function plot(K::FilledSet;c::String="",bare::Bool=false)
    # plot figure using imshow based in input preferences
    figure()
    isempty(c) && (c = K.meta.cmap)
    isempty(c) ? PyPlot.imshow(K.meta.iter ? K.iter : K.mix, extent=K.meta.∂) :
        PyPlot.imshow(K.meta.iter ? K.iter : K.mix, cmap=c, extent=K.meta.∂)
    tight_layout()
    if !bare
        # determine if plot is Iteration, Roots, or Limit
        typeof(K.meta.iter ? K.iter : K.mix) == Matrix{UInt8} ? t = L"iter. " :
            K.meta.m==1 ? t = L"roots" : t = L"limit"
        # annotate title using LaTeX
        ttext = "f:z\\mapsto $(rdpm(Algebra.latex(K.meta.E))),\\,"
        if K.meta.newt
            PyPlot.title(latexstring("$ttext m = $(K.meta.m), ")*t)
            # annotate y-axis with Newton's method
            ylabel(L"Fatou\,set:\,"*L"z\,↦\,z-m\,×\,f(z)\,/\,f\,'(z)")
        else
            PyPlot.title(latexstring(ttext)*t)
        end
        tight_layout()
        colorbar()
    end
end

function orbit(E,f::Function,bi::Matrix{Float64},orb::Int=0,depth::Int=1,incr::Int=384; plt::Function=plot)
    x,N,N2,orbit,bis = real_orb(E,f,bi,orb,depth,incr)
    # prepare for next figure
    figure()
    # plot background lines
    plot(x[:],N[:,1],"k--",x[:],N[:,2])
    # plot orbit cobweb path
    plot(orbit[:,1],orbit[:,2],"r")
    # plot f^2,f^3,f^4,...
    for h ∈ 3:depth+1
        plot(x,N[:,h],lw=1)
    end
    if ~(orb == 0)
        plt(range(bi[1],stop=bi[2],length=length(N2)),N2[:],"gray",marker="x",linestyle=":",lw=1)
        funs = [latexstring("\\phi(x_{0:$orb})")]
        funt = ", IC: \$ x_0 = $(bis[3])\$, \$ n\\in0:$orb\$"
    else
        funs = []
        funt = ""
    end
    # trim graph
    d=1.07
    xlim(bi[1],bi[2])
    ylim(minimum([d*minimum(N[:,2]),0]),maximum([d*maximum(N[:,2]),0]))
    # set title
    fune = rdpm(Algebra.latex(E))
    PyPlot.title(latexstring("\$ x \\mapsto $fune\$$funt"))
    # set legend
    legend(vcat([L"$y=x$",L"$\phi(x)$",L"(x_n,\phi(x_n))"],[latexstring("\\phi^{$x}(x)") for x ∈ 2:depth],funs))
    tight_layout()
    return
end


text:
#Cubic
doubleRotatedDensity = rotateStructure(rotateStructure(densityCube, 1.5, 1.5, 1.5), -1.5, -1.5, -1.5)
@test similarity(densityCube, doubleRotatedDensity) > 0.85


text:
module TimeSeriesA2P2

#Use this to print
#=
using Weave
codeName = "TimeSeriesA2P2"
weave(Pkg.dir("$(pwd())\\$(codeName).jl"),
  informat="script",
  out_path = "$(pwd())\\output\\$(codeName)_Appendix.html",
  doctype = "md2html")
=#

if pwd() ∉ LOAD_PATH
  push!(LOAD_PATH,pwd())
end

using DataFrames, Distributions, StatsBase, GZip, JLD,
  Gadfly, CTIO, CTReg, JuMP, NLopt, HCubature, StaticArrays, CTNumerical,
  ParallelAccelerator, Measures, Formatting

const DATA_PATH = pwd() * "\\data" #path where we store the data
const OUTPUT_PATH = pwd() * "\\output"
const YAHOO_DATE_FORMAT = "m/d/yyyy"
const SP500_NAME = "GSPC"
const DBV_NAME = "DBV"

const FOOTER_NAME = "footer.tex" #these are just headers for the tex tables
const HEADER_NAME = "header.tex"
const DECIMALS = 3
const PlotContainer = Union{Plot,Gadfly.Compose.Context,Vector{Plot}} #holds graphs
const OVERRIDE_STAR_STRINGS = ["\\ostar","\\dstar","\\tstar"]

#asignment specific constants
const FACTOR_NAME = "FMFactors"
const ASSET_NAME = "PORTSBM55"
const CUTOFF_DATE = 1963_07
const FACTOR_SYMS = [:Mkt_RF, :SMB, :HML]
const RFR_SYM = :RF
const MKT_SYM = :Mkt

BLAS.set_num_threads(1)

  n2s(x::T where T<:Real) = num2Str(x, DECIMALS, Ints=true, scaleHurdle=99.)

#pre-process the asset data
function preProcessMonthlyFA()::Void

  #read in the data (note the non-standard NA values)
  assetDF::DataFrame = readtable("$DATA_PATH\\$ASSET_NAME.csv", nastrings=["-999", "-99.99"])
  factorDF::DataFrame = readtable("$DATA_PATH\\$FACTOR_NAME.csv", nastrings=["-999", "-99.99"])

  factorDF[:,:Mkt] = factorDF[:,:Mkt_RF] .+ factorDF[:,:RF]

  #merge the dataframes
  DF::DataFrame = join(factorDF, assetDF, on = :DATE, kind=:inner)

  #write out a JLS binary file
  stream::IOStream = open("$DATA_PATH\\$(FACTOR_NAME)_$(ASSET_NAME).jls", "w")
  serialize(stream, DF)
  close(stream)
  return nothing
end

#gets the dataframe as a serial file
function getDF(pathStr::String, refreshData::Bool, preProcessFunc::Function)::DataFrame

  #re-process the data if needed or desired
  if !isfile(pathStr) || refreshData
    preProcessFunc()
  end

  #read the binary file
  stream::IOStream = open(pathStr)
  DF::DataFrame = deserialize(stream)
  close(stream)

  return DF
end

mutable struct PCA
  λ::Vector{Float64}
  ν::Vector{Vector{Float64}}
  svd::Base.LinAlg.SVD
  labels::Vector{Symbol}
  labelsInd::Dict
  normConstant::Float64
end

function PCA(Σ::Matrix{Float64},
    labels::Vector{Symbol} = ((i::Int)->Symbol("v$i")).(1:size(Σ,1)))::PCA
  #setup and factorize
  M::Int, N::Int = size(Σ)
  labelsInd::Dict = Dict(labels[i]=>i for i::Int ∈ 1:M)

  Σsvd::Base.LinAlg.SVD = svdfact(Σ)

  #recover the eigenvalues
  λ::Vector{Float64} = Σsvd[:S]
  normConstant::Float64 = sum(λ)
  λ ./= normConstant

  #recover the eigenvectors
  ν::Vector{Vector{Float64}} =
    ((i::Int)->(Σsvd[:U][:,i]./sum((abs2).(Σsvd[:U][:,i])))).(1:N)

  return PCA(λ, ν, Σsvd,labels,labelsInd,normConstant)
end

function PCA(DF::DataFrame, labels::Vector{Symbol})
  D::Matrix{Float64} = Matrix{Float64}(size(DF[:,labels]))
  D .= Matrix{Float64}(DF[:,labels])
  D .-= mean(D,1)
  return PCA(BLAS.gemm('T','N',1.0, D,D), labels) #(This is simply D'D)
end

#adds principal componenets to an existing dataframe
function PC(pca::PCA, DF::DataFrame;
    labels::Vector{Symbol}=pca.labels,
    dataLabels::Vector{Symbol}=labels,
    PCLabels::Vector{Symbol} = ((i::Int)->(Symbol("PC$(i)"))).(1:length(labels)),
    )::Tuple{DataFrame, Vector{Symbol}}

    D::Matrix{Float64} = DF[:,dataLabels]
    for i::Int ∈ 1:length(PCLabels)
      DF[:,PCLabels[i]] = BLAS.gemv('N', 1.0, D, pca.ν[pca.labelsInd[labels[i]]])
    end

    return (DF, PCLabels)
end

#Calculates the summary statistics
function summaryA2P2(FADF::DataFrame, PCLabels::Vector{Symbol}, cols::Vector{Symbol})
  #label the table
  colNames::Vector{Vector{String}} =
    [((s::Symbol)->replace(String(s),"_","")).(cols)]
  numCols = length(colNames[end])

  rowNames::Vector{String} = ((s::Symbol)->replace(String(s),"_","")).(PCLabels)

  numRows::Int = length(rowNames)
  content::Vector{Vector{String}} = [Vector{String}(numCols) for i::Int ∈ 1:numRows]

  #run the tests and record the results
  for r::Int ∈ 1:numRows, c::Int ∈ 1:numCols
    content[r][c] = n2s(cor(FADF[:,PCLabels[r]], FADF[cols[c]]))
  end

  #write the table
  summaryTable::String = texTable( "Correlations of Principal Components",
    """See Tex File""", #caption
    colNames, #colNames
    Vector{String}(),#contentRowNames
    Vector{Matrix{String}}(), #content
    rowNames, #descRowNames
    content, #descContent
    Vector{String}(), #notes
  )

  return summaryTable
end

#problem specific visualizations
function visualizeA2P2(pca::PCA, FADF::DataFrame, PCLabels::Vector{Symbol})::Void

  #convenience containers
  K::Int = length(pca.λ)
  plots::Vector{PlotContainer} = Vector{PlotContainer}()
  plotNames::Vector{String} = Vector{String}()

  Gadfly.push_theme(:default)

  #some basic data processing for the eigenvalues
  λLabels::Vector{Symbol} = ((i::Int)->(Symbol("λ$i"))).(1:K)
  λDataLabels::Vector{String} = (n2s).(pca.λ)
  λCum::Vector{Float64} = ((i::Int)->sum(pca.λ[1:i])).(1:K)
  λCumDataLabels::Vector{String} = (n2s).(λCum)

  push!(plotNames, "eigenvalues")
  push!(plots, plot(y=pca.λ, x=λLabels, label=λDataLabels, Geom.bar, Geom.label(position=:above),
    Guide.xlabel("λ"),Guide.ylabel("λᵢ / sum(λ)"),Coord.Cartesian(ymin=0., ymax=1.05),
    Guide.title("Eigenvalues λ")))

  push!(plotNames, "cumulative eig")
  push!(plots, plot(y=λCum, x=λLabels, label=λCumDataLabels,
    Geom.point, Geom.line, Geom.label(position=:above),
    Guide.xlabel("λ"),Guide.ylabel("Cumulative λᵢ / sum(λ)"),Coord.Cartesian(ymin=0., ymax=1.05),
    Guide.title("Cumulative Eigenvalues λ")))

  PC3DF::DataFrame = melt(DataFrame(labels=pca.labels,
    ν1=pca.ν[1], ν2=pca.ν[2], ν3=pca.ν[3]),:labels)
  push!(plotNames, "eigenvectors")
  push!(plots, plot(PC3DF, y=:value, x=:labels, color=:variable,
    Geom.bar(position=:dodge),# Geom.line, Geom.point,
    Guide.xlabel("λ"), Guide.ylabel("νᵢ/sum(ν.²)"),
    Guide.title("Eigenvectors ν")))

  push!(plotNames, "eigenvectors ν₁ vs ν₂")
  push!(plots,
    plot(x=pca.ν[1], y=pca.ν[2], label=(string).(pca.labels),
      Geom.point,Geom.label,
      Guide.xlabel("ν₁"), Guide.ylabel("ν₂"),
      Coord.Cartesian(ymin=-0.5, ymax=0.5, xmin=-0.5, xmax=0.5),
      Guide.title("Eigenvectors ν₁ vs ν₂")))
  push!(plotNames, "eigenvectors ν₂ vs ν₃")
  push!(plots,
    plot( x=pca.ν[2], y=pca.ν[3], label=(string).(pca.labels),
      Geom.point,Geom.label,# Geom.line, Geom.point,
      Guide.xlabel("ν₂"), Guide.ylabel("ν₃"),
      Coord.Cartesian(ymin=-0.5, ymax=0.5, xmin=-0.5, xmax=0.5),
      Guide.title("Eigenvectors ν₂ vs ν₃")))

  push!(plotNames, "eigenvectors ν₃ vs ν₁")
  push!(plots,
    plot(x=pca.ν[3], y=pca.ν[1], label=(string).(pca.labels),
      Geom.point,Geom.label,# Geom.line, Geom.point,
      Guide.xlabel("ν₃"), Guide.ylabel("ν₁"),
      Coord.Cartesian(ymin=-0.5, ymax=0.5, xmin=-0.5, xmax=0.5),
      Guide.title("Eigenvectors ν₃ vs ν₁")))

  #write the graphs
  for i::Int ∈ 1:length(plots)
    draw(SVG("$OUTPUT_PATH\\$(plotNames[i]).svg", 9inch, 7inch), plots[i])
  end

  return nothing

end

#problem specific code
function A2P2Script(refreshData::Bool = false)::Void
  FADF::DataFrame = getDF("$DATA_PATH\\$(FACTOR_NAME)_$(ASSET_NAME).jls",
    refreshData, preProcessMonthlyFA)
  FADF = FADF[FADF[:,:DATE].≥196307,:]

  #rename the columns to a common pattern
  rename!(FADF, [:SMALL_LoBM, :SMALL_HiBM, :BIG_LoBM, :BIG_HiBM],
    [:ME1_BM1, :ME1_BM5, :ME5_BM1, :ME5_BM5])

  #get portfolio handles
  FSyms::Vector{Symbol} = FACTOR_SYMS
  RFRSym::Symbol = RFR_SYM
  RSyms::Vector{Symbol} =
    names(FADF[:,setdiff(names(FADF), [:DATE; RFR_SYM; MKT_SYM; FACTOR_SYMS])])

  #do the eigenvalue decomposition
  Σpca::PCA = PCA(FADF, RSyms)

  #get the principal components
  (FADF, PCLabels::Vector{Symbol}) = PC(Σpca, FADF)
  visualizeA2P2(Σpca, FADF, PCLabels)

  summaryTable::String = summaryA2P2(FADF,PCLabels[1:5],[RFR_SYM; MKT_SYM; FACTOR_SYMS])

  writeTables2File([summaryTable],
    HEADER_NAME, FOOTER_NAME, path=OUTPUT_PATH,
    outName = "TimeSeries A2P2 corTable.tex")

  return nothing
end

#uncomment to run
@time begin
  A2P2Script()
end

end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
# Vehicle ownership model, based on the NHTS

# comment out when using --track-allocations so not so many mem files
using Distributed; addprocs()
@everywhere begin
    using Pkg
    Pkg.activate(Base.source_dir())
end

using CSV, CodecZlib, Profile, DiscreteChoiceModels, Dagger, Optim, Tables, Infiltrator, DataFrames

# to save space in the repo, data are gzipped. helper function to read a
# gzipped csv "straight" into JuliaDB (via a temporary file)
function read_gzipped_csv(T, file; cs=30000)
    open(GzipDecompressorStream, file) do is
        mktemp() do path, os
            write(os, read(is))
            if T == DTable
                DTable(Tables.rowtable(CSV.File(path)), cs)
            else
                T(CSV.File(path))
            end
        end
    end
end

function read_data(T, basepath; cs=30000)
    data = read_gzipped_csv(T, joinpath(basepath, "hhpub.csv.gz"), cs=cs)

    # topcode hhvehcnt
    if T == DTable
        data = map(r -> merge((hhveh_topcode=min(r.HHVEHCNT, 4),), pairs(r)), data)
    else
        data.hhveh_topcode = min.(data.HHVEHCNT, 4)
    end

    data
end

function run_model(data)
    multinomial_logit(
        @utility(begin
            0 ~ α0
            1 ~ α1 + β1homeown * (HOMEOWN == 1) + @dummy_code(β1, HH_RACE, [2, 3, 4, 5, 6, 97]) + β1_hhsize * HHSIZE
            2 ~ α2 + β2homeown * (HOMEOWN == 1) + @dummy_code(β2, HH_RACE, [2, 3, 4, 5, 6, 97]) + β2_hhsize * HHSIZE
            3 ~ α3 + β3homeown * (HOMEOWN == 1) + @dummy_code(β3, HH_RACE, [2, 3, 4, 5, 6, 97]) + β3_hhsize * HHSIZE
            4 ~ α4plus + β4plushomeown * (HOMEOWN == 1) + @dummy_code(β4plus, HH_RACE, [2, 3, 4, 5, 6, 97]) + β4plus_hhsize * HHSIZE

            α0 = 0, fixed
        end),
        :hhveh_topcode,
        data;
        verbose=:medium,
        method=Newton()
    )
end

function main()
    data = read_data(DataFrame, "test/data/nhts", cs=30000)
    model = @time run_model(data)

    println(summary(model))

end

rechunk(table, chunksize) = DTable(collect(table), chunksize)

main()


text:
module PlasmoPlots

using Statistics
using Plasmo
using Plots
using NetworkLayout
using GeometryBasics:Point2f0,Point
using Colors
using LightGraphs

export plot, spy, layout_plot, matrix_plot

include("layout_plot.jl")

include("matrix_plot.jl")

end


text:
using VectorizedRNG: local_rng

@inline function Random.rand!(A::AbstractStrideArray, args::Vararg{Any,K}) where {K}
  rand!(local_rng(), A, args...)
end
@inline function Random.randn!(A::AbstractStrideArray, args::Vararg{Any,K}) where {K}
  randn!(local_rng(), A, args...)
end
## ignore type...
@inline Random.rand!(A::AbstractStrideArray, ::Type{T}) where {T} = rand!(local_rng(), A)

function rand_expr(expr, args...)
  array = :(StrideArray(undef))
  for i = 2:length(expr.args)
    aᵢ = expr.args[i]
    if aᵢ isa Integer
      push!(array.args, StaticInt(Int(aᵢ)))
    elseif Meta.isexpr(aᵢ, :$, 1)
      push!(array.args, Expr(:call, GlobalRef(StrideArrays, :StaticInt), (only(aᵢ.args))))
    else
      push!(array.args, esc(aᵢ))
    end
  end
  f! = Symbol(expr.args[1], :!)
  # we call PtrArray so that `rand!` doesn't force heap allocation.
  call = Expr(:call, f!, Expr(:call, :PtrArray, :array))
  for arg ∈ args
    push!(call.args, esc(arg))
  end
  quote
    array = $array
    b = array.data
    # we add GC.@preserve to save the array.
    GC.@preserve b $call
    array
  end
end

"""
    @StrideArray rand(Float32, 3, 4)
    @StrideArray randn($(7>>1),  4)
    @StrideArray rand(7>>1,  4)

Creates a random `StrideArray`.
The default element type is `Float64`.
Dimensions will be statically sized if specified by an integer literal, or if interpolated.
"""
macro StrideArray(expr, args...)
  rand_expr(expr, args...)
end


text:
# metro_params.jl
# Select parameters and save them into JSON files
# Authored by Arthur Delarue on 10/28/15

using JSON

PARAMS = Dict(
"PROB" => 0.7,						# determines number of nodes with data (probability of no data for a given node pair)

"MODEL" => "relaxed_smooth",			# relaxed/strict _ smooth/nothing/random
"LAST_SMOOTH" => false, 				# true if last iteration is smooth, false otherwise. Ignore if model is smooth already

"MAX_ROUNDS" => 50,						# max number of iterations

"MIN_RIDES" => 1,						# min number of rides

"DYNAMIC_CONSTRAINTS" => true, 			# true for dynamic constraints, false otw
"GLOBAL_CONSTRAINT_UPDATE" => true,		# true if we remove worst paths globally, false otw

"SAMPLE_SIZE" => 1000, 					# starting number of constraints
"NUM_OD_ADDED" => 2500,					# number of (O,D) pairs to add
"NUM_OD_REMOVED" => 1000,				# number of (O,D) pairs to remove
"UPDATE_EVERY_N_ITERATIONS" => 2, 		# add number of (O,D) above every $N iterations
"ADD_IF_ABOVE_PERCENTILE" => 0.7,		# only add rides that have an error above this percentile
"REMOVE_IF_BELOW_PERCENTILE" => 0.3,	# only remove rides that have an error below this percentile

"TURN_COST" => 2.0,	 					# turning cost initial value
"TURN_COST_AS_VARIABLE" => false, 		# true if LP updates turning cost, false otw

"DELTA_BOUND" => [0.6, 0.5,0.4],		# relaxation base (sketchy)

"MAX_NUM_PATHS_PER_OD" => 3, 			# at least 1 please

"START_SIMPLE" => false 				# true if initial simple LP is used, false otherwise
)

function create_JSON_file(params::Dict{ASCIIString,Any}, fileName::AbstractString)
	file = open(fileName, "w")
	write(file, json(params))
	close(file)
end

function add_JSON_file(params::Dict{ASCIIString,Any})
	i = 0
	while isfile("JSONparams/params-$i.json")
		i += 1
	end
	create_JSON_file(params, "JSONparams/params-$i.json")
end

add_JSON_file(PARAMS)

max_stars_repo_name: sebmart/NYC-timing

max_issues_repo_name: sebmart/NYC-timing

max_forks_repo_name: sebmart/NYC-timing


text:
@testset "save and load" begin
    a = Variable(1.0)
    b = Variable(2.0)
    init(sess)
    save(sess, "temp.mat", [a,b])
    op1 = assign(a, 3.0)
    op2 = assign(b, 3.0)
    run(sess, op1)
    run(sess, op2)

    sess1 = Session()
    load(sess1, "temp.mat", [a,b])
    @test run(sess1, [a,b])≈[1.0;2.0]
    rm("temp.mat", force=true)
end

@testset "psave and pload" begin
py"""
def psave_test_function(x):
    return x+1
"""
    a = py"psave_test_function"
    psave(a, "temp.pkl")
    b = pload("temp.pkl")
    @test b(2)==3
    rm("temp.pkl", force=true)
end

@testset "diary" begin
    d = Diary("test")
    b = constant(1.0)
    a = scalar(b)
    for i = 1:100
        write(d, i, run(sess, a))
    end
    # activate(d)
    save(d, "mydiary")
    e = Diary("test2")
    load(e, "mydiary")
    # activate(e)

    try;rm("mydiary", recursive=true, force=true);catch;end
end


text:
id_me = myid()
id_other = addprocs(1)[1]

using MessageUtils
using Base.Test

macro dbgtest(ex)
    println("Testing $ex")
    @eval @test $ex
end

include("test_channel.jl")
include("test_ts.jl")
include("test_kvs.jl")
include("test_ct.jl")


text:
# module Ecto.Changeset

module Relation

import ....Ecto

function primary_keys!(schema::Ecto.Schema.t)
    pks = schema[:primary_key]
    isempty(pks) ? throw(Ecto.NoPrimaryKeyFieldError(schema)) : pks
end

end # Ecto.Changeset.Relation


text:
using Luxor, Test

function crescent_tests(fname)
    Drawing(800, 800, fname)
    origin()
    background("midnightblue")

    tiles = Tiler(800, 800, 6, 6)
    for (pos, n) in tiles
        pts = crescent(pos, rescale(n, 36, 1, -tiles.tilewidth/2, tiles.tilewidth/2), tiles.tilewidth/2, vertices=true)
        sethue("antiquewhite")
        poly(pts, :fill)
        sethue("purple")
        crescent(pos, 15 , pos + (20, 0), 20, :stroke)
    end

    # test for null crescent, no intersections
    pts = crescent(O, 10, O, 10, vertices=true)
    @test pts == nothing
    @test finish() == true
    println("...finished test: output in $(fname)")
end

crescent_tests("crescent-tests.svg")


text:
struct ProjCGWork{T}
    r::Vector{T}
    g::Vector{T}
    d::Vector{T}
    rp::Vector{T}
    gp::Vector{T}
    Ad::Vector{T}
    Utr::Vector{T}
end

ProjCGWork(n::Int, m::Int) = ProjCGWork([Array{Float64}(undef, n) for i in 1:6]..., Array{Float64}(undef, m))

#=
PROJCG does projected conjugate gradient with orthonormal projection matrix
U and PSD matrix A. That is, the following system is solved:

[ A  U ] [ x ] = [ b ]
[ U' 0 ] [ λ ]   [ c ]

Returns a direction of negative curvature of A if found, storing the direction
in x and setting λ to all NaNs.

Usage:
  [x, lambda, resvec, Uvec, dvec] = projcg(A, b, U, c, tol, maxit)

Inputs:
A, b, U, c - as defined above
tol [1e-6] - absolute tolerance for convergence of 2-norm of residual on the
  constraint manifold
maxit [size(A,1) - size(U,2)] - maximum number of iterations if tol is not met

Overwrites:
x, λ - solution

Output:
i - number of iterations taken
nr - final tolerance achieved for projected residual, (I - U*U')*(A*x - b)
=#
function projcg!(x, λ, A, U, b, c;
    tol::Float64=1e-6, maxit::Int=length(b)+length(c), work::ProjCGWork{T}=ProjCGWork(length(b), length(c))) where T

n = length(b)
m = length(c)

# work variables
r = view(work.r, 1:n)
g = view(work.g, 1:n)
d = view(work.d, 1:n)
rp = view(work.rp, 1:n)
gp = view(work.gp, 1:n)
Ad = view(work.Ad, 1:n)
Utr = view(work.Utr, 1:m)

mul!(x, U, c)                       # satisfies U' x = c
r .= b
mul!(r, A, x, 1.0, -1.0)            # r = Ax - b
g .= r
mul!(Utr, U', r)
mul!(g, U, Utr, -1.0, 1.0)          # g = r - U*(U'*r);
r .= g
d .= -1.0.*g                        # multiplication avoids allocation

nr0 = norm(g);

# resvec = [nr0; zeros(maxit, 1)];

i = 0
nr = Inf

while i < min(maxit, n+m)
    i += 1

    mul!(Ad, A, d)                  # Ad = A*d
    dAd = dot(d, Ad);

    if dAd <= 0
        # resvec = resvec(1:i);
        x .= d / norm(d)
        λ .= NaN
        return i, Inf
    end

    rg = dot(r, g);

    # rg cannot be negative
    if rg <= 0
        break
    end

    α = rg / dAd
    x .+= α.*d
    rp .= r .+ α.*Ad

    gp .= rp
    mul!(Utr, U', rp)
    mul!(gp, U, Utr, -1.0, 1.0)      # gp = rp - U*(U'*rp);
    β = dot(rp, gp) / rg
    d .= β.*d .- gp
    g .= gp
    r .= gp

    nr = norm(g);
    # resvec(i+1) = nr;

    if nr < tol
        # resvec = resvec(1:i+1);

        break
    end
end

# solve for Lagrange multipliers
r .= b
mul!(r, A, x, -1.0, 1.0)            # r = b - A*x

mul!(λ, U', r);          # λ = U'*(b - A*x)

return i, nr
end


text:
using Revise

using RCall

const N = Int(1e7); const K = 100;

R"""
memory.limit(2e9)
library(data.table)
N=$N; K=$K
set.seed(1)
DT <- data.table(
  id1 = sample(sprintf("id%03d",1:K), N, TRUE),      # large groups (char)
  v3 =  sample(round(runif(100,max=100),4), N, TRUE) # numeric e.g. 23.5749
)
cat("GB =", round(sum(gc()[,2])/1024, 3), "\n")
pt = proc.time()
DT[, sum(v3), keyby=id1]
x = timetaken(pt)
rm(DT); gc()
"""

@rget x;

using FastGroupBy, BenchmarkTools, DataFrames

srand(1);
df = DataFrame(id_str_small = rand(["id"*dec(k,3) for k = 1:K],N) , val = rand(round.(rand(K)*100,4), N));
x1 = @elapsed sumby!(df, :id_str_small, :val);

x1/parse(Float64, x[1:end-3])


text:
"""
Takes a time series and the noise envelope, and returns the values of the time
series that are "significant" i.e. outside the noise envelope
"""
function process_residual(mean, std)
    low = mean .- std
    high = mean .+ std
    real_ind = .!((low .< 0.0) .& (high .> 0.0))

    ctr = zeros(size(mean))
    ctr[real_ind] = mean[real_ind]

    return ctr
end

export process_residual


text:
struct Circles
    n::Int64
end
Circles() = Circles(10)

function draw_item(circles::Circles, canvas::Gtk.GtkCanvas, scene::Scene)
    ctx = Gtk.getgc(canvas)
    h = height(canvas)
    w = width(canvas)

    ti = scene.tNow
    for k=1:circles.n
        x = (cos(ti*3) * ((w/2.0)-30.0)) + (w/2.0) +5.0
        y = (sin(ti*4) * ((h/2.0)-40.0)) + (h/2.0) +15.0

        new_path(ctx)
        set_line_width(ctx, 10.0);
        set_source_rgb(ctx, 0.5, mod(ti,1.0), 0.7)
        arc(ctx, x, y, 20.0, 0.0, 6.28);
        stroke(ctx);
        ti = ti + 0.05
    end
    return nothing
end


text:
struct Cubic{BC<:BoundaryCondition} <: DegreeBC{3}
    bc::BC
end

(deg::Cubic)(gt::GridType) = Cubic(deg.bc(gt))

"""
    Cubic(bc::BoundaryCondition)

Indicate that the corresponding axis should use quadratic interpolation.

# Extended help

Assuming uniform knots with spacing 1, the `i`th piece of cubic spline
implemented here is defined as follows.

    y_i(x) = cm p(x-i) + c q(x-i) + cp q(1- (x-i)) + cpp p(1 - (x-i))

where

    p(δx) = 1/6 * (1-δx)^3
    q(δx) = 2/3 - δx^2 + 1/2 δx^3

and the values `cX` for `X ∈ {m, _, p, pp}` are the pre-filtered coefficients.

For future reference, this expands out to the following polynomial:

    y_i(x) = 1/6 cm (1+i-x)^3 + c (2/3 - (x-i)^2 + 1/2 (x-i)^3) +
             cp (2/3 - (1+i-x)^2 + 1/2 (1+i-x)^3) + 1/6 cpp (x-i)^3

When we derive boundary conditions we will use derivatives `y_0'(x)` and
`y_0''(x)`
"""
Cubic

function positions(deg::Cubic, ax, x)
    xf = floorbounds(x, ax)
    xf -= ifelse(xf > last(ax)-1, oneunit(xf), zero(xf))
    δx = x - xf
    expand_index(deg, fast_trunc(Int, xf), ax, δx), δx
end

expand_index(::Cubic{BC}, xi::Number, ax::AbstractUnitRange, δx) where BC = xi-1
expand_index(::Cubic{Periodic{GT}}, xi::Number, ax::AbstractUnitRange, δx) where GT<:GridType =
    (modrange(xi-1, ax), modrange(xi, ax), modrange(xi+1, ax), modrange(xi+2, ax))

function value_weights(::Cubic, δx)
    x3, xcomp3 = cub(δx), cub(1-δx)
    (SimpleRatio(1,6) * xcomp3,
     SimpleRatio(2,3) - sqr(δx) + SimpleRatio(1,2)*x3,
     SimpleRatio(2,3) - sqr(1-δx) + SimpleRatio(1,2)*xcomp3,
     SimpleRatio(1,6) * x3)
end

function gradient_weights(::Cubic, δx)
    x2, xcomp2 = sqr(δx), sqr(1-δx)
    (-SimpleRatio(1,2) * xcomp2,
     -2*δx + SimpleRatio(3,2)*x2,
     +2*(1-δx) - SimpleRatio(3,2)*xcomp2,
     SimpleRatio(1,2) * x2)
end

hessian_weights(::Cubic, δx) = (1-δx, 3*δx-2, 3*(1-δx)-2, δx)

# ------------ #
# Prefiltering #
# ------------ #

padded_axis(ax::AbstractUnitRange, ::BSpline{<:Cubic}) = first(ax)-1:last(ax)+1
padded_axis(ax::AbstractUnitRange, ::BSpline{Cubic{Periodic{GT}}}) where GT<:GridType = ax

# # Due to padding we can extend the bounds
# lbound(ax, ::BSpline{Cubic{BC}}, ::OnGrid) where BC = first(ax) - 0.5
# ubound(ax, ::BSpline{Cubic{BC}}, ::OnGrid) where BC = last(ax) + 0.5

"""
`Cubic`: continuity in function value, first and second derivatives yields

    2/3 1/6
    1/6 2/3 1/6
        1/6 2/3 1/6
           ⋱  ⋱   ⋱
"""
function inner_system_diags(::Type{T}, n::Int, ::Cubic) where {T}
    du = fill(convert(T, SimpleRatio(1, 6)), n-1)
    d = fill(convert(T, SimpleRatio(2, 3)), n)
    dl = copy(du)
    dl, d, du
end

"""
`Cubic{Flat}` `OnGrid` amounts to setting `y_1'(x) = 0` at `x = 1`.
Applying this condition yields

    -cm + cp = 0
"""
function prefiltering_system(::Type{T}, ::Type{TC}, n::Int,
                             degree::Cubic{Flat{OnGrid}}) where {T,TC}
    dl, d, du = inner_system_diags(T, n, degree)
    d[1] = d[end] = -oneunit(T)
    du[1] = dl[end] = zero(T)

    # Now Woodbury correction to set `[1, 3], [n, n-2] ==> 1`
    specs = WoodburyMatrices.sparse_factors(T, n, (1, 3, oneunit(T)), (n, n-2, oneunit(T)))

    Woodbury(lut!(dl, d, du), specs...), zeros(TC, n)
end

"""
`Cubic{Flat}`, `OnCell` amounts to setting `y_1'(x) = 0` at `x = 1/2`.
Applying this condition yields

    -9/8 cm + 11/8 c - 3/8 cp + 1/8 cpp = 0

or, equivalently,

    -9 cm + 11 c -3 cp + 1 cpp = 0

(Note that we use `y_1'(x)` although it is strictly not valid in this domain; if we
were to use `y_0'(x)` we would have to introduce new coefficients, so that would not
close the system. Instead, we extend the outermost polynomial for an extra half-cell.)
"""
function prefiltering_system(::Type{T}, ::Type{TC}, n::Int,
                             degree::Cubic{Flat{OnCell}}) where {T,TC}
    dl, d, du = inner_system_diags(T,n,degree)
    d[1] = d[end] = -9
    du[1] = dl[end] = 11

    # now need Woodbury correction to set :
    #    - [1, 3] and [n, n-2] ==> -3
    #    - [1, 4] and [n, n-3] ==> 1
    specs = WoodburyMatrices.sparse_factors(T, n,
                                  (1, 3, T(-3)),
                                  (n, n-2, T(-3)),
                                  (1, 4, oneunit(T)),
                                  (n, n-3, oneunit(T))
                                  )

    Woodbury(lut!(dl, d, du), specs...), zeros(TC, n)
end

"""
`Cubic{Line}` `OnCell` amounts to setting `y_1''(x) = 0` at `x = 1/2`.
Applying this condition yields

    3 cm -7 c + 5 cp -1 cpp = 0

(Note that we use `y_1'(x)` although it is strictly not valid in this domain; if we
were to use `y_0'(x)` we would have to introduce new coefficients, so that would not
close the system. Instead, we extend the outermost polynomial for an extra half-cell.)
"""
function prefiltering_system(::Type{T}, ::Type{TC}, n::Int,
                             degree::Cubic{Line{OnCell}}) where {T,TC}
    dl,d,du = inner_system_diags(T,n,degree)
    d[1] = d[end] = 3
    du[1] = dl[end] = -7

    # now need Woodbury correction to set :
    #    - [1, 3] and [n, n-2] ==> -3
    #    - [1, 4] and [n, n-3] ==> 1
    specs = WoodburyMatrices.sparse_factors(T, n,
                                  (1, 3, T(5)),
                                  (n, n-2, T(5)),
                                  (1, 4, -oneunit(T)),
                                  (n, n-3, -oneunit(T))
                                  )

    Woodbury(lut!(dl, d, du), specs...), zeros(TC, n)
end

"""
`Cubic{Line}` `OnGrid` amounts to setting `y_1''(x) = 0` at `x = 1`. Applying this
condition gives:

    1 cm -2 c + 1 cp = 0
"""
function prefiltering_system(::Type{T}, ::Type{TC}, n::Int,
                             degree::Cubic{Line{OnGrid}}) where {T,TC}
    dl,d,du = inner_system_diags(T,n,degree)
    d[1] = d[end] = 1
    du[1] = dl[end] = -2

    # now need Woodbury correction to set :
    #    - [1, 3] and [n, n-2] ==> 1
    specs = WoodburyMatrices.sparse_factors(T, n,
                                  (1, 3, oneunit(T)),
                                  (n, n-2, oneunit(T)),
                                  )

    Woodbury(lut!(dl, d, du), specs...), zeros(TC, n)
end

"""
`Cubic{Periodic}` `OnGrid` closes the system by looking at the coefficients themselves
as periodic, yielding

    c0 = c(N+1)

where `N` is the number of data points.
"""
function prefiltering_system(::Type{T}, ::Type{TC}, n::Int,
                             degree::Cubic{<:Periodic}) where {T,TC}
    dl, d, du = inner_system_diags(T,n,degree)

    specs = WoodburyMatrices.sparse_factors(T, n,
                                  (1, n, du[1]),
                                  (n, 1, dl[end])
                                  )

    Woodbury(lut!(dl, d, du), specs...), zeros(TC, n)
end

"""
`Cubic{Free}` `OnGrid` and `Cubic{Free}` `OnCell` amount to requiring an extra
continuous derivative at the second-to-last cell boundary; this means
`y_1'''(2) = y_2'''(2)`, yielding

    1 cm -3 c + 3 cp -1 cpp = 0
"""
function prefiltering_system(::Type{T}, ::Type{TC}, n::Int,
                             degree::Cubic{<:Free}) where {T,TC}
    dl, d, du = inner_system_diags(T,n,degree)

    specs = WoodburyMatrices.sparse_factors(T, n,
                                  (1, n, du[1]),
                                  (n, 1, dl[end])
                                  )

    Woodbury(lut!(dl, d, du), specs...), zeros(TC, n)
end


text:
"""
```julia
    LogisticSoftMaxLikelihood()
```

The multiclass likelihood with a logistic-softmax mapping: :
```math
p(y=i|{fₖ}₁ᴷ) = σ(fᵢ)/∑ₖ σ(fₖ)
```
where `σ` is the logistic function.
This likelihood has the same properties as [softmax](https://en.wikipedia.org/wiki/Softmax_function).
---

For the analytical version, the likelihood is augmented multiple times. More details can be found in the paper [Multi-Class Gaussian Process Classification Made Conjugate: Efficient Inference via Data Augmentation](https://arxiv.org/abs/1905.09670)
"""
struct LogisticSoftMaxLikelihood{T<:Real} <: MultiClassLikelihood{T}
    Y::AbstractVector{BitVector} #Mapping from instances to classes
    class_mapping::AbstractVector{Any} # Classes labels mapping
    ind_mapping::Dict{Any,Int} # Mapping from label to index
    y_class::AbstractVector{Int64} # GP Index for each sample
    c::AbstractVector{AbstractVector{T}} # Second moment of fₖ
    α::AbstractVector{T} # Variational parameter of Gamma distribution
    β::AbstractVector{T} # Variational parameter of Gamma distribution
    θ::AbstractVector{AbstractVector{T}} # Variational parameter of Polya-Gamma distribution
    γ::AbstractVector{AbstractVector{T}} # Variational parameter of Poisson distribution
    function LogisticSoftMaxLikelihood{T}() where {T<:Real}
        new{T}()
    end
    function LogisticSoftMaxLikelihood{T}(Y::AbstractVector{<:BitVector},
    class_mapping::AbstractVector, ind_mapping::Dict{<:Any,<:Int},y_class::AbstractVector{<:Int}) where {T<:Real}
        new{T}(Y,class_mapping,ind_mapping,y_class)
    end
    function LogisticSoftMaxLikelihood{T}(Y::AbstractVector{<:BitVector},
    class_mapping::AbstractVector, ind_mapping::Dict{<:Any,<:Int},y_class::AbstractVector{<:Int},
    c::AbstractVector{<:AbstractVector{<:Real}}, α::AbstractVector{<:Real},
    β::AbstractVector, θ::AbstractVector{<:AbstractVector},γ::AbstractVector{<:AbstractVector{<:Real}}) where {T<:Real}
        new{T}(Y,class_mapping,ind_mapping,y_class,c,α,β,θ,γ)
    end
end

function LogisticSoftMaxLikelihood()
    LogisticSoftMaxLikelihood{Float64}()
end

function pdf(l::LogisticSoftMaxLikelihood,f::AbstractVector)
    logisticsoftmax(f)
end

function pdf(l::LogisticSoftMaxLikelihood,y::Integer,f::AbstractVector)
    logisticsoftmax(f)[y]
end

function Base.show(io::IO,model::LogisticSoftMaxLikelihood{T}) where T
    print(io,"Logistic-Softmax Likelihood")
end

function init_likelihood(likelihood::LogisticSoftMaxLikelihood{T},inference::Inference{T},nLatent::Integer,nSamplesUsed::Integer,nFeatures::Integer) where T
    if inference isa AnalyticVI || inference isa GibbsSampling
        c = [ones(T,nSamplesUsed) for i in 1:nLatent]
        α = nLatent*ones(T,nSamplesUsed)
        β = nLatent*ones(T,nSamplesUsed)
        θ = [abs.(rand(T,nSamplesUsed))*2 for i in 1:nLatent]
        γ = [abs.(rand(T,nSamplesUsed)) for i in 1:nLatent]
        LogisticSoftMaxLikelihood{T}(likelihood.Y,likelihood.class_mapping,likelihood.ind_mapping,likelihood.y_class,c,α,β,θ,γ)
    else
        return likelihood
    end
end

## Local Updates##
function local_updates!(l::LogisticSoftMaxLikelihood,y,μ::NTuple{N,<:AbstractVector},Σ::NTuple{N,<:AbstractVector}) where {T,N}
    l.c .= broadcast((Σ,μ)->sqrt.(Σ+abs2.(μ)),Σ,μ)
    for _ in 1:2
        l.γ .= broadcast(
            (β,c,μ,ψα)->0.5 / β * exp.(ψα) .* safe_expcosh.(-0.5*μ, 0.5*c),                                    [l.β],l.c,μ,[digamma.(l.α)])
        l.α .= 1.0.+(l.γ...)
    end
    l.θ .= broadcast((y,γ,c)->0.5*(y+γ)./c.*tanh.(0.5.*c),
                                    y,l.γ,l.c)
    return nothing
end

function sample_local!(l::LogisticSoftMaxLikelihood{T},y::AbstractVector,f) where {T}
    l.γ .= broadcast(f->rand.(Poisson.(0.5*l.α.*safe_expcosh.(-0.5*f,0.5*f))), f)
    l.α .= rand.(Gamma.(one(T).+(l.γ...),1.0./l.β))
    pg = PolyaGammaDist()
    set_ω!(l,broadcast((y,γ,f)->draw.([pg],y.+γ,μ),y,l.γ,f))
    return nothing
end

## Global Gradient Section ##

@inline ∇E_μ(l::LogisticSoftMaxLikelihood,::AOptimizer,y::AbstractVector) where {T} = 0.5.*(y.-l.γ)
@inline ∇E_Σ(l::LogisticSoftMaxLikelihood,::AOptimizer,y::AbstractVector) where {T} = 0.5.*l.θ

## ELBO Section ##
function expec_log_likelihood(l::LogisticSoftMaxLikelihood{T},i::AnalyticVI,y,μ,Σ) where {T}
    tot = -length(y)*logtwo
    tot += -sum(sum(l.γ.+y))*logtwo
    tot +=  0.5*sum(broadcast((θ,γ,y,μ,Σ)->dot(μ,(y-γ))-dot(θ,abs2.(μ))-dot(θ,Σ),l.θ,l.γ,y,μ,Σ))
    return tot
end

AugmentedKL(l::LogisticSoftMaxLikelihood,y::AbstractVector) = PolyaGammaKL(l,y) + PoissonKL(l) + GammaEntropy(l)

PolyaGammaKL(l::LogisticSoftMaxLikelihood,y) =  sum(broadcast(PolyaGammaKL,y.+l.γ,l.c,l.θ))

PoissonKL(l::LogisticSoftMaxLikelihood) =     sum(broadcast(PoissonKL,l.γ,[l.α./l.β],[digamma.(l.α).-log.(l.β)]))

##  Compute the equivalent of KL divergence between an improper prior p(λ) (``1_{[0,\\infty]}``) and a variational Gamma distribution ##
function GammaEntropy(l::LogisticSoftMaxLikelihood)
    return -sum(l.α)+sum(log,l.β[1])-sum(lgamma,l.α)-dot(1.0.-l.α,digamma.(l.α))
end

## Numerical Gradient Section ##

function grad_samples(model::AbstractGP{T,<:LogisticSoftMaxLikelihood,<:NumericalVI},samples::AbstractMatrix{T},index::Integer) where {T}
    class = model.likelihood.y_class[index]::Int64
    grad_μ = zeros(T,model.nLatent)
    grad_Σ = zeros(T,model.nLatent)
    g_μ = similar(grad_μ)
    nSamples = size(samples,1)
    @views @inbounds for i in 1:nSamples
        σ = logistic.(samples[i,:])
        samples[i,:]  .= logisticsoftmax(samples[i,:])
        s = samples[i,class]
        g_μ .= grad_logisticsoftmax(samples[i,:],σ,class)/s
        grad_μ += g_μ
        grad_Σ += diaghessian_logisticsoftmax(samples[i,:],σ,class)/s - abs2.(g_μ)
    end
    for k in 1:model.nLatent
        model.inference.ν[k][index] = -grad_μ[k]/nSamples
        model.inference.λ[k][index] = grad_Σ[k]/nSamples
    end
end

function log_like_samples(model::AbstractGP{T,<:LogisticSoftMaxLikelihood},samples::AbstractMatrix,index::Integer) where {T}
    class = model.likelihood.y_class[index]
    nSamples = size(samples,1)
    loglike = zero(T)
    for i in 1:nSamples
        σ = logistic.(samples[i,:])
        loglike += log(σ[class])-log(sum(σ))
    end
    return loglike/nSamples
end

function grad_logisticsoftmax(s::AbstractVector{<:Real},σ::AbstractVector{<:Real},i::Integer)
    s[i]*(δ.(i,eachindex(σ)).-s).*(1.0.-σ)
end

function diaghessian_logisticsoftmax(s::AbstractVector{<:Real},σ::AbstractVector{<:Real},i::Integer)
    s[i]*(1.0.-σ).*(
    abs2.(δ.(i,eachindex(σ))-s).*(1.0.-σ)
    -s.*(1.0.-s).*(1.0.-σ)
    -σ.*(δ.(i,eachindex(σ))-s))
end

function hessian_logisticsoftmax(s::AbstractVector{T},σ::AbstractVector{T},i::Integer) where {T<:Real}
    m = length(s)
    hessian = zeros(T,m,m)
    @inbounds for j in 1:m
        for k in 1:m
            hessian[j,k] = (1-σ[j])*s[i]*(
            (δ(i,k)-s[k])*(1.0-σ[k])*(δ(i,j)-s[j])
            -s[j]*(δ(j,k)-s[k])*(1.0-σ[k])
            -δ(k,j)*σ[j]*(δ(i,j)-s[j]))
        end
    end
    return hessian
end


text:
# Note: we define type aliases outside of the @muladd block to avoid Revise breaking when code
# inside the @muladd block is edited. See https://github.com/trixi-framework/Trixi.jl/issues/801
# for more details.

# `DGMulti` refers to both multiple DG types (polynomial/SBP, simplices/quads/hexes) as well as
# the use of multi-dimensional operators in the solver.
const DGMulti{NDIMS, ElemType, ApproxType, SurfaceIntegral, VolumeIntegral} =
  DG{<:RefElemData{NDIMS, ElemType, ApproxType}, Mortar, SurfaceIntegral, VolumeIntegral} where {Mortar}

# Type aliases. The first parameter is `ApproxType` since it is more commonly used for dispatch.
const DGMultiWeakForm{ApproxType, ElemType} =
  DGMulti{NDIMS, ElemType, ApproxType, <:SurfaceIntegralWeakForm, <:VolumeIntegralWeakForm} where {NDIMS}

const DGMultiFluxDiff{ApproxType, ElemType} =
  DGMulti{NDIMS, ElemType, ApproxType, <:SurfaceIntegralWeakForm, <:VolumeIntegralFluxDifferencing} where {NDIMS}

const DGMultiFluxDiffSBP{ApproxType, ElemType} =
  DGMulti{NDIMS, ElemType, ApproxType, <:SurfaceIntegralWeakForm, <:VolumeIntegralFluxDifferencing} where {NDIMS, ApproxType<:Union{SBP, AbstractDerivativeOperator}}

const DGMultiSBP{ApproxType, ElemType} =
  DGMulti{NDIMS, ElemType, ApproxType, SurfaceIntegral, VolumeIntegral} where {NDIMS, ElemType, ApproxType<:Union{SBP, AbstractDerivativeOperator}, SurfaceIntegral, VolumeIntegral}

# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

# these are necessary for pretty printing
polydeg(dg::DGMulti) = dg.basis.N
Base.summary(io::IO, dg::DG) where {DG <: DGMulti} = print(io, "DGMulti(polydeg=$(polydeg(dg)))")
Base.real(rd::RefElemData{NDIMS, Elem, ApproxType, Nfaces, RealT}) where {NDIMS, Elem, ApproxType, Nfaces, RealT} = RealT

"""
    DGMulti(; polydeg::Integer,
              element_type::AbstractElemShape,
              approximation_type=Polynomial(),
              surface_flux=flux_central,
              surface_integral=SurfaceIntegralWeakForm(surface_flux),
              volume_integral=VolumeIntegralWeakForm(),
              RefElemData_kwargs...)

Create a discontinuous Galerkin method which uses
- approximations of polynomial degree `polydeg`
- element type `element_type` (`Tri()`, `Quad()`, `Tet()`, and `Hex()` currently supported)

Optional:
- `approximation_type` (default is `Polynomial()`; `SBP()` also supported for `Tri()`, `Quad()`,
  and `Hex()` element types).
- `RefElemData_kwargs` are additional keyword arguments for `RefElemData`, such as `quad_rule_vol`.
  For more info, see the [StartUpDG.jl docs](https://jlchan.github.io/StartUpDG.jl/dev/).
"""
function DGMulti(; polydeg=nothing,
                   element_type::AbstractElemShape,
                   approximation_type=Polynomial(),
                   surface_flux=flux_central,
                   surface_integral=SurfaceIntegralWeakForm(surface_flux),
                   volume_integral=VolumeIntegralWeakForm(),
                   kwargs...)

  # call dispatchable constructor
  DGMulti(element_type, approximation_type, volume_integral, surface_integral;
          polydeg=polydeg, kwargs...)
end

# dispatchable constructor for DGMulti to allow for specialization
function DGMulti(element_type::AbstractElemShape,
                 approximation_type,
                 volume_integral,
                 surface_integral;
                 polydeg::Integer,
                 kwargs...)

  rd = RefElemData(element_type, approximation_type, polydeg; kwargs...)
  return DG(rd, nothing #= mortar =#, surface_integral, volume_integral)
end

########################################
#            DGMultiMesh
########################################

# now that `DGMulti` is defined, we can define constructors for `DGMultiMesh` which use `dg::DGMulti`

function DGMultiMesh(dg::DGMulti, geometric_term_type, md::MeshData{NDIMS}, boundary_faces) where {NDIMS}
  return DGMultiMesh{NDIMS, geometric_term_type, typeof(md), typeof(boundary_faces)}(md, boundary_faces)
end

# TODO: DGMulti, v0.5. These constructors which use `rd::RefElemData` are now redundant and can be removed.
function DGMultiMesh(vertex_coordinates::NTuple{NDIMS, Vector{Tv}}, EToV::Array{Ti,2}, rd::RefElemData;
                     is_on_boundary = nothing,
                     periodicity=ntuple(_->false, NDIMS), kwargs...) where {NDIMS, Tv, Ti}

  Base.depwarn("`DGMultiMesh` constructor with `rd::RefElemData` is deprecated. Use the constructor with `dg::DGMulti` instead.",
               :DGMultiMesh)
  if haskey(kwargs, :is_periodic)
    # TODO: DGMulti, v0.5. Remove deprecated keyword
    Base.depwarn("keyword argument `is_periodic` is now `periodicity`.", :DGMultiMesh)
    periodicity=kwargs[:is_periodic]
  end

  md = MeshData(vertex_coordinates, EToV, rd)
  if NDIMS==1
    md = StartUpDG.make_periodic(md, periodicity...)
  else
    md = StartUpDG.make_periodic(md, periodicity)
  end
  boundary_faces = StartUpDG.tag_boundary_faces(md, is_on_boundary)
  return DGMultiMesh{NDIMS, typeof(rd.elementType), typeof(md), typeof(boundary_faces)}(md, boundary_faces)
end

function DGMultiMesh(triangulateIO, rd::RefElemData{2, Tri}, boundary_dict::Dict{Symbol, Int})

  vertex_coordinates, EToV = StartUpDG.triangulateIO_to_VXYEToV(triangulateIO)
  md = MeshData(vertex_coordinates, EToV, rd)
  boundary_faces = StartUpDG.tag_boundary_faces(triangulateIO, rd, md, boundary_dict)
  return DGMultiMesh{2, typeof(rd.elementType), typeof(md), typeof(boundary_faces)}(md, boundary_faces)
end

# TODO: DGMulti, v0.5. Remove deprecated constructor
@deprecate VertexMappedMesh(args...; kwargs...) DGMultiMesh(args...; kwargs...)

# Mesh types used internally for trait dispatch
struct Cartesian end
struct VertexMapped end # where element geometry is determined by vertices.
struct Curved end

# type parameters for dispatch using `DGMultiMesh`
abstract type GeometricTermsType end
struct Affine <: GeometricTermsType end # mesh produces constant geometric terms
struct NonAffine <: GeometricTermsType end # mesh produces non-constant geometric terms

# choose MeshType based on the constructor and element type
GeometricTermsType(mesh_type, dg::DGMulti) = GeometricTermsType(mesh_type, dg.basis.elementType)
GeometricTermsType(mesh_type::Cartesian, element_type::AbstractElemShape) = Affine()
GeometricTermsType(mesh_type::TriangulateIO, element_type::Tri) = Affine()
GeometricTermsType(mesh_type::VertexMapped, element_type::Union{Tri, Tet}) = Affine()
GeometricTermsType(mesh_type::VertexMapped, element_type::Union{Quad, Hex}) = NonAffine()
GeometricTermsType(mesh_type::Curved, element_type::AbstractElemShape) = NonAffine()

# other potential constructor types to add later: Bilinear, Isoparametric{polydeg_geo}, Rational/Exact?
# other potential mesh types to add later: Polynomial{polydeg_geo}?

"""
  DGMultiMesh(vertex_coordinates, EToV, dg::DGMulti{NDIMS};
              is_on_boundary=nothing,
              periodicity=ntuple(_->false, NDIMS)) where {NDIMS, Tv}

- `vertex_coordinates` is a tuple of vectors containing x,y,... components of the vertex coordinates
- `EToV` is a 2D array containing element-to-vertex connectivities for each element
- `dg::DGMulti` contains information associated with to the reference element (e.g., quadrature,
  basis evaluation, differentiation, etc).
- `is_on_boundary` specifies boundary using a `Dict{Symbol, <:Function}`
- `periodicity` is a tuple of booleans specifying if the domain is periodic `true`/`false` in the
  (x,y,z) direction.
"""
# TODO: DGMulti v0.5. Standardize order of arguments, pass in `dg` first
function DGMultiMesh(vertex_coordinates, EToV, dg::DGMulti{NDIMS};
                     is_on_boundary=nothing,
                     periodicity=ntuple(_->false, NDIMS), kwargs...) where {NDIMS}
  if haskey(kwargs, :is_periodic)
    # TODO: DGMulti, v0.5. Remove deprecated keyword
    Base.depwarn("keyword argument `is_periodic` is now `periodicity`.", :DGMultiMesh)
    periodicity=kwargs[:is_periodic]
  end

  md = MeshData(vertex_coordinates, EToV, dg.basis)
  if NDIMS == 1
    md = StartUpDG.make_periodic(md, periodicity...)
  else
    md = StartUpDG.make_periodic(md, periodicity)
  end
  boundary_faces = StartUpDG.tag_boundary_faces(md, is_on_boundary)
  return DGMultiMesh(dg, GeometricTermsType(VertexMapped(), dg), md, boundary_faces)
end

"""
    DGMultiMesh(triangulateIO, dg::DGMulti{2, Tri}, boundary_dict::Dict{Symbol, Int})

- `triangulateIO` is a `TriangulateIO` mesh representation
- `dg::DGMulti` contains information associated with to the reference element (e.g., quadrature,
  basis evaluation, differentiation, etc).
- `boundary_dict` is a `Dict{Symbol, Int}` which associates each integer `TriangulateIO` boundary
  tag with a `Symbol`.
"""
# TODO: DGMulti v0.5, standardize order of arguments (`dg` first)
function DGMultiMesh(triangulateIO, dg::DGMulti{2, Tri}, boundary_dict::Dict{Symbol, Int};
                     periodicity=(false, false))
  vertex_coordinates, EToV = StartUpDG.triangulateIO_to_VXYEToV(triangulateIO)
  md = MeshData(vertex_coordinates, EToV, dg.basis)
  md = StartUpDG.make_periodic(md, periodicity)
  boundary_faces = StartUpDG.tag_boundary_faces(triangulateIO, dg.basis, md, boundary_dict)
  return DGMultiMesh(dg, GeometricTermsType(TriangulateIO(), dg), md, boundary_faces)
end

# TODO: DGMulti. Make `cells_per_dimension` a non-keyword argument for easier dispatch.
"""
    DGMultiMesh(dg::DGMulti; cells_per_dimension,
                coordinates_min=(-1.0, -1.0), coordinates_max=(1.0, 1.0),
                is_on_boundary=nothing,
                periodicity=ntuple(_ -> false, NDIMS))

Constructs a Cartesian [`DGMultiMesh`](@ref) with element type `dg.basis.elementType`. The domain is
the tensor product of the intervals `[coordinates_min[i], coordinates_max[i]]`.
- `is_on_boundary` specifies boundary using a `Dict{Symbol, <:Function}`
- `periodicity` is a tuple of `Bool`s specifying periodicity = `true`/`false` in the (x,y,z) direction.
"""
function DGMultiMesh(dg::DGMulti{NDIMS}; cells_per_dimension,
                     coordinates_min=ntuple(_ -> -one(real(dg)), NDIMS),
                     coordinates_max=ntuple(_ -> one(real(dg)), NDIMS),
                     is_on_boundary=nothing,
                     periodicity=ntuple(_ -> false, NDIMS), kwargs...) where {NDIMS}

  if haskey(kwargs, :is_periodic)
    # TODO: DGMulti. Deprecate `is_periodic` in version 0.5
    Base.depwarn("keyword argument `is_periodic` is now `periodicity`.", :DGMultiMesh)
    periodicity=kwargs[:is_periodic]
  end

  vertex_coordinates, EToV = StartUpDG.uniform_mesh(dg.basis.elementType, cells_per_dimension...)
  domain_lengths = coordinates_max .- coordinates_min
  for i in 1:NDIMS
    @. vertex_coordinates[i] = 0.5 * (vertex_coordinates[i] + 1) * domain_lengths[i] + coordinates_min[i]
  end

  md = MeshData(vertex_coordinates, EToV, dg.basis)
  if NDIMS == 1
    md = StartUpDG.make_periodic(md, periodicity...)
  else
    md = StartUpDG.make_periodic(md, periodicity)
  end
  boundary_faces = StartUpDG.tag_boundary_faces(md, is_on_boundary)
  return DGMultiMesh(dg, GeometricTermsType(Cartesian(), dg), md, boundary_faces)
end

"""
    DGMultiMesh(dg::DGMulti{NDIMS}, cells_per_dimension, mapping;
                is_on_boundary=nothing,
                periodicity=ntuple(_ -> false, NDIMS), kwargs...) where {NDIMS}

Constructs a `Curved()` [`DGMultiMesh`](@ref) with element type `dg.basis.elementType`.
- `mapping` is a function which maps from a reference [-1, 1]^NDIMS domain to a mapped domain,
   e.g., `xy = mapping(x, y)` in 2D.
- `is_on_boundary` specifies boundary using a `Dict{Symbol, <:Function}`
- `periodicity` is a tuple of `Bool`s specifying periodicity = `true`/`false` in the (x,y,z) direction.
"""
function DGMultiMesh(dg::DGMulti{NDIMS}, cells_per_dimension, mapping;
                     is_on_boundary=nothing,
                     periodicity=ntuple(_ -> false, NDIMS), kwargs...) where {NDIMS}

  vertex_coordinates, EToV = StartUpDG.uniform_mesh(dg.basis.elementType, cells_per_dimension...)
  md = MeshData(vertex_coordinates, EToV, dg.basis)
  md = NDIMS==1 ? StartUpDG.make_periodic(md, periodicity...) : StartUpDG.make_periodic(md, periodicity)

  @unpack xyz = md
  for i in eachindex(xyz[1])
    new_xyz = mapping(getindex.(xyz, i)...)
    setindex!.(xyz, new_xyz, i)
  end
  md_curved = MeshData(dg.basis, md, xyz...)

  # interpolate geometric terms to both volume and face cubature points
  @unpack rstxyzJ = md_curved
  @unpack Vq, Vf = dg.basis
  rstxyzJ_interpolated = map(x -> [Vq; Vf] * x, rstxyzJ)
  md_curved = @set md_curved.rstxyzJ = rstxyzJ_interpolated

  boundary_faces = StartUpDG.tag_boundary_faces(md_curved, is_on_boundary)
  return DGMultiMesh(dg, GeometricTermsType(Curved(), dg), md_curved, boundary_faces)
end

# Todo: DGMulti. Add traits for dispatch on affine/curved meshes here.

# Matrix type for lazy construction of physical differentiation matrices
# Constructs a lazy linear combination of B = ∑_i coeffs[i] * A[i]
struct LazyMatrixLinearCombo{Tcoeffs, N, Tv, TA <: AbstractMatrix{Tv}} <: AbstractMatrix{Tv}
  matrices::NTuple{N, TA}
  coeffs::NTuple{N, Tcoeffs}
  function LazyMatrixLinearCombo(matrices, coeffs)
    @assert all(matrix -> size(matrix) == size(first(matrices)), matrices)
    new{typeof(first(coeffs)), length(matrices), eltype(first(matrices)), typeof(first(matrices))}(matrices, coeffs)
  end
end
Base.eltype(A::LazyMatrixLinearCombo) = eltype(first(A.matrices))
Base.IndexStyle(A::LazyMatrixLinearCombo) = IndexCartesian()
Base.size(A::LazyMatrixLinearCombo) = size(first(A.matrices))

@inline function Base.getindex(A::LazyMatrixLinearCombo{<:Real, N}, i, j) where {N}
  val = zero(eltype(A))
  for k in Base.OneTo(N)
    val = val + A.coeffs[k] * getindex(A.matrices[k], i, j)
  end
  return val
end

# `SimpleKronecker` lazily stores a Kronecker product `kron(ntuple(A, NDIMS)...)`.
# This object also allocates some temporary storage to enable the fast computation
# of matrix-vector products.
struct SimpleKronecker{NDIMS, TA, Ttmp}
  A::TA
  tmp_storage::Ttmp # temporary array used for Kronecker multiplication
end

# constructor for SimpleKronecker which requires specifying only `NDIMS` and
# the 1D matrix `A`.
function SimpleKronecker(NDIMS, A, eltype_A=eltype(A))
  @assert size(A, 1) == size(A, 2) # check if square
  tmp_storage=[zeros(eltype_A, ntuple(_ -> size(A, 2), NDIMS)...) for _ in 1:Threads.nthreads()]
  return SimpleKronecker{NDIMS, typeof(A), typeof(tmp_storage)}(A, tmp_storage)
end

# Computes `b = kron(A, A) * x` in an optimized fashion
function LinearAlgebra.mul!(b_in, A_kronecker::SimpleKronecker{2}, x_in)

  @unpack A = A_kronecker
  tmp_storage = A_kronecker.tmp_storage[Threads.threadid()]
  n = size(A, 2)

  # copy `x_in` to `tmp_storage` to avoid mutating the input
  @assert length(tmp_storage) == length(x_in)
  for i in eachindex(tmp_storage)
    tmp_storage[i] = x_in[i]
  end
  x = reshape(tmp_storage, n, n)
  b = reshape(b_in, n, n)

  @turbo for j in 1:n, i in 1:n
    tmp = zero(eltype(x))
    for ii in 1:n
      tmp = tmp + A[i, ii] * x[ii, j]
    end
    b[i, j] = tmp
  end

  @turbo for j in 1:n, i in 1:n
    tmp = zero(eltype(x))
    for jj in 1:n
      tmp = tmp + A[j, jj] * b[i, jj]
    end
    x[i, j] = tmp
  end

  @turbo for i in eachindex(b_in)
    b_in[i] = x[i]
  end

  return nothing
end

# Computes `b = kron(A, A, A) * x` in an optimized fashion
function LinearAlgebra.mul!(b_in, A_kronecker::SimpleKronecker{3}, x_in)

  @unpack A = A_kronecker
  tmp_storage = A_kronecker.tmp_storage[Threads.threadid()]
  n = size(A, 2)

  # copy `x_in` to `tmp_storage` to avoid mutating the input
  for i in eachindex(tmp_storage)
    tmp_storage[i] = x_in[i]
  end
  x = reshape(tmp_storage, n, n, n)
  b = reshape(b_in, n, n, n)

  @turbo for k in 1:n, j in 1:n, i in 1:n
    tmp = zero(eltype(x))
    for ii in 1:n
      tmp = tmp + A[i, ii] * x[ii, j, k]
    end
    b[i, j, k] = tmp
  end

  @turbo for k in 1:n, j in 1:n, i in 1:n
    tmp = zero(eltype(x))
    for jj in 1:n
      tmp = tmp + A[j, jj] * b[i, jj, k]
    end
    x[i, j, k] = tmp
  end

  @turbo for k in 1:n, j in 1:n, i in 1:n
    tmp = zero(eltype(x))
    for kk in 1:n
      tmp = tmp + A[k, kk] * x[i, j, kk]
    end
    b[i, j, k] = tmp
  end

  return nothing
end

end # @muladd


text:
isrev(metnet::MetNet, ider::IDER_TYPE) = (indx = rxnindex(metnet, ider);
            metnet.lb[indx] < 0.0 && metnet.ub[indx] > 0.0)

isblock(metnet::MetNet, ider::IDER_TYPE) = (indx = rxnindex(metnet, ider);
    metnet.lb[indx] == 0.0 && metnet.ub[indx] == 0.0)

isopen(metnet::MetNet, ider::IDER_TYPE) = !isblock(metnet::MetNet, ider::IDER_TYPE)

isfwd_bounded(metnet::MetNet, ider::IDER_TYPE) = (indx = rxnindex(metnet, ider);
    metnet.lb[indx] >= 0.0 && metnet.ub[indx] > 0.0)

isbkwd_bounded(metnet::MetNet, ider::IDER_TYPE) = (indx = rxnindex(metnet, ider);
    metnet.lb[indx] < 0.0 && metnet.ub[indx] <= 0.0)

isfwd_defined(metnet::MetNet, ider::IDER_TYPE) = (indx = rxnindex(metnet, ider);
    length(rxn_reacts(metnet, indx)) > 0)

isbkwd_defined(metnet::MetNet, ider::IDER_TYPE) = (indx = rxnindex(metnet, ider);
    length(rxn_prods(metnet, indx)) > 0)

isfixxed(metnet::MetNet, ider::IDER_TYPE) = (indx = rxnindex(metnet, ider);
    metnet.lb[indx] == metnet.ub[indx] != 0.0)

revs(metnet::MetNet) = findall((metnet.lb .< 0.0) .& (metnet.ub .> 0.0))
revscount(metnet::MetNet) = length(revs(metnet))

blocks(metnet::MetNet) = findall((metnet.lb .== 0.0) .& (metnet.ub .== 0.0))
blockscount(metnet::MetNet) = length(blocks(metnet))

fwds_bounded(metnet::MetNet) = findall((metnet.lb .>= 0.0) .& (metnet.ub .> 0.0))
fwds_boundedcount(metnet::MetNet) = length(fwds_bounded(metnet))

bkwds_bounded(metnet::MetNet) = findall((metnet.lb .< 0.0) .& (metnet.ub .<= 0.0))
bkwds_boundedcount(metnet::MetNet) = length(bkwds_bounded(metnet))

fixxeds(metnet::MetNet) = findall((metnet.lb .== metnet.ub .!= 0.0))
fixxedscount(metnet::MetNet) = length(fixxeds(metnet))

allfwd(metnet::MetNet) = fwds_boundedcount(metnet) == rxnscount(metnet)

function is_exchange(metnet::MetNet, ider::IDER_TYPE)
    reacts = rxn_reacts(metnet, ider)
    prods = rxn_prods(metnet, ider)
    return xor(isempty(reacts), isempty(prods))
end
exchanges(metnet::MetNet) = findall(x -> is_exchange(metnet, x), metnet.rxns)
exchangescount(metnet::MetNet) = length(exchanges(metnet))


text:
using OCToolbox
using QuantumOpticsBase
using Test

b = SpinBasis(1//2)
sx = sigmax(b)
sy = sigmay(b)
sz = sigmaz(b)
example_control = [0.810635963493793
, 0.46484987867321426
, 0.34970273151618936
, 0.7916035280985718
, 0.3613862392184721
, 0.619927796151712
, 0.6967410801849956
, 0.055413818931025594
, 0.5452231997511754
, 0.3901904774054561]

@testset "Test_utils.jl" begin
    @test commutator(sx, sx) == 0.0*sx
    @test commutator(sx, sy) == 2.0im*sz
    @test count(x->x==0.5, clip_controls(example_control, 0.5)) == 5

end


text:
# Test the solvehook's options
@test_logs (:warn, "Problem requires use of active set constraint. Setting active_constr = true") solve(m) == :Optimal
@test solve(m, active_constr = true) == :Optimal && abs(getflexibilityindex(m) - 3.6) <= 1e-2
@test solve(m, active_constr = true, diag = true) == :Optimal && abs(getflexibilityindex(m) - 3.6) <= 1e-2
@test solve(m, active_constr = true, conic_δ = true) == :Optimal && abs(getflexibilityindex(m) - 3.6) <= 1e-2
@test solve(m, active_constr = true) == :Optimal && abs(getflexibilityindex(m) - 3.6) <= 1e-2
@test_logs (:warn, "The slack upper bound (U) value of 1000 is too low for this problem. Thus, U is being set to 10000.0 !") solve(m, active_constr = true, U = 1000) == :Optimal

# Test with other uncertainty sets
dev = ones(flex_data.numRVs) * 10
setuncertaintyset(m, :Hyperbox, [[dev]; [dev]])
@test solve(m, active_constr = true) == :Optimal && abs(getflexibilityindex(m) - 0.5) <= 1e-2
setuncertaintyset(m, :PNorm, 1)
@test solve(m, active_constr = true) == :Optimal && abs(getflexibilityindex(m) - 6.667) <= 1e-2
setuncertaintyset(m, :PNorm, 2)
@test solve(m, active_constr = true) == :Optimal && abs(getflexibilityindex(m) - 6.325) <= 1e-2
setuncertaintyset(m, :PNorm, Inf)
@test solve(m, active_constr = true) == :Optimal && abs(getflexibilityindex(m) - 5) <= 1e-2

# Check safeguards
dev = ones(flex_data.numRVs - 1) * 10
setuncertaintyset(m, :Hyperbox, [[dev]; [dev]])
@test_throws ErrorException solve(m, active_constr = true) == :Optimal
setuncertaintyset(m, :Ellipsoid, Matrix(I, 2, 2))
@test_throws ErrorException solve(m, active_constr = true) == :Optimal
setuncertaintyset(m, :Ellipsoid, covar)
@test solve(m, active_constr = true) == :Optimal && abs(getflexibilityindex(m) - 3.6) <= 1e-2

# Test the solution time
@test getsolutiontime(m) != nothing