text:
module MLBasePlotting

export plotperf, area_under_pr, area_under_roc

using MLBase
using Gadfly

function plotperf(sortedTruth::AbstractVector; curveType="pr", name="", resolution=600)
    methods = Dict()
    methods["predictor"] = (sortedTruth, -collect(1:length(sortedTruth)))
    plotperf(methods, curveType=curveType, name=name, resolution=resolution)
end
function plotperf(truth::AbstractVector, predictor::AbstractVector; curveType="pr", name="", resolution=600)
    methods = Dict()
    methods["predictor"] = (truth, predictor)
    plotperf(methods, curveType=curveType, name=name, resolution=resolution)
end
function plotperf(methods; curveType="pr", name="", resolution=600)

    if curveType == "roc"
        warn("The ROC curve type is not yet fully debugged.")
        xlabel = "False Positive Rate"
        ylabel = "True Positive Rate"
        xmap = false_positive_rate
        ymap = true_positive_rate
        if name == "" name = "ROC" end
    elseif curveType == "pr"
        xlabel = "Recall"
        ylabel = "Precision"
        xmap = recall
        ymap = precision
        if name == "" name = "Precision/Recall" end
    end

    # generate the random predictors
    xmean = zeros(resolution)
    ymean = zeros(resolution)
    layers = Any[]
    truth = methods[first(keys(methods))][1]
    numRandom = 10
    for i in 1:numRandom
        predictor = invperm(sortperm(rand(length(truth))))
        rocData = MLBase.roc(round(Int64, truth), float(predictor), resolution)
        vals = collect(map(x->(xmap(x), -ymap(x)), rocData))
        sort!(vals)
        xvals = map(x->x[1], vals)
        yvals = map(x->-x[2], vals)
        xmean .+= xvals
        ymean .+= yvals
        aucValue = @sprintf("%0.03f", MLBasePlotting.area_under_curve(xvals, yvals))
        push!(layers, layer(
            x=[0.0; xvals], y=[yvals[1]; yvals],
            Geom.line,
            Theme(default_color=colorant"lightgrey")
        ))
    end
    xmean /= numRandom
    ymean /= numRandom
    push!(layers, layer(
        x=xmean, y=ymean,
        Geom.line,
        Theme(default_color=colorant"grey")
    ))

    # plot the mean of the random predictors
    labels = ASCIIString[]
    xdata = Float64[]
    ydata = Float64[]
    aucValue = @sprintf("%0.03f", area_under_curve([0.0; xmean], [ymean[1]; ymean]))
    append!(labels, collect(repeat(["AUC = $aucValue, random"], inner=[length(xmean)+1])))
    append!(xdata, [0.0; xmean])
    append!(ydata, [ymean[1]; ymean])

    # plot the actual passed data
    for (key,(truth,predictor)) in methods
        rocData = MLBase.roc(round(Int64, truth), float(invperm(sortperm(predictor))), resolution)
        vals = collect(map(x->(xmap(x), -ymap(x)), rocData))
        sort!(vals)
        xvals = map(x->x[1], vals)
        yvals = map(x->-x[2], vals)
        aucValue = @sprintf("%0.03f", MLBasePlotting.area_under_curve([0.0; xvals], [yvals[1]; yvals]))
        append!(labels, collect(repeat(["AUC = $aucValue, $key"], inner=[length(xvals)+1])))
        append!(xdata, [0.0; xvals])
        append!(ydata, [yvals[1]; yvals])
    end

    plot(
        Guide.title("$name"),
        Guide.XLabel(xlabel),
        Guide.YLabel(ylabel),
        Guide.colorkey("Methods"),
        layer(
            x=xdata, y=ydata, color=labels,
            Geom.line
        ),
        layers...,
        Scale.color_discrete_manual(["grey", "blue", "red", "green", "purple", "pink", "orange"][1:length(methods)+1]...)
    )
end

function area_under_pr(truth::AbstractVector, predictor::AbstractVector; resolution=4000)
    rocData = MLBase.roc(round(Int64, truth), float(invperm(sortperm(predictor))), resolution)
    vals = collect(map(x->(recall(x), -precision(x)), rocData))
    sort!(vals)
    xvals = map(x->x[1], vals)
    yvals = map(x->-x[2], vals)
    area_under_curve([0.0; xvals], [yvals[1]; yvals]) # make sure we extend all the way to zero
end
function area_under_roc(truth::AbstractVector, predictor::AbstractVector; resolution=4000)
    rocData = MLBase.roc(round(Int64, truth), float(invperm(sortperm(predictor))), resolution)
    vals = collect(map(x->(false_positive_rate(x), -true_positive_rate(x)), rocData))
    sort!(vals)
    xvals = map(x->x[1], vals)
    yvals = map(x->-x[2], vals)
    area_under_curve([0.0; xvals], [yvals[1]; yvals]) # make sure we extend all the way to zero
end

# handles NaN values by using previous valid values
function area_under_curve(x, y) # must be sorted by increasing x
    area = 0.0
    lastVal = NaN
    for i in 2:length(x)
        v = (y[i-1]+y[i])/2 * (x[i]-x[i-1])
        if !isnan(v)
            area += v
            lastVal = v
        elseif !isnan(lastVal)
            area += lastVal
        end
    end
    area
end

end # module


text:
using CurveProximityQueries
using ConvexProximityQueries
using StaticArrays
using Plots
using Random: seed!

function movecp(t; k=3)
   x1 = 3*cos(k*t)*cos(t)
   y1 = cos(k*t)*sin(t)
   x2 = 3*(sin(t)+1.0)
   y2 = 3*cos(t)
   SVector{2}(x1, y1), SVector{2}(x2, y2)
end

anim = @animate for i=1:360
    t = i/180*π
    b = [SVector{2}(0.0, 0.0),
         movecp(t)...,
         SVector{2}(3.0, 0.0)]
    B = Bernstein(b...)

    seed!(1)
    A = randpoly([1.5, 0.0], t; scale=1.0)

    plot(A)
    plot!(xlims = (-1.5,4.5), ylims = (-2.0,2.0))

    @show t
    if collision_detection(A, B)
        plot!(B, linecolor = RGB(213/255,99/255,92/255), linewidth = 4)
    else
        AB = closest_points(A, B)
        plot!(AB, linewidth = 3)
        plot!(B, linecolor = RGB(106/255,92/225,213/255), linewidth = 4)
    end
end

gif(anim, "logo1.gif", fps = 30)

anim = @animate for i=1:360
    t = i/180*π
    b = [SVector{2}(0.0, 0.0),
         movecp(t)...,
         SVector{2}(3.0, 0.0)]
    B = Bernstein(b...)

    a = [SVector{2}(0.0, 1.5),
         movecp(t; k=2)...,
         SVector{2}(3.0, 1.5)]
    A = Bernstein(a...)

    plot(xlims = (-1.5,4.5), ylims = (-2.0,2.0))

    @show t
    if collision_detection(A, B; atol=1e-7)
        plot!(B, linecolor = RGB(213/255,99/255,92/255), linewidth = 4)
        plot!(A, linecolor = RGB(213/255,99/255,92/255), linewidth = 4)
    else
        AB = closest_points(A, B; atol=1e-7)
        plot!(AB, linewidth = 3)
        plot!(B, linecolor = RGB(106/255,92/225,213/255), linewidth = 4)
        plot!(A, linecolor = RGB(120/255,213/255,104/255), linewidth = 4)
    end
end

gif(anim, "logo2.gif", fps = 30)

ext: jl


text:
# this function is used to get the
# suffix of the glUniform function
# for the correct type
_type_suffix(::Type{<:Bool}) = "i"
_type_suffix(::Type{<:Signed}) = "i"
_type_suffix(::Type{<:Unsigned}) = "ui"
_type_suffix(::Type{<:AbstractFloat}) = "f"

"""
        uniform!(location, values...) -> Nothing

Set the uniform `vec` variable at `location`.

# Examples

```julia
myuniform = getuniform(program, "someUniform")
# set float vec3
uniform!(myuniform, 1.0, 1.0, 0.5)
# set float vec4
uniform!(myuniform, Cfloat[1, 2, 3, 4])
# set a 4x4 matrix
uniform!(myuniform, rand(Cfloat, 4, 4))
```
"""
function uniform! end

@generated function uniform!(location::GLint, values::Vararg{T, N}) where {N, T <: Real}
    suffix = _type_suffix(T)
    glFunc = Symbol("glUniform$(N)$suffix")
    return :( $(glFunc)(location, values...) )
end

@generated function uniform!(location::GLint, vector::SVector{N, T}) where {N, T <: Real}
    suffix = _type_suffix(T)
    glFunc = Symbol("glUniform$(N)$(suffix)v")
    return :( $(glFunc)(location, 1, vector) )
end

@generated function uniform!(location::GLint, matrix::SMatrix{N, M, T}, transposed::Bool = false) where {N, M, T <: Real}
    glFunc = Symbol(
        if N == M
            "glUniformMatrix$(N)fv"
        else
            "glUniformMatrix$(N)x$(M)fv"
        end
    )
    return :( $(glFunc)(location, 1, transposed, Cfloat[matrix...]) )
end

function uniform!(location::GLint, vector::AbstractVector{T}) where {T <: Real}
    N = length(vector)
    uniform!(location, SVector{N, T}(vector...))
end

function uniform!(location::GLint, matrix::AbstractMatrix{T}, transposed::Bool = false) where {T <: Real}
    N, M = size(matrix)
    uniform!(location, SMatrix{N, M, T}(matrix), transposed)
end
"""
        getuniform(program, name) -> GLint

Gets the location of the uniform variable
identified by `name`.
"""
getuniform(program::GLuint, name::SymString) = glGetUniformLocation(program, string(name))
"""
        @uniforms foo, bar, ... = program

Get the location of `foo`, `bar` and `etc` from `program`,
and assign each uniform to the corresponding variable.

This macro transform the following expression:

```julia
@uniforms foo, bar, qux: baz = program
```

Into:
```julia
foo = getuniform(program, :foo)
bar = getuniform(program, :bar)
qux = getuniform(program, :baz)
```

You can also rename the variables as following:

```julia
@uniforms a: foo, b: bar = program
```

The variable `a` contains the location of `foo`,
and `b` contains the loctation of `bar`.
"""
macro uniforms(ex)
    @assert ex.head == :(=)
    program = esc(ex.args[2])
    if !Meta.isexpr(ex.args[1], :tuple)
        args = [ex.args[1]]
    else
        args = ex.args[1].args
    end
    return Expr(:block,
        map(args) do arg
            if arg isa Symbol
                alias, name = arg, arg
            elseif Meta.isexpr(arg, :(=))
                alias, name = arg.args
            elseif Meta.isexpr(arg, :call)
                alias, name = arg.args[2:3]
            end
            return Expr(:(=), esc(alias), Expr(:call, :getuniform, program, QuoteNode(name)))
        end...
    )
end

macro uniforms!(ex)
    args = ex.args
    for i in eachindex(args)
        if Meta.isexpr(args[i], :(:=))
            name, value = args[i].args
            args[i] = Expr(:call, :uniform!, name, value)
        end
    end
    return esc(Expr(:block, args...))
end


text:
# Line geometry connects (x, y) coordinates with lines.

struct LineGeometry <: Gadfly.GeometryElement
    default_statistic::Gadfly.StatisticElement

    # Do not reorder points along the x-axis.
    preserve_order::Bool

    order::Int

    tag::Symbol
end

function LineGeometry(default_statistic=Gadfly.Stat.identity();
                      preserve_order=false, order=2, tag=empty_tag)
    LineGeometry(default_statistic, preserve_order, order, tag)
end

const line = LineGeometry

function contour(; levels=15, samples=150, preserve_order=true)
    return LineGeometry(Gadfly.Stat.contour(levels=levels, samples=samples),
                                            preserve_order=preserve_order)
end

# Only allowing identity statistic in paths b/c I don't think any
# any of the others will work with `preserve_order=true` right now
path() = LineGeometry(preserve_order=true)

density(; bandwidth::Real=-Inf) =
    LineGeometry(Gadfly.Stat.density(bandwidth=bandwidth))

density2d(; bandwidth::Tuple{Real,Real}=(-Inf,-Inf), levels=15) =
    LineGeometry(Gadfly.Stat.density2d(bandwidth=bandwidth, levels=levels); preserve_order=true)

smooth(; method::Symbol=:loess, smoothing::Float64=0.75) =
    LineGeometry(Gadfly.Stat.smooth(method=method, smoothing=smoothing), order=5)

step(; direction::Symbol=:hv) = LineGeometry(Gadfly.Stat.step(direction=direction))

default_statistic(geom::LineGeometry) = geom.default_statistic

element_aesthetics(::LineGeometry) = [:x, :y, :color, :group]

# Render line geometry.
#
# Args:
#   geom: line geometry.
#   theme: the plot's theme.
#   aes: aesthetics.
#
# Returns:
#   A compose Form.
#
function render(geom::LineGeometry, theme::Gadfly.Theme, aes::Gadfly.Aesthetics)
    Gadfly.assert_aesthetics_defined("Geom.line", aes, :x, :y)
    Gadfly.assert_aesthetics_equal_length("Geom.line", aes,
                                          element_aesthetics(geom)...)

    default_aes = Gadfly.Aesthetics()
    default_aes.color = discretize_make_ia(RGBA{Float32}[theme.default_color])
    aes = inherit(aes, default_aes)

    ctx = context(order=geom.order)
    XT, YT, CT = eltype(aes.x), eltype(aes.y), eltype(aes.color)
    XYT = Tuple{XT, YT}

    line_style = Gadfly.get_stroke_vector(theme.line_style)

    if aes.group != nothing
        GT = eltype(aes.group)

        if !geom.preserve_order
            p = sortperm(aes.x)
            aes_group = aes.group[p]
            aes_color = length(aes.color) > 1 ? aes.color[p] : aes.color
            aes_x = aes.x[p]
            aes_y = aes.y[p]
        else
            aes_group = copy(aes.group)
            aes_color = copy(aes.color)
            aes_x = copy(aes.x)
            aes_y = copy(aes.y)
        end

        # organize x, y pairs into lines
        if length(aes_group) > length(aes_color)
            p = sortperm(aes_group)
        elseif length(aes_color) > length(aes_group)
            p = sortperm(aes_color, lt=Gadfly.color_isless)
        else
            p = sortperm(collect((Tuple{GT, CT}),zip(aes_group, aes_color)),
                         lt=Gadfly.group_color_isless)
        end
        permute!(aes_group, p)
        permute!(aes_color, p)
        permute!(aes_x, p)
        permute!(aes_y, p)

        points = Vector{XYT}[]
        points_colors = CT[]
        points_groups = GT[]

        first_point = true
        for (i, (x, y, c, g)) in enumerate(zip(aes_x, aes_y, cycle(aes_color), aes_group))
            if !isconcrete(x) || !isconcrete(y)
                first_point = true
                continue
            end

            if i > 1 && (c != points_colors[end] || g != points_groups[end])
                first_point = true
            end

            if first_point
                push!(points, XYT[])
                push!(points_colors, c)
                push!(points_groups, g)
                first_point = false
            end

            push!(points[end], (x, y))
        end

        classes = [svg_color_class_from_label(aes.color_label([c])[1])
                   for (c, g) in zip(points_colors, points_groups)]

        ctx = compose!(ctx, (context(), Compose.line(points,geom.tag),
                        stroke(points_colors),
                        strokedash(line_style),
                        svgclass(classes)),
                      svgclass("geometry"))

    elseif length(aes.color) == 1 &&
            !(isa(aes.color, IndirectArray) && count(!ismissing, aes.color.values) > 1)
        T = (Tuple{eltype(aes.x), eltype(aes.y)})
        points = T[(x, y) for (x, y) in zip(aes.x, aes.y)]
        geom.preserve_order || sort!(points, by=first)
        ctx = compose!(ctx, (context(), Compose.line([points],geom.tag),
                       stroke(aes.color[1]),
                       strokedash(line_style)),
                       svgclass("geometry"))
    else
        if !geom.preserve_order
            p = sortperm(aes.x)
            aes_color = aes.color[p]
            aes_x = aes.x[p]
            aes_y = aes.y[p]
        else
            aes_color = copy(aes.color)
            aes_x = copy(aes.x)
            aes_y = copy(aes.y)
        end

        # organize x, y pairs into lines
        p = sortperm(aes_color, lt=Gadfly.color_isless)
        permute!(aes_color, p)
        permute!(aes_x, p)
        permute!(aes_y, p)

        points = Vector{XYT}[]
        points_colors = CT[]

        first_point = true
        for (i, (x, y, c)) in enumerate(zip(aes_x, aes_y, aes_color))
            if !isconcrete(x) || !isconcrete(y)
                first_point = true
                continue
            end

            if isempty(points_colors) || c != points_colors[end]
                first_point = true
            end

            if first_point
                push!(points, XYT[])
                push!(points_colors, c)
                first_point = false
            end

            push!(points[end], (x, y))
        end

        classes = [svg_color_class_from_label(aes.color_label([c])[1])
                   for c in points_colors]

        ctx = compose!(ctx, (context(), Compose.line(points,geom.tag),
                        stroke(points_colors),
                        strokedash(line_style),
                        svgclass(classes)),
                      svgclass("geometry"))
    end

    return compose!(ctx, fill(nothing), linewidth(theme.line_width))
end


text:
struct NLEProblem{T}
    F::T
end
function value(nleq::NLEProblem, x)
    nleq.F(x)
end

struct NLEOptions{T, Tmi}
    f_abstol::T
    f_reltol::T
    maxiter::Tmi
end
NLEOptions(; f_abstol=1e-8, f_reltol=1e-8, maxiter=10^4) = NLEOptions(f_abstol, f_reltol, maxiter)

include("secant.jl")
include("newton.jl")
#include("inplace_loop.jl")
#include("outofplace_loop.jl")


text:
#export IndexGraph

#using LightGraphs
"""
The data structure that stores a graph and a mapping from indices to their
representations in the graph
"""
mutable struct IndexGraph{T}
  idx::Dict{T,Int64}
  graph::Graph
end

function IndexGraph{T}(nodes::AbstractArray{T}, edges::AbstractArray{Tuple{T,T}})
  idx = Dict{T, Int64}(n => i for (i,n) in enumerate(nodes))
  graph = Graph(length(idx))
  for e in edges
    if !haskey(idx, e[1]) || !haskey(idx,e[2])
      error("Edges may only contain nodes in the graph")
    else
      add_edge!(graph, idx[e[1]], idx[e[2]])
    end
  end
  IndexGraph(idx, graph)
end

function IndexGraph{T}(nodes::AbstractArray{T})
  idx = Dict{T, Int64}(n => i for (i,n) in enumerate(nodes))
  graph = Graph(length(idx))
  IndexGraph(idx, graph)
end

function IndexGraph(g::Graph)
  idx = Dict{Int64, Int64}(i => i for i in 1:nv(g))
  graph = g
  IndexGraph(idx, graph)
end

"""
Helper function to embed an index graph into multidimensional losses.
"""
function embed_graph(ig::IndexGraph, yidxs::Array)
  nodes = 1:yidxs[end][end]
  newedges = Tuple{Int64,Int64}[]
  for edge in edges(ig.graph)
    push!(newedges, (yidxs[edge.src], yidxs[edge.dst]))
  end
  IndexGraph(nodes, newedges)
end


text:
import LinearAlgebra: cholesky, diag, Diagonal, Hermitian
import Printf: @printf
import Random
import Statistics: mean, median, std
using Test
import RandomCorrelationMatrices

# This is not a numerical test. The code just generates some examples.
# However, if it fails to run, thats a good sign something is wrong.
Random.seed!(12345)
println("Random correlation matrix")
n = 4
for η in [2, 8, 32, 128]
    r = RandomCorrelationMatrices.randcormatrix(n, η)
    @show η
    for i in 1:n
        for j in 1:n
            @printf("%5.2f ", r[i,j])
        end
        println()
    end

    ranges = Float64[]
    for s in 1:1000
        r = RandomCorrelationMatrices.randcormatrix(n, η)
        r0 = r - Matrix(Diagonal((diag(r))))
        push!(ranges, maximum(vec(r0)) - minimum(vec(r0)))
    end
    @show mean(ranges)
    @show median(ranges)
end

@testset "Random covariance matrix" begin
    Random.seed!(12345)
    σ = [100, 200, 300, 400, 500]
    Σ = RandomCorrelationMatrices.randcovmatrix(5, 10, σ)
    L = cholesky(Hermitian(Σ)).L
    samples = 100_000
    results = zeros(5, samples)
    for i in 1:samples
        results[:,i] = L*randn(5)
    end
    means = mean(results, dims=2)
    @test all(abs.(means) .<= 5)
    stddevs = std(results, dims=2)
    @test all(abs.(stddevs - σ) .<= 5)
    @show σ
    @show round.(stddevs, digits=1)
end


text:
tex = Input("f(x) = \\int_{-\\infty}^\\infty
           \\hat f(\\xi)\\,e^{2 \\pi i \\xi x}
           \\,d\\xi")
mode = Input(false)

function main(window)
    push!(window.assets, "latex")
    push!(window.assets, "widgets")

    lift(tex, mode) do t, m
        vbox(h1("LaTeX"),
             hbox("LaTeX support is via ", hskip(1em), latex("\\KaTeX")),
             textinput(t) >>> tex,
             hbox("Show as a block", hskip(1em), checkbox(value=m) >>> mode),
             vskip(1em),
             latex(t, block=m)) |> pad(1em) |> maxwidth(800px)
    end
end


text:
module Fields

export Field, Initialize_Field

struct Field
   name::String

   spin::Int64
   boost::Int64
   falloff::Int64

   mv::Int64

   nx::Int64
   ny::Int64

   n  ::Array{ComplexF64,2}
   tmp::Array{ComplexF64,2}
   np1::Array{ComplexF64,2}

   k::Array{ComplexF64,2}

   rad_d1::Array{ComplexF64,2}
   rad_d2::Array{ComplexF64,2}

   raised ::Array{ComplexF64,2}
   lowered::Array{ComplexF64,2}
   sph_lap::Array{ComplexF64,2}

   edth       ::Array{ComplexF64,2}
   edth_prime ::Array{ComplexF64,2}
   thorn      ::Array{ComplexF64,2}
   thorn_prime::Array{ComplexF64,2}

   Field(;
         name,
         spin, boost, falloff, mv,
         nx, ny
        ) = new(name,
                spin, boost, falloff, mv,
                nx, ny,
                zeros(ComplexF64,nx,ny),
                zeros(ComplexF64,nx,ny),
                zeros(ComplexF64,nx,ny),
                zeros(ComplexF64,nx,ny),
                zeros(ComplexF64,nx,ny),
                zeros(ComplexF64,nx,ny),
                zeros(ComplexF64,nx,ny),
                zeros(ComplexF64,nx,ny),
                zeros(ComplexF64,nx,ny),
                zeros(ComplexF64,nx,ny),
                zeros(ComplexF64,nx,ny),
                zeros(ComplexF64,nx,ny),
                zeros(ComplexF64,nx,ny)
               )
end

function Initialize_Field(;
      name::String,
      spin::Int64,
      boost::Int64,
      falloff::Int64,
      Mvals::Vector{Int64},
      nx::Int64,
      ny::Int64)
   return Dict([
      (mv,Field(name=name,spin=spin,boost=boost,falloff=falloff,mv=mv,nx=nx,ny=ny))
      for mv in Mvals
     ])
end

end


text:
@testset "solve" begin
    @testset "basic" begin
        @system SSolve(Controller) begin
            a ~ preserve(parameter)
            b ~ preserve(parameter)
            c ~ preserve(parameter)
            x(a, b, c) => begin
                a*x^2 + b*x + c
            end ~ solve
        end
        s1 = instance(SSolve, config=:SSolve => (a=1, b=2, c=1))
        @test s1.x' == -1
        s2 = instance(SSolve, config=:SSolve => (a=1, b=-2, c=1))
        @test s2.x' == 1
        s3 = instance(SSolve, config=:SSolve => (a=1, b=2, c=-3))
        @test s3.x' ∈ [1, -3]
        s4 = instance(SSolve, config=:SSolve => (a=1, b=-2, c=-3))
        @test s4.x' ∈ [-1, 3]
    end

    @testset "unit" begin
        @system SSolveUnit(Controller) begin
            a ~ preserve(u"m^-1", parameter)
            b ~ preserve(parameter)
            c ~ preserve(u"m", parameter)
            x(a, b, c) => begin
                a*x^2 + b*x + c
            end ~ solve(u"m")
        end
        s = instance(SSolveUnit, config=:SSolveUnit => (a=1, b=2, c=1))
        @test s.x' == -1u"m"
    end

    @testset "unit with scale" begin
        @system SSolveUnitScale(Controller) begin
            a ~ preserve(u"m^-1", parameter)
            b ~ preserve(u"cm/m", parameter) #HACK: caution with dimensionless unit!
            c ~ preserve(u"cm", parameter)
            x(a, b, c) => begin
                a*x^2 + b*x + c
            end ~ solve(u"cm")
        end
        s = instance(SSolveUnitScale, config=:SSolveUnitScale => (a=0.01, b=2, c=1))
        @test s.x' == -100u"cm"
    end

    @testset "linear" begin
        @system SSolveLinear(Controller) begin
            x => (2x ⩵ 1) ~ solve
        end
        s = instance(SSolveLinear)
        @test s.x' == 0.5
    end

    @testset "quadratic single" begin
        @system SSolveQuadraticSingle(Controller) begin
            x => (x^2 ⩵ 1) ~ solve
        end
        s = instance(SSolveQuadraticSingle)
        @test s.x' == 1
    end

    @testset "quadratic double" begin
        @system SSolveQuadraticDouble(Controller) begin
            l ~ preserve(parameter)
            u ~ preserve(parameter)
            x => (x^2 + 2x - 3) ~ solve(lower=l, upper=u)
        end
        s1 = instance(SSolveQuadraticDouble, config=:SSolveQuadraticDouble => (l=0, u=Inf))
        @test s1.x' == 1
        s2 = instance(SSolveQuadraticDouble, config=:SSolveQuadraticDouble => (l=-Inf, u=0))
        @test s2.x' == -3
    end

    @testset "cubic" begin
        @system SSolveCubic(Controller) begin
            l ~ preserve(parameter)
            u ~ preserve(parameter)
            x => ((x-5)*(x-15)*(x-25)) ~ solve(lower=l, upper=u)
        end
        s0 = instance(SSolveCubic, config=:SSolveCubic => (l=-Inf, u=0))
        @test s0.x' == 0
        s1 = instance(SSolveCubic, config=:SSolveCubic => (l=-Inf, u=10))
        @test s1.x' ≈ 5
        s2 = instance(SSolveCubic, config=:SSolveCubic => (l=10, u=20))
        @test s2.x' ≈ 15
        s3 = instance(SSolveCubic, config=:SSolveCubic => (l=20, u=30))
        @test s3.x' ≈ 25
        s4 = instance(SSolveCubic, config=:SSolveCubic => (l=30, u=Inf))
        @test s4.x' == 30
    end

    @testset "equal" begin
        @system SSolveEqual(Controller) begin
            x1 => (x1^2 + 2x1 - 3) ~ solve
            x2 => (x2^2 + 2x2 - 3 ⩵ 0) ~ solve
        end
        s = instance(SSolveEqual)
        @test s.x1' == s.x2'
    end
end


text:
using LinearAlgebra

"""
Functions to ensure gate call counts match expected number
"""
function ncu_gate_calls_default(n::Int)
    @assert n >= 2

    if n == 2
        return 5
    else
        return 3*ncu_gate_calls_default(n-1) + 2
    end
end

function ncu_gate_calls_opt(n::Int)
    @assert n >= 3
    if n == 3
        return 13
    else
        return 3*ncu_gate_calls_opt(n-1) + 2
    end
end

@testset "Test Grover algorithm" begin
    @testset "Test diffusion operator" begin
        cct = QXZoo.Circuit.Circ(5)
    end

    @testset "Test bitstring oracle operator" begin
        cct1 = QXZoo.Circuit.Circ(9)
        cct2 = QXZoo.Circuit.Circ(9)
        QXZoo.Grover.bitstring_ncu!(cct1, 11, collect(1:5), 9, DefaultGates.GateSymbols.c_x, collect(6:8))

        cct2 << DefaultGates.x(3)
        cct2 << DefaultGates.x(5)
        cct2 << DefaultGates.x(9)
        QXZoo.NCU.apply_ncu!(cct2, [5,8], [], 9, DefaultGates.GateSymbols.c_x)
        QXZoo.NCU.apply_ncu!(cct2, [4,7], [], 8, DefaultGates.GateSymbols.c_x)
        QXZoo.NCU.apply_ncu!(cct2, [3,6], [], 7, DefaultGates.GateSymbols.c_x)
        QXZoo.NCU.apply_ncu!(cct2, [1,2], [], 6, DefaultGates.GateSymbols.c_x)
        QXZoo.NCU.apply_ncu!(cct2, [3,6], [], 7, DefaultGates.GateSymbols.c_x)
        QXZoo.NCU.apply_ncu!(cct2, [4,7], [], 8, DefaultGates.GateSymbols.c_x)
        QXZoo.NCU.apply_ncu!(cct2, [5,8], [], 9, DefaultGates.GateSymbols.c_x)
        QXZoo.NCU.apply_ncu!(cct2, [4,7], [], 8, DefaultGates.GateSymbols.c_x)
        QXZoo.NCU.apply_ncu!(cct2, [3,6], [], 7, DefaultGates.GateSymbols.c_x)
        QXZoo.NCU.apply_ncu!(cct2, [1,2], [], 6, DefaultGates.GateSymbols.c_x)
        QXZoo.NCU.apply_ncu!(cct2, [3,6], [], 7, DefaultGates.GateSymbols.c_x)
        QXZoo.NCU.apply_ncu!(cct2, [4,7], [], 8, DefaultGates.GateSymbols.c_x)
        cct2 << DefaultGates.x(9)
        cct2 << DefaultGates.x(5)
        cct2 << DefaultGates.x(3)

        # X @ (0b00100, 0b10000)*2
        @test cct1.circ_ops.len == cct2.circ_ops.len

        @test cct1.gate_set == cct2.gate_set

        h1 = cct1.circ_ops.head
        h2 = cct2.circ_ops.head
        for i in 1:cct1.circ_ops.len
            @test h1.data === h2.data
            h1 = h1.next
            h2 = h2.next
        end
    end

    @testset "Test phase oracle operator" begin
        cct1 = QXZoo.Circuit.Circ(6)
        cct2 = QXZoo.Circuit.Circ(9)
        cct3 = QXZoo.Circuit.Circ(9)
        QXZoo.Grover.bitstring_phase_oracle!(cct1, 11, collect(1:5), 6)
        QXZoo.Grover.bitstring_phase_oracle!(cct2, 11, collect(1:5), 9, collect(6:8))

        cct3 << DefaultGates.x(3)
        cct3 << DefaultGates.x(5)
        cct3 << DefaultGates.x(9)
        cct3 << DefaultGates.h(9)
        QXZoo.NCU.apply_ncu!(cct3, [5,8], [], 9, DefaultGates.GateSymbols.c_x)
        QXZoo.NCU.apply_ncu!(cct3, [4,7], [], 8, DefaultGates.GateSymbols.c_x)
        QXZoo.NCU.apply_ncu!(cct3, [3,6], [], 7, DefaultGates.GateSymbols.c_x)
        QXZoo.NCU.apply_ncu!(cct3, [1,2], [], 6, DefaultGates.GateSymbols.c_x)
        QXZoo.NCU.apply_ncu!(cct3, [3,6], [], 7, DefaultGates.GateSymbols.c_x)
        QXZoo.NCU.apply_ncu!(cct3, [4,7], [], 8, DefaultGates.GateSymbols.c_x)
        QXZoo.NCU.apply_ncu!(cct3, [5,8], [], 9, DefaultGates.GateSymbols.c_x)
        QXZoo.NCU.apply_ncu!(cct3, [4,7], [], 8, DefaultGates.GateSymbols.c_x)
        QXZoo.NCU.apply_ncu!(cct3, [3,6], [], 7, DefaultGates.GateSymbols.c_x)
        QXZoo.NCU.apply_ncu!(cct3, [1,2], [], 6, DefaultGates.GateSymbols.c_x)
        QXZoo.NCU.apply_ncu!(cct3, [3,6], [], 7, DefaultGates.GateSymbols.c_x)
        QXZoo.NCU.apply_ncu!(cct3, [4,7], [], 8, DefaultGates.GateSymbols.c_x)
        cct3 << DefaultGates.h(9)
        cct3 << DefaultGates.x(9)
        cct3 << DefaultGates.x(5)
        cct3 << DefaultGates.x(3)

        # X @ (0b00100, 0b10000)*2 + 2 Hadamard
        @test cct1.circ_ops.len == ncu_gate_calls_opt(5) + 6
        @test cct2.gate_set == cct3.gate_set

        h2 = cct2.circ_ops.head
        h3 = cct3.circ_ops.head
        for i in 1:cct2.circ_ops.len
            @test h2.data === h3.data
            h2 = h2.next
            h3 = h3.next
        end
    end

    @testset "Test diffusion operator" begin
        cct1 = QXZoo.Circuit.Circ(9)
        cct2 = QXZoo.Circuit.Circ(9)
        QXZoo.Grover.apply_diffusion!(cct1, collect(1:5), 9, collect(6:8))

        for i in 1:5
            cct2 << DefaultGates.h(i)
            cct2 << DefaultGates.x(i)
        end
        cct2 << DefaultGates.h(9)
        cct2 << DefaultGates.x(9)

        cct2 << DefaultGates.h(9)
        QXZoo.NCU.apply_ncu!(cct2, [5,8], [], 9, DefaultGates.GateSymbols.c_x)
        QXZoo.NCU.apply_ncu!(cct2, [4,7], [], 8, DefaultGates.GateSymbols.c_x)
        QXZoo.NCU.apply_ncu!(cct2, [3,6], [], 7, DefaultGates.GateSymbols.c_x)
        QXZoo.NCU.apply_ncu!(cct2, [1,2], [], 6, DefaultGates.GateSymbols.c_x)
        QXZoo.NCU.apply_ncu!(cct2, [3,6], [], 7, DefaultGates.GateSymbols.c_x)
        QXZoo.NCU.apply_ncu!(cct2, [4,7], [], 8, DefaultGates.GateSymbols.c_x)
        QXZoo.NCU.apply_ncu!(cct2, [5,8], [], 9, DefaultGates.GateSymbols.c_x)
        QXZoo.NCU.apply_ncu!(cct2, [4,7], [], 8, DefaultGates.GateSymbols.c_x)
        QXZoo.NCU.apply_ncu!(cct2, [3,6], [], 7, DefaultGates.GateSymbols.c_x)
        QXZoo.NCU.apply_ncu!(cct2, [1,2], [], 6, DefaultGates.GateSymbols.c_x)
        QXZoo.NCU.apply_ncu!(cct2, [3,6], [], 7, DefaultGates.GateSymbols.c_x)
        QXZoo.NCU.apply_ncu!(cct2, [4,7], [], 8, DefaultGates.GateSymbols.c_x)
        cct2 << DefaultGates.h(9)

        for i in 1:5
            cct2 << DefaultGates.x(i)
            cct2 << DefaultGates.h(i)
        end
        cct2 << DefaultGates.x(9)
        cct2 << DefaultGates.h(9)

        @test cct1.circ_ops.len == cct2.circ_ops.len

        @test cct1.gate_set == cct2.gate_set

        h1 = cct1.circ_ops.head;
        h2 = cct2.circ_ops.head;
        for i in 1:cct1.circ_ops.len
            @test h1.data === h2.data
            h1 = h1.next
            h2 = h2.next
        end
    end

    @testset "Test Grover combined" begin
        cct = QXZoo.Circuit.Circ(5)
        indices = collect(1:cct.num_qubits)
        QXZoo.Grover.run_grover!(cct, indices, 11)
    end
end


text:
function _fetch(repodir::String, url::String; verbose = false)

    !_check_gitdir(repodir) && return false

    rhash = _remote_HEAD_hash(url)
    isempty(rhash) && return false

    # do fetch
    _run("git -C $(repodir) fetch 2>&1"; verbose, ignorestatus = true)

    # check success
    # TODO: add a ceck here
    return true
end
_fetch(gl::GitLink; verbose = false) = _fetch(repo_dir(gl), remote_url(gl); verbose)


text:
using StackOverflow
using Test

qHolder = getrecentquestionsfortag()
@test_nowarn qHolder[1].answer_count

@test_nowarn getrecentbadgerecipients()

siteinfoholder = getsiteinfo()
@test_nowarn dump(siteinfoholder)


text:
#=
############################################################################################
using Soss: Soss, ConditionalModel
import Soss: Soss, predict, simulate

############################################################################################
"Wrapper functions to work with Soss Models"
#=
Note: Some improvements to work on:
    -> In Objective (if multiple sampler are used to target only subset of parameter):
		When Objective is formed, make posterior such that only tagged parameter are evaluated
			-> would need to change old model.id in that case
    -> Make Methods for:
        simulate
        generate
        Predict -> predict(m(), (namedtuplevals))
	-> For Sequential Estimation methods (SMC and beyond)
        - Need to update data in SOSS model separately (when created as input?)
        - Need a way to separate log-prior and log-likelihood
=#

############################################################################################
# ModelWrapper part
"""
$(SIGNATURES)
Best guess for Soss Model parameter, excluding data and hyperparameter. Not exported.

# Examples
```julia
```

"""
function _guess_soss_param(soss_posterior::M) where {M<:Soss.ConditionalModel}
    ## All parameter from posterior.model.dists, except data (Hyperparameter should be fixed and have no model.dists entry)
    param = setdiff(keys(soss_posterior.model.dists), keys(soss_posterior.obs))
    ## Check if all param Symbols have values assigned
    ArgCheck.@argcheck all(haskey(soss_posterior.argvals, sym) for sym in param) "Not all posterior Soss model parameter have initial value assigned, please create initial value for all of them"
    ## Return NamedTuple with initial parameter
    return subset(soss_posterior.argvals, tuple(param...))
end

############################################################################################
function ModelWrapper(
    soss_posterior::M, flattendefault::F=FlattenDefault()
) where {M<:Soss.ConditionalModel,F<:FlattenDefault}
    ## Check if Posterior Soss Model provided
    ArgCheck.@argcheck !isempty(soss_posterior.obs) "No posterior Soss model provided, please use: posterior = MyModel | (MyDataName = MyDataValues,)"
    ## Guess all parameter
    val = _guess_soss_param(soss_posterior)
    ## Create prior struct from val
    _prior = subset(soss_posterior.model.dists, val)
    _prior_eval = NamedTuple{keys(_prior)}(eval(_prior[sym]) for sym in keys(_prior))
    ## Create Param NamedTuple as Safetye check that all defined Soss parameter are consistent with Param syntax
    #!NOTE: This is not ideal because all information was already there, but temporarily converting (val, prior) to Param struct guarantees that ModelWrapper can handle user input.
    params = NamedTuple{keys(val)}(
        Param(val[iter], _prior_eval[iter]) for iter in eachindex(val)
    )
    ## Return ModelWrapper
    return ModelWrapper(soss_posterior, params, flattendefault)
end

#=
#!NOTE: Could use predict(posterior, vals) directly as it returns data dimension. However, simualate needed for any algorithms and definitons slightly different in packages, so we leave it blank for now.
function simulate(model::ModelWrapper{M}) where {M<:Soss.ConditionalModel}
end
=#

############################################################################################
# Objective part

#!TODO: Make posterior such that only tagged parameter are evaluated
function Objective(model::ModelWrapper{M}) where {M<:Soss.ConditionalModel}
    return Objective(model::ModelWrapper{M}, nothing, Tagged(model))
end
function Objective(
    model::ModelWrapper{M}, tagged::T
) where {M<:Soss.ConditionalModel,T<:Tagged}
    return Objective(model::ModelWrapper{M}, nothing, tagged)
end

function (objective::Objective{<:ModelWrapper{M}})(
    θ::NamedTuple
) where {M<:Soss.ConditionalModel}
    return Soss.logdensity(objective.model.id(θ))
end

#!TODO: Need a way to update logposterior (model.id) with new data
#=

=#

#=
#NOTE: Waiting for: https://github.com/cscherrer/Soss.jl/issues/301
function predict(_rng::Random.AbstractRNG, objective::Objective{<:ModelWrapper{M}}) where {M<:Soss.ConditionalModel}
	return nothing
end
function generate(_rng::Random.AbstractRNG, objective::Objective{<:ModelWrapper{M}}) where {M<:Soss.ConditionalModel}
	return nothing
end
=#

############################################################################################
# Export
export ModelWrapper
=#


text:
# 2D diffusion problem

# Packages and inclusions
using ModelingToolkit,DiffEqOperators,LinearAlgebra,Test,OrdinaryDiffEq, DomainSets
using ModelingToolkit: Differential

# Tests
@testset "Test 00: Dt(u(t,x,y)) ~ Dxx(u(t,x,y)) + Dyy(u(t,x,y))" begin

    # Variables, parameters, and derivatives
    @parameters t x y
    @variables u(..)
    Dxx = Differential(x)^2
    Dyy = Differential(y)^2
    Dt = Differential(t)
    t_min= 0.
    t_max = 2.0
    x_min = 0.
    x_max = 2.
    y_min = 0.
    y_max = 2.

    # Analytic solution
    analytic_sol_func(t,x,y) = exp(x+y)*cos(x+y+4t)

    # Equation
    eq  = Dt(u(t,x,y)) ~ Dxx(u(t,x,y)) + Dyy(u(t,x,y))

    # Initial and boundary conditions
    bcs = [u(t_min,x,y) ~ analytic_sol_func(t_min,x,y),
           u(t,x_min,y) ~ analytic_sol_func(t,x_min,y),
           u(t,x_max,y) ~ analytic_sol_func(t,x_max,y),
           u(t,x,y_min) ~ analytic_sol_func(t,x,y_min),
           u(t,x,y_max) ~ analytic_sol_func(t,x,y_max)]

    # Space and time domains
    domains = [t ∈ Interval(t_min,t_max),
            x ∈ Interval(x_min,x_max),
            y ∈ Interval(y_min,y_max)]

    # Space and time domains
    @named pdesys = PDESystem([eq],bcs,domains,[t,x,y],[u(t,x,y)])

    dx = 0.1; dy = 0.2

    # Test against exact solution
    Nx = floor(Int64, (x_max - x_min) / dx) + 1
    Ny = floor(Int64, (y_max - y_min) / dy) + 1
    @variables u[1:Nx,1:Ny](t)
    r_space_x = x_min:dx:x_max
    r_space_y = y_min:dy:y_max
    asf = reshape([analytic_sol_func(t_max,r_space_x[i],r_space_y[j]) for j in 1:Ny for i in 1:Nx],(Nx,Ny))
    # Method of lines discretization
    order = 2
    discretization = MOLFiniteDifference([x=>dx,y=>dy],t;centered_order=order)
    prob = ModelingToolkit.discretize(pdesys,discretization)

    # Solution of the ODE system
    sol = solve(prob,Tsit5())

    # Test against exact solution
    sol′ = reshape([sol[u[(i-1)*Ny+j]][end] for i in 1:Nx for j in 1:Ny],(Nx,Ny))
    @test asf ≈ sol′ atol=0.4

    #Plot
    #using Plots
    #heatmap(sol′)
    #savefig("MOL_Linear_Diffusion_2D_Test00.png")
end

@testset "Test 01: Dt(u(t,x,y)) ~ Dx( a(x,y,u) * Dx(u(t,x,y))) + Dy( a(x,y,u) * Dy(u(t,x,y)))" begin

    # Variables, parameters, and derivatives
    @parameters t x y
    @variables u(..)
    Dx = Differential(x)
    Dy = Differential(y)
    Dt = Differential(t)
    t_min= 0.
    t_max = 2.0
    x_min = 0.
    x_max = 2.
    y_min = 0.
    y_max = 2.

    # Analytic solution
    analytic_sol_func(t,x,y) = exp(x+y)*cos(x+y+4t)

    # Equation
    eq = Dt(u(t,x,y)) ~ Dx( (u(t,x,y)^2 / exp(x+y)^2 + sin(x+y+4t)^2)^0.5 * Dx(u(t,x,y))) +
                        Dy( (u(t,x,y)^2 / exp(x+y)^2 + sin(x+y+4t)^2)^0.5 * Dy(u(t,x,y)))

    # Initial and boundary conditions
    bcs = [u(t_min,x,y) ~ analytic_sol_func(t_min,x,y),
           u(t,x_min,y) ~ analytic_sol_func(t,x_min,y),
           u(t,x_max,y) ~ analytic_sol_func(t,x_max,y),
           u(t,x,y_min) ~ analytic_sol_func(t,x,y_min),
           u(t,x,y_max) ~ analytic_sol_func(t,x,y_max)]

    # Space and time domains
    domains = [t ∈ Interval(t_min,t_max),
            x ∈ Interval(x_min,x_max),
            y ∈ Interval(y_min,y_max)]

    # Space and time domains
    @named pdesys = PDESystem([eq],bcs,domains,[t,x,y],[u(t,x,y)])

    # Method of lines discretization
    dx = 0.1; dy = 0.2
    discretization = MOLFiniteDifference([x=>dx,y=>dy],t)
    prob = ModelingToolkit.discretize(pdesys,discretization)

    # Solution of the ODE system
    sol = solve(prob,Rosenbrock23())

    # Test against exact solution
    Nx = floor(Int64, (x_max - x_min) / dx) + 1
    Ny = floor(Int64, (y_max - y_min) / dy) + 1
    @variables u[1:Nx,1:Ny](t)
    sol′ = reshape([sol[u[(i-1)*Ny+j]][end] for i in 1:Nx for j in 1:Ny],(Nx,Ny))
    r_space_x = x_min:dx:x_max
    r_space_y = y_min:dy:y_max
    asf = reshape([analytic_sol_func(t_max,r_space_x[i],r_space_y[j]) for j in 1:Ny for i in 1:Nx],(Nx,Ny))
    m = max(asf...)
    @test asf / m ≈ sol′ / m  atol=0.4 # TODO: use lower atol when g(x) is improved in MOL_discretize.jl

    #Plot
    #using Plots
    #heatmap(sol′)
    #savefig("MOL_NonLinear_Diffusion_2D_Test01.png")

end


text:
using ComponentArrays
using DifferentialEquations
using UnPack: @unpack

tspan = (0.0, 20.0)

## Lorenz system
function lorenz!(D, u, p, t; f=0.0)
    @unpack σ, ρ, β = p
    @unpack x, y, z = u

    D.x = σ*(y - x)
    D.y = x*(ρ - z) - y - f
    D.z = x*y - β*z
    return nothing
end
function lorenz_jac!(D, u, p, t)
    @unpack σ, ρ, β = p
    @unpack x, y, z = u

    D[:x,:x] = -σ
    D[:x,:y] = σ

    D[:y,:x] = ρ
    D[:y,:y] = -1
    D[:y,:z] = -x

    D[:z,:x] = y
    D[:z,:y] = x
    D[:z,:z] = -β
    return nothing
end

lorenz_p = (σ=10.0, ρ=28.0, β=8/3)
lorenz_ic = ComponentArray(x=0.0, y=0.0, z=0.0)
lorenz_fun = ODEFunction(lorenz!, jac=lorenz_jac!)
lorenz_prob = ODEProblem(lorenz_fun, lorenz_ic, tspan, lorenz_p)

## Lotka-Volterra system
function lotka!(D, u, p, t; f=0.0)
    @unpack α, β, γ, δ = p
    @unpack x, y = u

    D.x =  α*x - β*x*y + f
    D.y = -γ*y + δ*x*y
    return nothing
end
function lotka_jac!(D, u, p, t)
    @unpack α, β, γ, δ = p
    @unpack x, y = u

    D[:x,:x] = α - β*y
    D[:x,:y] = -β*x

    D[:y,:x] = δ*y
    D[:y,:y] = -γ + δ*x
    return nothing
end

lotka_p = (α=2/3, β=4/3, γ=1.0, δ=1.0)
lotka_ic = ComponentArray(x=1.0, y=1.0)
lotka_fun = ODEFunction(lotka!, jac=lotka_jac!)
lotka_prob = ODEProblem(lotka_fun, lotka_ic, tspan, lotka_p)

## Composed Lorenz and Lotka-Volterra system
function composed!(D, u, p, t)
    c = p.c #coupling parameter
    @unpack lorenz, lotka = u

    lorenz!(D.lorenz, lorenz, p.lorenz, t, f=c*lotka.x)
    lotka!(D.lotka, lotka, p.lotka, t, f=c*lorenz.x)
    return nothing
end
function composed_jac!(D, u, p, t)
    c = p.c
    @unpack lorenz, lotka = u

    lorenz_jac!(@view(D[:lorenz,:lorenz]), lorenz, p.lorenz, t)
    lotka_jac!(@view(D[:lotka,:lotka]), lotka, p.lotka, t)

    @view(D[:lorenz,:lotka])[:y,:x] = -c
    @view(D[:lotka,:lorenz])[:x,:x] = c
    return nothing
end

comp_p = (lorenz=lorenz_p, lotka=lotka_p, c=0.01)
comp_ic = ComponentArray(lorenz=lorenz_ic, lotka=lotka_ic)
comp_fun = ODEFunction(composed!, jac=composed_jac!)
comp_prob = ODEProblem(comp_fun, comp_ic, tspan, comp_p)

## Solve problem
# We can solve the composed system...
comp_sol = solve(comp_prob, Rodas5())

# ...or we can unit test one of the component systems
lotka_sol = solve(lotka_prob, Rodas5())


text:
if haskey(ENV, "DOCKER_BUILD")
    include("build-docker.jl")
else
    include("build-standard.jl")
end


text:
lesser(x::TimeOrderedGreenFunction) = x.L

function lesser(x::TimeOrderedConvolution{TimeOrderedGreenFunction, TimeOrderedGreenFunction})
    a, b = x.A, x.B
    a.R * b.L + a.L * b.A
end

function lesser(x::TimeOrderedConvolution{TimeOrderedGreenFunction, TimeOrderedConvolution{TimeOrderedGreenFunction,TimeOrderedGreenFunction}})
    a, b, c = x.A, x.B.A, x.B.B
    lesser3(a, b, c)
end

function lesser(x::TimeOrderedConvolution{TimeOrderedConvolution{TimeOrderedGreenFunction,TimeOrderedGreenFunction}, TimeOrderedGreenFunction})
    a, b, c = x.A.A, x.A.B, x.B
    lesser3(a, b, c)
end

function lesser3(a, b, c)
    x = a.L * b.A * c.A
    x - x' + a.R * b.L * c.A
end

greater(x::TimeOrderedGreenFunction) = x.G

function greater(x::TimeOrderedConvolution{TimeOrderedGreenFunction, TimeOrderedGreenFunction})
    a, b = x.A, x.B
    a.R * b.G + a.G * b.A
end

function greater(x::TimeOrderedConvolution{TimeOrderedGreenFunction, TimeOrderedConvolution{TimeOrderedGreenFunction,TimeOrderedGreenFunction}})
    a, b, c = x.A, x.B.A, x.B.B
    greater3(a, b, c)
end

function greater(x::TimeOrderedConvolution{TimeOrderedConvolution{TimeOrderedGreenFunction,TimeOrderedGreenFunction}, TimeOrderedGreenFunction})
    a, b, c = x.A.A, x.A.B, x.B
    greater3(a, b, c)
end

function greater3(a, b, c)
    x = a.G * b.A * c.A
    x - x' + a.R * b.G * c.A
end


text:
# Generate 2D vertices in a rectangular Lattice
function _generate_2d_vertices!(vertices::Vector{Tensors.Vec{2,T}}, nx, ny, LL, LR, UR, UL) where {T}
  for i in 0:ny-1
    ratio_bounds = i / (ny-1)

    x0 = LL[1] * (1 - ratio_bounds) + ratio_bounds * UL[1]
    x1 = LR[1] * (1 - ratio_bounds) + ratio_bounds * UR[1]

    y0 = LL[2] * (1 - ratio_bounds) + ratio_bounds * UL[2]
    y1 = LR[2] * (1 - ratio_bounds) + ratio_bounds * UR[2]

    for j in 0:nx-1
        ratio = j / (nx-1)
        x = x0 * (1 - ratio) + ratio * x1
        y = y0 * (1 - ratio) + ratio * y1
        push!(vertices, Tensors.Vec{2}((x, y)))
    end
end
end

# Shortcuts
unitSquareMesh2(cellType, nel::NTuple{2,Int}) = rectangle_mesh2(cellType, nel, Tensors.Vec{2}((0.0,0.0)), Tensors.Vec{2}((1.0,1.0)))

#########################
# Triangle Cells 2D   #
#########################
"""
rectangle_mesh2(::Type{TriangleCell}, nel::NTuple{2,Int}, LL::Vec{2,T}, UR::Vec{2,T})
Generate a rectangular mesh with triangular cells, where `LL` is the low left vertex
and `UR` is the upper right one. `nel` is a tuple with the number of partions to be
used in each dimension.
"""
function rectangle_mesh2(::Type{TriangleCell}, nel::NTuple{2,Int}, LL::Tensors.Vec{2,T}, UR::Tensors.Vec{2,T}) where {T}
  LR = Tensors.Vec{2}((UR[1],LL[2]))
  UL = Tensors.Vec{2}((LL[1],UR[2]))
  nel_x = nel[1]; nel_y = nel[2]; nel_tot = 2*nel_x*nel_y
  n_vertices_x = nel_x + 1; n_vertices_y = nel_y + 1
  n_vertices = n_vertices_x * n_vertices_y

  # Generate vertices
  vertices = Tensors.Vec{2,Float64}[]
  _generate_2d_vertices!(vertices, n_vertices_x, n_vertices_y, LL, LR, UR, UL)
  typeof(vertices)

  # Generate cells
  vertex_array = reshape(collect(1:n_vertices), (n_vertices_x, n_vertices_y))
  offsets = [1]
  k = 1
  indices = []
  for j in 1:nel_y, i in 1:nel_x
      k = k + 3
      push!(offsets,k)
      push!(indices,(vertex_array[i,j], vertex_array[i+1,j], vertex_array[i,j+1])...) # ◺
      k = k + 3
      push!(offsets,k)
      push!(indices, (vertex_array[i+1,j], vertex_array[i+1,j+1], vertex_array[i,j+1])...) # ◹
  end

  geometry = Dict((2,0) => MeshConectivity(Tuple(indices), Tuple(offsets)))
  entities = (n_vertices, nel_x+nel_y+3*nel_x*nel_y, nel_tot)

  # Compute boundaries
  cell_array = reshape(collect(1:nel_tot),(2, nel_x, nel_y))
  boundary = Tuple{Int,Int}[[(cl, 1) for cl in cell_array[1,:,1]];
                               [(cl, 1) for cl in cell_array[2,end,:]];
                               [(cl, 2) for cl in cell_array[2,:,end]];
                               [(cl, 3) for cl in cell_array[1,1,:]]]
  offset = 0
  edgesets = Dict{String,Set{NTuple{2,Int}}}()
  edgesets["bottom"] = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[1,:,1]))   .+ offset]); offset += length(cell_array[1,:,1])
  edgesets["right"]  = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[2,end,:])) .+ offset]); offset += length(cell_array[2,end,:])
  edgesets["top"]    = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[2,:,end])) .+ offset]); offset += length(cell_array[2,:,end])
  edgesets["left"]   = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[1,1,:]))   .+ offset]); offset += length(cell_array[1,1,:])
  edgesets["boundary"] = union(edgesets["bottom"],edgesets["right"],edgesets["top"],edgesets["left"])

  return PolytopalMesh2(entities,vertices,geometry,Dict(1=>edgesets))
end

#########################
# Rectangle Cells 2D   #
#########################
"""
rectangle_mesh2(::Type{RectangleCell}, nel::NTuple{2,Int}, LL::Vec{2,T}, UR::Vec{2,T}) where {T}
Generate a rectangular mesh with triangular cells, where `LL` is the low left vertex
and `UR` is the upper right one. `nel` is a tuple with the number of partions to be
used in each dimension.
"""
function rectangle_mesh2(::Type{RectangleCell}, nel::NTuple{2,Int}, LL::Tensors.Vec{2,T}, UR::Tensors.Vec{2,T}) where {T}
  LR = Tensors.Vec{2}((UR[1],LL[2]))
  UL = Tensors.Vec{2}((LL[1],UR[2]))
  nel_x = nel[1]; nel_y = nel[2]; nel_tot = nel_x*nel_y
  n_vertices_x = nel_x + 1; n_vertices_y = nel_y + 1
  n_vertices = n_vertices_x * n_vertices_y

  # Generate vertices
  vertices = Tensors.Vec{2,Float64}[]
  _generate_2d_vertices!(vertices, n_vertices_x, n_vertices_y, LL, LR, UR, UL)

  # Generate cells
  vertex_array = reshape(collect(1:n_vertices), (n_vertices_x, n_vertices_y))
  offsets = [1]
  k = 1
  indices = []
  for j in 1:nel_y, i in 1:nel_x
    k = k + 4
    push!(offsets,k)
    push!(indices,(vertex_array[i,j], vertex_array[i+1,j], vertex_array[i+1,j+1], vertex_array[i,j+1])...)
  end

  geometry = Dict((2,0) => MeshConectivity(Tuple(indices), Tuple(offsets)))
  entities = (n_vertices, nel_x+nel_y+2*nel_x*nel_y, nel_tot)

  # Cell faces
  cell_array = reshape(collect(1:nel_tot),(nel_x, nel_y))
  boundary = Tuple{Int,Int}[[(cl, 1) for cl in cell_array[:,1]];
                            [(cl, 2) for cl in cell_array[end,:]];
                            [(cl, 3) for cl in cell_array[:,end]];
                            [(cl, 4) for cl in cell_array[1,:]]]
  # Cell face sets
  offset = 0
  edgesets = Dict{String, Set{Tuple{Int,Int}}}()
  edgesets["bottom"] = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[:,1]))   .+ offset]); offset += length(cell_array[:,1])
  edgesets["right"]  = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[end,:])) .+ offset]); offset += length(cell_array[end,:])
  edgesets["top"]    = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[:,end])) .+ offset]); offset += length(cell_array[:,end])
  edgesets["left"]   = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[1,:]))   .+ offset]); offset += length(cell_array[1,:])
  edgesets["boundary"] = union(edgesets["bottom"],edgesets["right"],edgesets["top"],edgesets["left"])
  return PolytopalMesh2(entities,vertices,geometry,Dict(1=>edgesets))
end

# #########################
# # Hexagonal Cells 2D
# #############################

function _gen_hexagon2(centroid,LL,UR,hex_width,hex_heigth)
  hex_verts = [
      Tensors.Vec{2}((centroid[1]-hex_width/2,centroid[2]-hex_heigth/4)),
      Tensors.Vec{2}((centroid[1],centroid[2]-hex_heigth/2)),
      Tensors.Vec{2}((centroid[1]+hex_width/2,centroid[2]-hex_heigth/4)),
      Tensors.Vec{2}((centroid[1]+hex_width/2,centroid[2]+hex_heigth/4)),
      Tensors.Vec{2}((centroid[1],centroid[2]+hex_heigth/2)),
      Tensors.Vec{2}((centroid[1]-hex_width/2,centroid[2]+hex_heigth/4)),
  ]
  filter_verts = filter(x->(x[1]>=LL[1] && x[1]<=UR[1]
                        && x[2]>=LL[2] && x[2] <= UR[2]),hex_verts)
  if size(filter_verts,1) == 2 #corners case
      x = centroid[1] + (centroid[1] == LL[1] ? hex_width/2 : -hex_width/2)
      new_vertex = Tensors.Vec{2}((x, centroid[2]))
      if centroid == LL
          filter_verts = (centroid,new_vertex,filter_verts...)
      elseif centroid[2] == LL[2] && centroid[1] == UR[1]
          filter_verts = (new_vertex,centroid,filter_verts...)
      elseif centroid[2] == UR[2] && centroid[1] == LL[1]
          filter_verts = (filter_verts...,new_vertex,centroid)
      elseif centroid == UR
          filter_verts = (filter_verts...,centroid,new_vertex)
      else
          throw("error on hexagon at corner $centroid")
      end
  elseif size(filter_verts,1) == 3 #top and bottom pentagon/triangle case
      nv1 = Tensors.Vec{2}((centroid[1]-hex_width/2, centroid[2]))
      nv2 = Tensors.Vec{2}((centroid[1]+hex_width/2, centroid[2]))
      if centroid[2] == LL[2]
          filter_verts = (nv1,nv2,filter_verts...)
      elseif centroid[2] == UR[2]
          filter_verts = (nv1,filter_verts...,nv2)
      elseif centroid[2] > UR[2]
          filter_verts = (filter_verts...,)
      else
          throw("error at bottom centroid $centroid")
      end
  end
  return filter_verts
end

function rectangle_mesh2(::Type{HexagonCell}, nel::NTuple{2,Int}, LL::Tensors.Vec{2,T}, UR::Tensors.Vec{2,T}) where {T}
  LR = Tensors.Vec{2}((UR[1],LL[2]))
  UL = Tensors.Vec{2}((LL[1],UR[2]))
  nel_x = nel[1]; nel_y = nel[2]

  # Generate vertices
  vertices = Tensors.Vec{2,Float64}[]
  hex_width = 1.0
  hex_heigth = 1.0
  w_scale = (LR[1] - LL[1])/nel_x
  h_scale = (UL[2] - LL[1])/nel_y
  LLn = Tensors.Vec{2}((0.0,0.0))
  URn = Tensors.Vec{2}((nel_x,nel_y))

  map_func(x::Tensors.Vec{2,T}) where T = Tensors.Vec{2}((x[1]*w_scale+LL[1],x[2]*h_scale+LL[2]))

  centroids = Tensors.Vec{2,T}[]
  _generate_2d_hex_centroids!(centroids,LL, nel_y, nel_x, hex_width, hex_heigth)
  used_vertices = Dict{Tensors.Vec{2,T},Int}()
  #Add cells
  offsets = [1]
  k = 1
  indices = []
  nextvert = 1 # next free vertex to use
  for c_i in 1:size(centroids,1)
          cell_verts = Int[]
          for vert in _gen_hexagon2(centroids[c_i],LLn,URn,hex_width,hex_heigth)
              nextvert = _push_cell_vertex!(vert,cell_verts, used_vertices,vertices,nextvert, map_func)
          end
          k = k + size(cell_verts,1)
          push!(offsets,k)
          push!(indices,cell_verts...)
  end
  nel_tot = size(centroids,1)
  n_vertices = size(vertices,1)
  geometry = Dict((2,0) => MeshConectivity(Tuple(indices), Tuple(offsets)))

  # Cell edges
  n_edges = 0
  edgesets = Dict{String,Set{NTuple{2,Int}}}()
  t_edge = mod1(nel_y,3)==1 ? 4 : 3
  t_idxs = (nel_tot-nel_x + (mod1(nel_y,3)==3 ? 0 : 1)):nel_tot
  l_bounds = Tuple{Int,Int}[]
  r_bounds = Tuple{Int,Int}[]
  c_idx = 1
  for i in 1:nel_y
      if mod1(i,3) == 1
          push!(l_bounds,(c_idx,4))
          push!(r_bounds,(c_idx+nel_x,2))
          ed_idx = (i == nel_y ? 5 : 6)
          push!(l_bounds,(c_idx+nel_x+1,ed_idx))
          push!(r_bounds,(c_idx+2*nel_x,3))
          c_idx = c_idx + nel_x*2+1
          n_edges = n_edges + 3 + 6*nel_x
      elseif mod1(i,3) == 2
          push!(l_bounds,(c_idx,4))
          push!(r_bounds,(c_idx+nel_x,2))
          c_idx = c_idx + nel_x + 1
          n_edges = n_edges + 1 + 6*nel_x
      else
          push!(l_bounds,(c_idx,6))
          push!(r_bounds,(c_idx+nel_x-1,3))
          if i == nel_y
              push!(l_bounds,(c_idx+nel_x,4))
              push!(r_bounds,(c_idx+2*nel_x,2))
          end
          c_idx = c_idx + nel_x
          n_edges = n_edges + 3 + 5*nel_x
      end
  end
  edgesets["bottom"] = Set{Tuple{Int,Int}}([(c_i,1) for c_i in 1:(nel_x+1)])
  edgesets["top"]    = Set{Tuple{Int,Int}}([(i,t_edge) for i in t_idxs])
  edgesets["right"]  = Set{Tuple{Int,Int}}([(x[1],x[2]) for x in r_bounds])
  edgesets["left"]   = Set{Tuple{Int,Int}}([(x[1],x[2]) for x in l_bounds])
  edgesets["boundary"] = union(edgesets["bottom"],edgesets["right"],edgesets["top"],edgesets["left"])

  entities = (n_vertices, n_edges, nel_tot)

  return PolytopalMesh2(entities,vertices,geometry,Dict(1=>edgesets))
end

##########################
# 3D Data
##########################
function _build_hexahedron_geometry!(geometry,indices, offsets, nel, nedges,nfaces)
    indices1 = [Int[] for _ in 1:nel]
    indices2 = [Int[] for _ in 1:nfaces]
    indices3 = [Int[] for _ in 1:nel]
    indices4 = [Int[] for _ in 1:nedges]
    facesDict = Dict{Set{Int},Int}()
    edgesDict = Dict{Set{Int},Int}()
    nextidx = 1 #first face index
    next_e_idx = 1 #first face index
    for cell in 1:nel
        nodes = indices[offsets[cell]:offsets[cell+1]-1]
        # Compute faces geometries
        Vi = [nodes[1:4],nodes[5:8],nodes[[2,3,6,7]],nodes[[1,5,8,4]],nodes[[3,4,8,7]],nodes[[1,2,6,5]]]
        for vi in Vi
            token = Base.ht_keyindex2!(facesDict, Set(vi))
            if token > 0 # reuse edge index
            reuse_idx = facesDict.vals[token]
            push!(indices1[cell], reuse_idx)
            else
            push!(indices1[cell], nextidx)
            indices2[nextidx] = vi
            Base._setindex!(facesDict, nextidx, Set(vi), -token)
            nextidx += 1
            end
        end
        # Compute edges geometries
        Vi = [nodes[1:2],nodes[2:3],nodes[3:4],nodes[[4,1]],nodes[5:6],nodes[6:7],
              nodes[7:8],nodes[[8,5]],nodes[[1,4]],nodes[[2,3]],nodes[[6,7]],nodes[[5,8]]]
        for vi in Vi
            token = Base.ht_keyindex2!(edgesDict, Set(vi))
            if token > 0 # reuse edge index
            reuse_idx = edgesDict.vals[token]
            push!(indices3[cell], reuse_idx)
            else
            push!(indices3[cell], next_e_idx)
            indices4[next_e_idx] = vi
            Base._setindex!(edgesDict, next_e_idx, Set(vi), -token)
            next_e_idx += 1
            end
        end
    end
    V = _pack_connectivity(indices1)
    push!(geometry,(3,2)=>V)
    V = _pack_connectivity(indices2)
    push!(geometry,(2,0)=>V)
    V = _pack_connectivity(indices3)
    push!(geometry,(3,1)=>V)
    V = _pack_connectivity(indices4)
    push!(geometry,(1,0)=>V)
  end

  function _build_tetrahedron_geometry!(geometry,indices, offsets, nel, nedges,nfaces)
    indices1 = [Int[] for _ in 1:nel]
    indices2 = [Int[] for _ in 1:nfaces]
    indices3 = [Int[] for _ in 1:nel]
    indices4 = [Int[] for _ in 1:nedges]
    facesDict = Dict{Set{Int},Int}()
    edgesDict = Dict{Set{Int},Int}()
    nextidx = 1 #first face index
    next_e_idx = 1 #first face index
    for cell in 1:nel
        nodes = indices[offsets[cell]:offsets[cell+1]-1]
        # Compute faces geometries
        Vi = [nodes[1:3],nodes[2:4],nodes[[1,3,4]],nodes[[1,2,4]]]
        for vi in Vi
            token = Base.ht_keyindex2!(facesDict, Set(vi))
            if token > 0 # reuse edge index
            reuse_idx = facesDict.vals[token]
            push!(indices1[cell], reuse_idx)
            else
            push!(indices1[cell], nextidx)
            indices2[nextidx] = vi
            Base._setindex!(facesDict, nextidx, Set(vi), -token)
            nextidx += 1
            end
        end
        # Compute edges geometries
        N = size(nodes,1)
        Vi = [[1,2],[2,3],[3,4],[4,1],[1,3],[2,4]]
        for vi in Vi
            token = Base.ht_keyindex2!(edgesDict, Set(vi))
            if token > 0 # reuse edge index
            reuse_idx = edgesDict.vals[token]
            push!(indices3[cell], reuse_idx)
            else
            push!(indices3[cell], next_e_idx)
            indices4[next_e_idx] = vi
            Base._setindex!(edgesDict, next_e_idx, Set(vi), -token)
            next_e_idx += 1
            end
        end
    end
    V = _pack_connectivity(indices1)
    push!(geometry,(3,2)=>V)
    V = _pack_connectivity(indices2)
    push!(geometry,(2,0)=>V)
    V = _pack_connectivity(indices3)
    push!(geometry,(3,1)=>V)
    V = _pack_connectivity(indices4)
    push!(geometry,(1,0)=>V)
  end

#####################################
# Hexahedron 3D
#####################################
function hyper_rectagle_mesh2(::Type{HexahedronCell}, nel::NTuple{3,Int}, left::Vec{3,T}=Vec{3}((0.0,0.0,0.0)), right::Vec{3,T}=Vec{3}((1.0,1.0,1.0))) where {T}
    nel_x = nel[1]; nel_y = nel[2]; nel_z = nel[3]; nel_tot = nel_x*nel_y*nel_z
    n_nodes_x = nel_x + 1; n_nodes_y = nel_y + 1; n_nodes_z = nel_z + 1
    n_nodes = n_nodes_x * n_nodes_y * n_nodes_z

    # Generate vertices
    coords_x = range(left[1], stop=right[1], length=n_nodes_x)
    coords_y = range(left[2], stop=right[2], length=n_nodes_y)
    coords_z = range(left[3], stop=right[3], length=n_nodes_z)
    nodes = Tensors.Vec{3,T}[]
    for k in 1:n_nodes_z, j in 1:n_nodes_y, i in 1:n_nodes_x
        push!(nodes, Tensors.Vec{3}((coords_x[i], coords_y[j], coords_z[k])))
    end

    # Generate cells
    node_array = reshape(collect(1:n_nodes), (n_nodes_x, n_nodes_y, n_nodes_z))
    offsets = [1]
    l = 1
    indices = Int[]
    for k in 1:nel_z, j in 1:nel_y, i in 1:nel_x
        l = l + 8
        push!(offsets, l)
        # cube = (1, 2, 3, 4, 5, 6, 7, 8)
        # left = (1, 4, 5, 8), right = (2, 3, 6, 7)
        # front = (1, 2, 5, 6), back = (3, 4, 7, 8)
        # bottom = (1, 2, 3, 4), top = (5, 6, 7, 8)
        push!(indices, (node_array[i,j,k], node_array[i+1,j,k], node_array[i+1,j+1,k], node_array[i,j+1,k],
                                node_array[i,j,k+1], node_array[i+1,j,k+1], node_array[i+1,j+1,k+1], node_array[i,j+1,k+1])...)
    end

    n_edges = (nel_x+nel_y+2*nel_x*nel_y)*(nel_z+1) + (nel_y+1)*nel_z*(nel_x+1)
    n_faces = nel_x*nel_y*(nel_z+1)+nel_y*nel_z*(nel_x+1)+nel_x*nel_z*(nel_y+1)
    entities = (n_nodes, n_edges, n_faces,nel_tot)

    geometry = Dict((3,0) => MeshConectivity(Tuple(indices), Tuple(offsets)))
    _build_hexahedron_geometry!(geometry, indices, offsets,nel_tot, n_edges, n_faces)

    # Cell faces
    cell_array = reshape(collect(1:nel_tot),(nel_x, nel_y, nel_z))
    boundary = Tuple{Int,Int}[[(cl, 1) for cl in cell_array[:,:,1][:]];
                              [(cl, 2) for cl in cell_array[:,1,:][:]];
                              [(cl, 3) for cl in cell_array[end,:,:][:]];
                              [(cl, 4) for cl in cell_array[:,end,:][:]];
                              [(cl, 5) for cl in cell_array[1,:,:][:]];
                              [(cl, 6) for cl in cell_array[:,:,end][:]]]

    # Cell face sets
    offset = 0
    facesets = Dict{String,Set{Tuple{Int,Int}}}()
    facesets["bottom"] = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[:,:,1][:]))   .+ offset]); offset += length(cell_array[:,:,1][:])
    facesets["front"]  = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[:,1,:][:]))   .+ offset]); offset += length(cell_array[:,1,:][:])
    facesets["right"]  = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[end,:,:][:])) .+ offset]); offset += length(cell_array[end,:,:][:])
    facesets["back"]   = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[:,end,:][:])) .+ offset]); offset += length(cell_array[:,end,:][:])
    facesets["left"]   = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[1,:,:][:]))   .+ offset]); offset += length(cell_array[1,:,:][:])
    facesets["top"]    = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[:,:,end][:])) .+ offset]); offset += length(cell_array[:,:,end][:])
    facesets["boundary"] = union(facesets["bottom"],facesets["right"],facesets["top"],facesets["left"],facesets["back"],facesets["front"])

    return PolytopalMesh2(entities,nodes,geometry,Dict(2=>facesets))
end

#########################33
# Tetrahedron
##########################33

function hyper_rectagle_mesh2(::Type{TetrahedronCell}, nel::NTuple{3,Int}, left::Vec{3,T}=Vec{3}((0.0,0.0,0.0)), right::Vec{3,T}=Vec{3}((1.0,1.0,1.0))) where {T}
    nodes_per_dim = nel .+ 1

    cells_per_cube = 6
    n_nodes = prod(nodes_per_dim)
    nel_tot = cells_per_cube * prod(nel)

    n_nodes_x, n_nodes_y, n_nodes_z = nodes_per_dim
    nel_x, nel_y, nel_z = nel

    # Generate nodes
    coords_x = range(left[1], stop=right[1], length=n_nodes_x)
    coords_y = range(left[2], stop=right[2], length=n_nodes_y)
    coords_z = range(left[3], stop=right[3], length=n_nodes_z)
    numbering = reshape(1:n_nodes, nodes_per_dim)

    # Pre-allocate the nodes & cells
    nodes = Vector{Tensors.Vec{3,T}}(undef, n_nodes)

    # Generate nodes
    node_idx = 1
    @inbounds for k in 1:n_nodes_z, j in 1:n_nodes_y, i in 1:n_nodes_x
        nodes[node_idx] = Tensors.Vec{3}((coords_x[i], coords_y[j], coords_z[k]))
        node_idx += 1
    end

    # Generate cells, case 1 from: http://www.baumanneduard.ch/Splitting%20a%20cube%20in%20tetrahedras2.htm
    # cube = (1, 2, 3, 4, 5, 6, 7, 8)
    # left = (1, 4, 5, 8), right = (2, 3, 6, 7)
    # front = (1, 2, 5, 6), back = (3, 4, 7, 8)
    # bottom = (1, 2, 3, 4), top = (5, 6, 7, 8)
    offsets = [1]
    l = 1
    indices = []
    @inbounds for k in 1:nel_z, j in 1:nel_y, i in 1:nel_x
        cell = (
            numbering[i  , j  , k],
            numbering[i+1, j  , k],
            numbering[i+1, j+1, k],
            numbering[i  , j+1, k],
            numbering[i  , j  , k+1],
            numbering[i+1, j  , k+1],
            numbering[i+1, j+1, k+1],
            numbering[i  , j+1, k+1]
        )
        l += 4; push!(offsets, l)
        push!(indices, (cell[1], cell[2], cell[4], cell[8])...)
        l += 4; push!(offsets, l)
        push!(indices, (cell[1], cell[5], cell[2], cell[8])...)
        l += 4; push!(offsets, l)
        push!(indices, (cell[2], cell[3], cell[4], cell[8])...)
        l += 4; push!(offsets, l)
        push!(indices, (cell[2], cell[7], cell[3], cell[8])...)
        l += 4; push!(offsets, l)
        push!(indices, (cell[2], cell[5], cell[6], cell[8])...)
        l += 4; push!(offsets, l)
        push!(indices, (cell[2], cell[6], cell[7], cell[8])...)
    end

    geometry = Dict((3,0) => MeshConectivity(Tuple(indices), Tuple(offsets)))
    n_edges = (nel_x+nel_y+3*nel_x*nel_y)*(nel_z+1) + (nel_z+2*nel_z*nel_y)*(nel_x+1) + nel_x*nel_z*(nel_y+1) + nel_y*nel_x*nel_z
    n_faces = 2*(nel_x*nel_y*(nel_z+1)+nel_y*nel_z*(nel_x+1)+nel_x*nel_z*(nel_y+1)) + 6*nel_x*nel_y*nel_z
    entities = (n_nodes, n_edges, n_faces,nel_tot)
    _build_tetrahedron_geometry!(geometry, indices, offsets,nel_tot, n_edges, n_faces)

    # Order the cells as c_nxyz[n, x, y, z] such that we can look up boundary cells
    c_nxyz = reshape(1:nel_tot, (cells_per_cube, nel...))

    @views le = [map(x -> (x,4), c_nxyz[1, 1, :, :][:])   ; map(x -> (x,2), c_nxyz[2, 1, :, :][:])]
    @views ri = [map(x -> (x,1), c_nxyz[4, end, :, :][:]) ; map(x -> (x,1), c_nxyz[6, end, :, :][:])]
    @views fr = [map(x -> (x,1), c_nxyz[2, :, 1, :][:])   ; map(x -> (x,1), c_nxyz[5, :, 1, :][:])]
    @views ba = [map(x -> (x,3), c_nxyz[3, :, end, :][:]) ; map(x -> (x,3), c_nxyz[4, :, end, :][:])]
    @views bo = [map(x -> (x,1), c_nxyz[1, :, :, 1][:])   ; map(x -> (x,1), c_nxyz[3, :, :, 1][:])]
    @views to = [map(x -> (x,3), c_nxyz[5, :, :, end][:]) ; map(x -> (x,3), c_nxyz[6, :, :, end][:])]

    facesets = Dict(
        "left" => Set(le),
        "right" => Set(ri),
        "front" => Set(fr),
        "back" => Set(ba),
        "bottom" => Set(bo),
        "top" => Set(to),
        "boundary" => Set(union(le,ri,fr,ba,bo,to))
    )

    return PolytopalMesh2(entities,nodes,geometry,Dict(2=>facesets))

end

max_stars_repo_name: Paulms/jFEMToools

max_issues_repo_name: Paulms/jFEMToools

max_forks_repo_name: Paulms/jFEMToools


text:
# This file contains some extra checks that don't fall into the other
# categories, for example because they're very platform-specific.

check_os_abi(::ObjectHandle, ::Platform, rest...; verbose::Bool = false, kwargs...) = true

function check_os_abi(oh::ELFHandle, ::FreeBSD, rest...; verbose::Bool = false, kwargs...)
    if oh.ei.osabi != ELF.ELFOSABI_FREEBSD
        # The dynamic loader should not have problems in this case, but the
        # linker may not appreciate.  Let the user know about this.
        if verbose
            msg = replace("""
            $(basename(path(oh))) has an ELF header OS/ABI value that is not set to FreeBSD
            ($(ELF.ELFOSABI_FREEBSD)), this may be an issue at link time"
            """, '\n' => ' ')
            @warn(strip(msg))
        end
        return false
    end
    return true
end


text:
"""
    Entry(plottype::PlotFunc, positional::Arguments, named::NamedArguments)

Define plottype as well as positional and named arguments for a single plot.
"""
struct Entry
    plottype::PlotFunc
    positional::Arguments
    named::NamedArguments
end

function Base.get(entry::Entry, key::Int, default)
    return get(entry.positional, key, default)
end

function Base.get(entry::Entry, key::Symbol, default)
    return get(entry.named, key, default)
end

function Base.append!(e1::Entry, e2::Entry)
    plottype = assert_equal(e1.plottype, e2.plottype)
    positional = map(append!, e1.positional, e2.positional)
    named = map(append_or_assertequal!, e1.named, e2.named)
    return Entry(plottype, positional, named)
end

# Use technique from https://github.com/JuliaPlots/AlgebraOfGraphics.jl/pull/289
# to encode all axis information without creating the axis.
struct AxisSpec
    type::Union{Type{Axis}, Type{Axis3}}
    position::Tuple{Int,Int}
    attributes::NamedArguments
end

extract_type(; type=Axis, options...) = type, NamedArguments(options)

function AxisSpec(position, options)
    type, attributes = extract_type(; pairs(options)...)
    return AxisSpec(type, Tuple(position), attributes)
end

struct AxisSpecEntries
    axis::AxisSpec
    entries::Vector{Entry}
    categoricalscales::MixedArguments
    continuousscales::MixedArguments
end

"""
    AxisEntries(axis::Union{Axis, Nothing}, entries::Vector{Entry}, categoricalscales, continuousscales)

Define all ingredients to make plots on an axis.
Each categorical scale should be a `CategoricalScale`, and each continuous
scale should be a `ContinuousScale`.
"""
struct AxisEntries
    axis::Union{Axis, Axis3}
    entries::Vector{Entry}
    categoricalscales::MixedArguments
    continuousscales::MixedArguments
end

function AxisEntries(ae::AxisSpecEntries, fig)
    ax = ae.axis.type(fig[ae.axis.position...]; pairs(ae.axis.attributes)...)
    AxisEntries(ax, ae.entries, ae.categoricalscales, ae.continuousscales)
end

function AxisEntries(ae::AxisSpecEntries, ax::Union{Axis,Axis3})
    AxisEntries(ax, ae.entries, ae.categoricalscales, ae.continuousscales)
end

function Makie.plot!(ae::AxisEntries)
    axis, entries = ae.axis, ae.entries
    for entry in entries
        plot!(entry.plottype, axis, entry.positional...; pairs(entry.named)...)
    end
    return ae
end

entries(grid::AbstractMatrix{AxisEntries}) = Iterators.flatten(ae.entries for ae in grid)
entries(grid::AbstractMatrix{AxisSpecEntries}) = Iterators.flatten(ae.entries for ae in grid)

struct FigureGrid
    figure::Figure
    grid::Matrix{AxisEntries}
end

Base.show(io::IO, fg::FigureGrid) = show(io, fg.figure)
Base.show(io::IO, m::MIME, fg::FigureGrid) = show(io, m, fg.figure)
Base.show(io::IO, ::MIME"text/plain", fg::FigureGrid) = print(io, "FigureGrid()")

Base.showable(mime::MIME{M}, fg::FigureGrid) where {M} = showable(mime, fg.figure)

Base.display(fg::FigureGrid) = display(fg.figure)

function FileIO.save(filename::String, fg::FigureGrid; kwargs...)
    return FileIO.save(FileIO.query(filename), fg; kwargs...)
end

function FileIO.save(file::FileIO.Formatted, fg::FigureGrid; kwargs...)
    return FileIO.save(file, fg.figure; kwargs...)
end

to_tuple(fg) = (fg.figure, fg.grid)

Base.iterate(fg::FigureGrid) = iterate(to_tuple(fg))
Base.iterate(fg::FigureGrid, i) = iterate(to_tuple(fg), i)


text:
#=
TERMINOLOGY

Client          The side of a connection which initiates.
Server          The side of a connection which accepts.
Origin          Sender
Destination     Receiver
HTS             HTTP server
JCE             Julia Client Echoing
BCE             Browser Client Echoing

Note we use 'inverse speed' below. In lack of better words, we call this speed.
This is more compatible with BenchmarkTools and more directly useful.
For some fun reasoning, refer https://www.lesswrong.com/posts/qNiH3RRTiXXyvMJRg/inverse-speed

Note 'latency' can also be defined in other ways. Ad hoc:

speed           [ns/b]  time / amount
messagesize     [b]     Datasize for an individual message (one or more frames)
clientlatency   [ns]    Time from start sending from client to having received on destination
serverlatency   [ns]    Time from start sending from server to having received on destination
clientspeed     [ns/b] = [s/GB]
                        = upspeed = clientlatency / messagesize
                        For JCE: Client sends the message, server records its having received time
                        For BCE: Server sends the message, client sends back two messages of different lengths.
                                 Clientspeed is calculated from server time records.
serverspeed     [ns/b] = [s/GB]
                        = Downspeed = serverlatency / messagesize
                        For JCE: Server sends the message, client measures its having received time
                        For BCE: Server sends the message, client sends back two messages of different lengths.
                                 Serverspeed is calculated from server time records.
clientbandwidth [ns/b] = [s/GB]
                          Average clientspeed, time per datasize
serverbandwidth [ns/b] = [s/GB]
                         Average serverspeed, time per datasize

=#

#=
-----  TODO  ----

Not used:
** delay           [ns]    Time usefully spent between received message and start of reply
** clientRTT       [ns]    Round-trip-time for a message initiated by client
** serverRTT       [ns]    Round-trip-time for a message initiated by server
** Up2down      ClientBandwidth / ServerBandwidth
** throughput   datasize / (RTT - Delay)
** serversize      [b]     Messagesize from server
** clientsize      [b]     Messagesize from client

Outlined benchmarks for optimizing WebSockets:

- Clientlatency() @ Bandwidth ≈ 0, Clientmessage = 1b
    Requires some tweaking of BenchMarkTools.Trial, postponed.
- Serverlatency() @ Bandwidth ≈ 0, Servermessage = 1b
    Requires some tweaking of BenchMarkTools.Trial, postponed.
- Maximized bandwidth(VSIZE) @ Up2Down ≈ 1
    Not using BenchmarkTools directly
- ClientRTT(VDELAY) @ Bandwidth ≈ 0, Clientmessage = 1b
- ServerRTT(VDELAY) @ Bandwidth ≈ 0, Servermessage = 1b
    The last two checks vs. possibly related issue #1608 on ZMQ

Outlined benchmarks for developing an application using WebSockets (postponed):
- ClientRTT(VSIZE) @ Bandwidth[0%, 50%, 100%]
    ClientRTT determines user responsiveness, e.g. mouse clicks in a browser
    This requires small messages or threaded read / write (async has to wait for
    the server finishing sending its message)
- Bandwidth(VUP2DOWN, VSIZE)
    This depends on the network in use and operating system
        [(up2down, msize)  for up2down in VUP2DOWN for msize in VSIZE]
- Serverlatency(VSIZE) @ Bandwidth[0%, 50%, 100%]
    The outliers (due to e.g. semi-random garbage collection) determines
    choice of message size and buffers for media streams
- Clientlatency() @ 100% bandwidth @ Up2Down
    The outliers (due to e.g. semi-random garbage collection) determines
    choice of message size and buffers for long-running calculations which
    connects to a server for distribution of results
=#
if !@isdefined(SRCPATH)
    const SRCPATH = Base.source_dir() == nothing ? Pkg.dir("WebSockets", "benchmark") : Base.source_dir()
    const LOGGINGPATH = realpath(joinpath(SRCPATH, "../logutils/"))
    SRCPATH ∉ LOAD_PATH && push!(LOAD_PATH, SRCPATH)
    LOGGINGPATH ∉ LOAD_PATH && push!(LOAD_PATH, LOGGINGPATH)
    include(joinpath(SRCPATH, "functions_open_browsers.jl"))
end
# Don't run this if the output files are recent?
include("benchmark_prepare.jl")
# include detailed benchmarks using the results from the above as nominal.


text:
# Practical Machine learning
# Decision Trees based learning - Random Forests example
# Chapter 5

using MLBase

# internal data structure for a random forest classifier
type Classifier
    n_samples::Int
    n_features::Int
    n_max_features::Int
    label_mapping::LabelMap
    outtype::DataType
    improvements::Vector{Float64}
    oob_error::Float64
    trees::Vector{Tree}

    function Classifier(rf, x, y)
        n_samples, n_features = size(x)

        if n_samples != length(y)
            throw(DimensionMismatch(""))
        end

        n_max_features = resolve_max_features(rf.max_features, n_features)
        @assert 0 < n_max_features <= n_features

        label_mapping = labelmap(y)
        outtype = eltype(y)
        improvements = zeros(Float64, n_features)
        trees = Array(Tree, rf.n_estimators)
        new(n_samples, n_features, n_max_features, label_mapping, outtype, improvements, NaN, trees)
    end
end

typealias RandomForestClassifier RandomForest{Classifier}

function RandomForestClassifier(;n_estimators::Int=10, max_features::Union{Integer,AbstractFloat,Symbol}=:sqrt, max_depth=nothing, min_samples_split::Int=2, criterion::Symbol=:gini)
    if !(is(criterion, :gini) || is(criterion, :entropy))
        error("criterion is invalid (got: $criterion)")
    end
    RandomForest{Classifier}(n_estimators, max_features, max_depth, min_samples_split, criterion)
end

function set_weight!(bootstrap::Vector{Int}, sample_weight::Vector{Float64})
    @assert length(bootstrap) == length(sample_weight)
    # initialize weight
    for i in 1:endof(sample_weight)
        sample_weight[i] = 0.0
    end
    # set weight
    for i in bootstrap
        sample_weight[i] += 1.0
    end
end

function fit{T<:TabularData}(rf::RandomForestClassifier, x::T, y::AbstractVector)
    learner = Classifier(rf, x, y)
    y_encoded = labelencode(learner.label_mapping, y)
    n_samples = learner.n_samples

    # pre-allocation
    bootstrap = Array(Int, n_samples)
    sample_weight = Array(Float64, n_samples)
    oob_error = 0.

    for b in 1:rf.n_estimators
        rand!(bootstrap, 1:n_samples)
        set_weight!(bootstrap, sample_weight)
        example = Trees.Example{T}(x, y_encoded, sample_weight)
        tree = Trees.Tree()
        Trees.fit(tree, example, rf.criterion, learner.n_max_features, rf.max_depth, rf.min_samples_split)
        learner.trees[b] = tree

        hit = 0
        miss = 0
        for s in 1:n_samples
            if sample_weight[s] != 0.0
                continue
            end

            # out-of-bag sample
            if Trees.predict(tree, vec(x[s, :])) == y_encoded[s]
                hit += 1
            else
                miss += 1
            end
        end
        oob_error += miss / (hit + miss)
    end

    set_improvements!(learner)
    learner.oob_error = oob_error / rf.n_estimators
    rf.learner = learner
    return
end

function predict{T<:TabularData}(rf::RandomForestClassifier, x::T)
    if is(rf.learner, nothing)
        error("not yet trained")
    end

    n_samples = size(x, 1)
    output = Array(Int, n_samples)
    n_labels = length(rf.learner.label_mapping)

    for i in 1:n_samples
        counts = zeros(Int, n_labels)
        for b in 1:rf.n_estimators
            tree = rf.learner.trees[b]
            vote = Trees.predict(tree, vec(x[i, :]))
            counts[vote] += 1
        end
        output[i] = indmax(counts)
    end

    # TODO: use labeldecode method when MLBase.jl is released
    # labeldecode(rf.learner.label_mapping, output)
    predicted = Array(rf.learner.outtype, n_samples)
    for (i, o) in enumerate(output)
        predicted[i] = rf.learner.label_mapping.vs[o]
    end
    predicted
end

function normalize!(v)
    s = sum(v)
    for i in 1:endof(v)
        v[i] = v[i] / s
    end
    return
end

function oob_error(rf::RandomForestClassifier)
    if is(rf.learner, nothing)
        error("not yet trained")
    end

    rf.learner.oob_error
end


text:
import REPL
using REPL.TerminalMenus
using Test

TerminalMenus.config(supress_output=true)

function simulateInput(expectedResult, menu::TerminalMenus.AbstractMenu, keys...)
    # If we cannot write to the buffer, skip the test
    !(:buffer in fieldnames(typeof(stdin))) && return true

    keydict =  Dict(:up => "\e[A",
                    :down => "\e[B",
                    :enter => "\r")

    for key in keys
        if isa(key, Symbol)
            write(stdin.buffer, keydict[key])
        else
            write(stdin.buffer, "$key")
        end
    end

    request(menu) == expectedResult
end

include("radio_menu.jl")
include("multiselect_menu.jl")

# Other test

# scroll must only accept symbols
@test_throws TypeError TerminalMenus.config(scroll=true)
# :foo is not a valid scroll option
@test_throws ArgumentError TerminalMenus.config(scroll=:foo)
# Test scroll wrap
TerminalMenus.config(scroll=:wrap)
@test TerminalMenus.CONFIG[:scroll_wrap] == true
# Updating some params shouldn't change other ones
TerminalMenus.config(charset=:ascii)
@test TerminalMenus.CONFIG[:scroll_wrap] == true
TerminalMenus.config(scroll=:nowrap)
@test TerminalMenus.CONFIG[:scroll_wrap] == false

max_stars_repo_name: djsegal/julia-fork

max_issues_repo_name: djsegal/julia-fork

max_forks_repo_name: djsegal/julia-fork


text:
function filter_CIR(δ, γ, σ, λ, data; silence = false, trim0 = false)

    times = keys(data) |> collect |> sort
    Λ_of_t = Dict{Float64, Array{Int64,1}}()
    wms_of_t = Dict{Float64, Array{Float64,1}}()
    θ_of_t = Dict{Float64, Float64}()

    filtered_θ, filtered_Λ, filtered_wms = update_CIR_params([1.], δ, γ/σ^2, λ, [0], data[0])

    Λ_of_t[0] = filtered_Λ
    wms_of_t[0] = filtered_wms # = 1.
    θ_of_t[0] = filtered_θ

    for k in 1:(length(times)-1)
        if (!silence)
            println("Step index: $k")
            println("Number of components: $(length(filtered_Λ))")
        end
        filtered_θ, filtered_Λ, filtered_wms = get_next_filtering_distribution(filtered_Λ, filtered_wms, filtered_θ, times[k], times[k+1], δ, γ, σ, λ, data[times[k+1]])
        if trim0
            filtered_Λ, filtered_wms = filtered_Λ[filtered_wms.>0], filtered_wms[filtered_wms.>0]
        end
        Λ_of_t[times[k+1]] = filtered_Λ
        wms_of_t[times[k+1]] = filtered_wms
        θ_of_t[times[k+1]] = filtered_θ
    end

    return Λ_of_t, wms_of_t, θ_of_t

end

function update_CIR_params(wms::Array{Ty,1}, δ::Real, θ::Real, λ::Real, Λ, y::Array{Tz,1}) where {Ty<:Number,Tz<:Integer}
    α = δ/2#Alternative parametrisation

    ny = sum(y)
    J = length(y)

    return θ + J*λ, Λ .+ ny, next_wms_from_wms_prime(wms, Λ, y, θ, α)
end

function get_next_filtering_distribution(current_Λ, current_wms, current_θ, current_time, next_time, δ, γ, σ, λ, next_y)
    predicted_θ, predicted_Λ, predicted_wms = predict_CIR_params(current_wms, δ, current_θ, γ, σ, current_Λ, next_time-current_time)
    filtered_θ, filtered_Λ, filtered_wms = update_CIR_params(predicted_wms, δ, predicted_θ, λ, predicted_Λ, next_y)

    return filtered_θ, filtered_Λ, filtered_wms
end

function predict_CIR_params(wms::Array{Ty,1}, δ::Ty, θ::Ty, γ::Ty, σ::Ty, Λ, t::Ty) where Ty<:Number

    p = γ/σ^2*1/(θ*exp(2*γ*t) + γ/σ^2 - θ)
    θ_new = p * θ * exp(2*γ*t)
    Λ_new = 0:maximum(Λ)

    return θ_new, Λ_new, next_wms_prime_from_wms(wms, Λ, t, θ, γ, σ)

end

max_stars_repo_name: konkam/DualOptimalFiltering

max_issues_repo_name: konkam/DualOptimalFiltering

max_forks_repo_name: konkam/DualOptimalFiltering


text:
using Test
using Forecast
import DataFrames: DataFrame

@testset "datasets: co2" begin

    co2_paper = co2()
    co2_df = co2(true)
    @test co2_paper isa DataFrame
    @test co2_df isa DataFrame
    @test size(co2_paper) == (4609,2)
    @test size(co2_df) == (17166, 17)

end

@testset "datasets: seaborne" begin

    sb = seaborne()
    sb_full = seaborne(true)
    @test sb isa DataFrame
    @test sb_full isa DataFrame
    @test size(sb) == (2199,4)
    @test size(sb_full) == (13209, 10)

end

@testset "datasets: air" begin

    a = air()
    @test a isa DataFrame
    @test size(a) == (144,2)

end

@testset "datasets: quakes" begin

    q = quakes()
    @test q isa DataFrame
    @test size(q) == (71,2)

end


text:
function initialize!(integrator, cache::ROCK2ConstantCache)
  integrator.kshortsize = 2
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
  integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal

  # Avoid undefined entries if k is an array of arrays
  integrator.fsallast = zero(integrator.fsalfirst)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
end

@muladd function perform_step!(integrator, cache::ROCK2ConstantCache, repeat_step=false)
  @unpack t, dt, uprev, u, f, p, fsalfirst = integrator
  @unpack ms, fp1, fp2, recf = cache
  maxeig!(integrator, cache)
  # The the number of degree for Chebyshev polynomial
  mdeg = Int(floor(sqrt((1.5 + dt * integrator.eigen_est)/0.811) + 1))
  if mdeg >= 200
    mdeg = 200
  end
  cache.mdeg = max(mdeg, 3) - 2
  cache.mdeg != cache.mdegprev && choosedeg!(cache)
  # recurrence
  # for the first stage
  temp1 = dt * recf[cache.recind][1]
  ci1 = t + temp1
  ci2 = t + temp1
  ci3 = t
  gprev2 = copy(uprev)
  gprev = uprev + temp1 * fsalfirst
  ms[cache.mdeg] < 2 && ( u = gprev )
  # for the second to the ms[cache.mdeg] th stages
  for i in 2:ms[cache.mdeg]
    μ, κ = recf[cache.recind + (i - 2)]
    ν = -1 - κ
    dtμ = dt*μ
    ci1 = dtμ - ν * ci2 - κ * ci3
    u = dtμ * u - ν * gprev - κ * gprev2
    i < ms[cache.mdeg] && (gprev2 = gprev; gprev = u)
    ci3 = ci2
    ci2 = ci1
  end # end if
  # two-stage finishing procedure.
  temp1 = dt * fp1[cache.mdeg]
  temp2 = dt * fp2[cache.mdeg]
  gprev2 = f(u, p, ci1)
  gprev = u + temp1 * gprev2
  ci1 += temp1
  u = f(gprev, p, ci1)
  temp3 = temp2 * (u - gprev2)
  u = gprev + temp1 * u + temp3
  # error estimate
  if integrator.opts.adaptive
    atmp = calculate_residuals(temp3, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm)
    integrator.EEst = integrator.opts.internalnorm(atmp)
  end
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast = f(u, p, t+dt)
  integrator.u = u
end

function initialize!(integrator, cache::ROCK2Cache)
  integrator.kshortsize = 2
  resize!(integrator.k, integrator.kshortsize)
  integrator.fsalfirst = cache.fsalfirst  # done by pointers, no copying
  integrator.fsallast = cache.k
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.f(integrator.fsalfirst, integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
end

@muladd function perform_step!(integrator, cache::ROCK2Cache, repeat_step=false)
  @unpack t, dt, uprev, u, f, p, fsalfirst = integrator
  @unpack k, k2, tmp, gprev2, gprev, atmp = cache
  @unpack ms, fp1, fp2, recf = cache.constantcache
  ccache = cache.constantcache
  maxeig!(integrator, cache)
  # The the number of degree for Chebyshev polynomial
  mdeg = Int(floor(sqrt((1.5 + dt * integrator.eigen_est)/0.811) + 1))
  if mdeg >= 200
    mdeg = 200
  end
  ccache.mdeg = max(mdeg, 3) - 2
  ccache.mdeg != ccache.mdegprev && choosedeg!(cache)
  # recurrence
  # for the first stage
  temp1 = dt * recf[ccache.recind][1]
  ci1 = t + temp1
  ci2 = t + temp1
  ci3 = t
  @. gprev2 = uprev
  @. gprev = uprev + temp1 * fsalfirst
  ms[ccache.mdeg] < 2 && ( @. u = gprev )
  # for the second to the ms[ccache.mdeg] th stages
  for i in 2:ms[ccache.mdeg]
    μ, κ = recf[cache.recind + (i - 2)]
    ν = κ - 1
    temp1 = dt * μ
    temp2 = 1 + κ
    temp3 = -κ
    ci1 = temp1 + temp2 * ci2 + temp3 * ci3
    @. u = temp1 * u + temp2 * gprev + temp3 * gprev2
    i < ms[ccache.mdeg] && (gprev2 .= gprev; gprev .= u)
    ci3 = ci2
    ci2 = ci1
  end # end if
  # two-stage finishing procedure.
  temp1 = dt * fp1[ccache.mdeg]
  temp2 = dt * fp2[ccache.mdeg]
  f(k, u, p, ci1)
  @. gprev = u + temp1 * k
  ci1 += temp1
  f(k2, gprev, p, ci1)
  @. tmp = temp2 * (k2 - k)
  @. u = gprev + temp1 * k2 + tmp
  # error estimate
  if integrator.opts.adaptive
    calculate_residuals!(atmp, tmp, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm)
    integrator.EEst = integrator.opts.internalnorm(atmp)
  end
  integrator.k[1] = integrator.fsalfirst
  f(integrator.fsallast, u, p, t+dt)
  integrator.k[2] = integrator.fsallast
  integrator.u = u
end

function initialize!(integrator, cache::ROCK4ConstantCache)
  integrator.kshortsize = 2
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
  integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
  # Avoid undefined entries if k is an array of arrays
  integrator.fsallast = zero(integrator.fsalfirst)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
end

@muladd function perform_step!(integrator, cache::ROCK4ConstantCache, repeat_step=false)
  @unpack t, dt, uprev, u, f, p, fsalfirst = integrator
  @unpack ms, fpa, fpb, fpbe, recf = cache
  maxeig!(integrator, cache)
  # The the number of degree for Chebyshev polynomial
  mdeg = Int(floor(sqrt((3.0 + dt * integrator.eigen_est)/0.353) + 1))
  if mdeg >= 152
    mdeg = 152
  end
  cache.mdeg = max(mdeg, 5) - 4
  cache.mdeg != cache.mdegprev && choosedeg!(cache)
  # recurrence
  # for the first stage
  temp1 = dt * recf[cache.recind][1]
  ci1 = t + temp1
  ci2 = t + temp1
  ci3 = t
  gprev2 = copy(uprev)
  gprev = uprev + temp1 * fsalfirst
  ms[cache.mdeg] < 2 && ( u = gprev )
  # for the second to the ms[cache.mdeg] th stages
  for i in 2:ms[cache.mdeg]
    μ, κ = recf[cache.recind + (i - 2)]
    ν = -1 - κ
    dtμ = dt*μ
    ci1 = dtμ - ν * ci2 - κ * ci3
    u = dtμ * u - ν * gprev - κ * gprev2
    i < ms[cache.mdeg] && (gprev2 = gprev; gprev = u)
    ci3 = ci2
    ci2 = ci1
  end
  # 4-stage finishing procedure.
  # Stage-1
  temp1 = dt * fpa[cache.mdeg][1]
  gprev = f(u, p, ci1)
  gprev3 = u + temp1 * gprev
  # Stage-2
  ci2 = ci1 + temp1
  temp1 = dt * fpa[cache.mdeg][2];
  temp2 = dt * fpa[cache.mdeg][3];
  gprev2 = f(gprev3, p, ci1)
  gprev4 = u + temp1 * gprev + temp2 * gprev2
  # Stage-3
  ci2 = ci1 + temp1 +temp2
  temp1 = dt * fpa[cache.mdeg][4]
  temp2 = dt * fpa[cache.mdeg][5]
  temp3 = dt * fpa[cache.mdeg][6]
  gprev3 = f(gprev4, p, ci2)
  gprev5 = u + temp1 * gprev + temp2 * gprev2 + temp3 * gprev3
  #Stage-4
  ci2 = ci1 + temp1 + temp2 + temp3
  temp1 = dt * fpb[cache.mdeg][1]
  temp2 = dt * fpb[cache.mdeg][2]
  temp3 = dt * fpb[cache.mdeg][3]
  temp4 = dt * fpb[cache.mdeg][4]
  gprev4 = f(gprev5, p, ci2)
  u = u + temp1 * gprev + temp2 * gprev2 + temp3 * gprev3 + temp4 * gprev4
  #Error estimate (embedded method of order 3)
  temp1 = dt * fpbe[cache.mdeg][1] - temp1
  temp2 = dt * fpbe[cache.mdeg][2] - temp2
  temp3 = dt * fpbe[cache.mdeg][3] - temp3
  temp4 = dt * fpbe[cache.mdeg][4] - temp4
  temp5 = dt * fpbe[cache.mdeg][5]
  gprev5 = f(u, p, t + dt)
  temp5 = temp1 * gprev + temp2 * gprev2 + temp3 * gprev3 + temp4 * gprev4 + temp5 * gprev5
  if integrator.opts.adaptive
    atmp = calculate_residuals(temp5, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm)
    integrator.EEst = integrator.opts.internalnorm(atmp)
  end
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast = f(u, p, t+dt)
  integrator.u = u
end

function initialize!(integrator, cache::ROCK4Cache)
  integrator.kshortsize = 2
  resize!(integrator.k, integrator.kshortsize)
  integrator.fsalfirst = cache.fsalfirst
  integrator.fsallast = cache.k
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.f(integrator.fsalfirst, integrator.uprev, integrator.p, integrator.t)
end

@muladd function perform_step!(integrator, cache::ROCK4Cache, repeat_step=false)
  @unpack t, dt, uprev, u, f, p, fsalfirst = integrator
  @unpack gprev, gprev2, gprev3, gprev4, gprev5, tmp, atmp, k, k2, k3, k4, k5 = cache
  @unpack ms, fpa, fpb, fpbe, recf = cache.constantcache
  ccache = cache.constantcache
  maxeig!(integrator, cache)
  # The the number of degree for Chebyshev polynomial
  mdeg = Int(floor(sqrt((3.0 + dt * integrator.eigen_est)/0.353) + 1))
  if mdeg >= 152
    mdeg = 152
  end
  ccache.mdeg = max(mdeg, 5) - 4
  ccache.mdeg != ccache.mdegprev && choosedeg!(cache)
  # recurrence
  # for the first stage
  temp1 = dt * recf[ccache.recind][1]
  ci1 = t + temp1
  ci2 = t + temp1
  ci3 = t
  @. gprev2 = uprev
  @. gprev = uprev + temp1 * fsalfirst
  ms[ccache.mdeg] < 2 && ( @. u = gprev )
  # for the second to the ms[ccache.mdeg] th stages
  for i in 2:ms[ccache.mdeg]
    μ, κ = recf[cache.recind + (i - 2)]
    ν = κ - 1
    temp1 = dt * μ
    temp2 = 1 + κ
    temp3 = -κ
    ci1 = temp1 + temp2 * ci2 + temp3 * ci3
    @. u = temp1 * u + temp2 * gprev + temp3 * gprev2
    i < ms[ccache.mdeg] && (gprev2 .= gprev; gprev .= u)
    ci3 = ci2
    ci2 = ci1
  end
  # 4-stage finishing procedure.
  # Stage-1
  temp1 = dt * fpa[ccache.mdeg][1]
  f(k, u, p, ci1)
  @. gprev3 = u + temp1 * k
  # Stage-2
  ci2 = ci1 + temp1
  temp1 = dt * fpa[ccache.mdeg][2];
  temp2 = dt * fpa[ccache.mdeg][3];
  f(k2, gprev3, p, ci1)
  @. gprev4 = u + temp1 * k + temp2 * k2
  # Stage-3
  ci2 = ci1 + temp1 +temp2
  temp1 = dt * fpa[ccache.mdeg][4]
  temp2 = dt * fpa[ccache.mdeg][5]
  temp3 = dt * fpa[ccache.mdeg][6]
  f(k3, gprev4, p, ci2)
  @. gprev5 = u + temp1 * k + temp2 * k2 + temp3 * k3
  #Stage-4
  ci2 = ci1 + temp1 + temp2 + temp3
  temp1 = dt * fpb[ccache.mdeg][1]
  temp2 = dt * fpb[ccache.mdeg][2]
  temp3 = dt * fpb[ccache.mdeg][3]
  temp4 = dt * fpb[ccache.mdeg][4]
  f(k4, gprev5, p, ci2)
  @. u = u + temp1 * k + temp2 * k2 + temp3 * k3 + temp4 * k4
  #Error estimate (embedded method of order 3)
  temp1 = dt * fpbe[ccache.mdeg][1] - temp1
  temp2 = dt * fpbe[ccache.mdeg][2] - temp2
  temp3 = dt * fpbe[ccache.mdeg][3] - temp3
  temp4 = dt * fpbe[ccache.mdeg][4] - temp4
  temp5 = dt * fpbe[ccache.mdeg][5]
  f(k5, u, p, t + dt)
  @. tmp = temp1 * k + temp2 * k2 + temp3 * k3 + temp4 * k4 + temp5 * k5
  if integrator.opts.adaptive
    calculate_residuals!(atmp, tmp, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm)
    integrator.EEst = integrator.opts.internalnorm(atmp)
  end
  integrator.k[1] = integrator.fsalfirst
  f(integrator.fsallast, u, p, t+dt)
  integrator.k[2] = integrator.fsallast
  integrator.u = u
end


text:
#=
Implementation of Genetic Algorithm

The constructor takes following keyword arguments:

- `populationSize`: The size of the population
- `crossoverRate`: The fraction of the population at the next generation, not including elite children, that is created by the crossover function.
- `mutationRate`: Probability of chromosome to be mutated
- `ɛ`/`epsilon`: Positive integer specifies how many individuals in the current generation are guaranteed to survive to the next generation. Floating number specifies fraction of population.
- `selection`: [Selection](@ref) function
- `crossover`: [Crossover](@ref) function (default: `identity`)
- `mutation`: [Mutation](@ref) function (default: `identity`)
- `metrics` is a collection of convergence metrics.
=#
struct GA <: AbstractOptimizer
    populationSize::Int
    crossoverRate::Float64
    mutationRate::Float64
    ɛ::Real
    selection
    crossover
    mutation
    metrics::ConvergenceMetrics

    GA(;
        populationSize::Int = 50,
        crossoverRate::Float64 = 0.8,
        mutationRate::Float64 = 0.1,
        ɛ::Real = 0,
        epsilon::Real = ɛ,
        selection = ((x, n) -> 1:n),
        crossover = identity,
        mutation = identity,
        metrics = ConvergenceMetric[AbsDiff(1e-12)],
    ) = new(
        populationSize,
        crossoverRate,
        mutationRate,
        epsilon,
        selection,
        crossover,
        mutation,
        metrics,
    )
end

population_size(method::GA) = method.populationSize
default_options(method::GA) = (iterations = 1000,)
summary(m::GA) =
    "GA[P=$(m.populationSize),x=$(m.crossoverRate),μ=$(m.mutationRate),ɛ=$(m.ɛ)]"
show(io::IO, m::GA) = print(io, summary(m))

mutable struct GAState <: AbstractOptimizerState
    iteration
    start_time
    stop_time
    metrics::ConvergenceMetrics # collection of convergence metrics

    eliteSize::Int
    population
    fittest
end

value(s::GAState) = objectives(s.fittest)
minimizer(s::GAState) = variables(s.fittest)
copy(s::GAState) = GAState(s.iteration, s.start_time, s.stop_time, copy(s.metrics), s.eliteSize, copy(s.population), copy(s.fittest))

"""Initialization of GA algorithm state"""
function initial_state(method::GA, objective, population, options)
    eliteSize = isa(method.ɛ, Int) ? method.ɛ : round(Int, method.ɛ * method.populationSize)

    fitness = value!(objective, population)
    _, idx = findmin(fitness)

    return GAState(0, 0, 0, copy(method.metrics), eliteSize, population, copy(population[idx]))
end

function update_state!(
    method::GA,
    state,
    objective,
    constraints,
    options
)
    parents = state.population

    populationSize = method.populationSize
    rng = options.rng
    offspring = copy(parents)

    # select offspring
    selected = method.selection(state.population, populationSize, rng = rng)

    # perform mating
    offspringSize = populationSize - state.eliteSize
    recombine!(offspring, parents, selected, method, offspringSize, rng = rng)

    # Elitism (copy population individuals before they pass to the offspring & get mutated)
    idxs = sortperm(vec(objectives(state.population)))
    for i = 1:state.eliteSize
        subs = offspringSize + i
        offspring[subs] = copy(parents[idxs[i]])
    end

    # perform mutation
    mutate!(offspring, method, constraints, rng = rng)

    # calculate fitness of the population
    evaluate!(objective, offspring, constraints)

    # select the best individual
    _, idx = findmin(objectives(offspring))
    state.fittest = offspring[idx]

    # replace population
    parents .= offspring

    return false
end

function recombine!(
    offspring,
    parents,
    selected,
    method,
    n = length(selected);
    rng::AbstractRNG = Random.default_rng(),
)
    mates = ((i, i == n ? i - 1 : i + 1) for i = 1:2:n)
    for (k, j) in mates
        pk, pj = selected[k], selected[j]
        p1, p2 = parents[pk], parents[pj]
        if rand(rng) < method.crossoverRate
            offspring[pk], offspring[pj] = method.crossover(p1, p2, rng = rng)
        else
            offspring[pk], offspring[pj] = p1, p2
        end
    end

end

function mutate!(population, method, constraints; rng::AbstractRNG = Random.default_rng())
    n = length(population)
    for i = 1:n
        if rand(rng) < method.mutationRate
            method.mutation(population[i], rng = rng)
        end

        (constraints, population[i])
    end
end

function evaluate!(objective, population, constraints)
    # calculate fitness of the population
    value!(objective, population)
    # apply penalty to fitness
    # penalty!(fitness, constraints, population)
end


text:
using InteractiveUtils

export allsubtypes, showsubtypes, pedigree

function allsubtypes(t::Type, result=Vector{Type}())
    push!(result, t)
    for st in subtypes(t)
	allsubtypes(st, result)
    end
    return result
end

function showsubtypes(t::Type, level=0)
    indent1 = "  "
    println("$(repeat(indent1, level))$t")
    for st in subtypes(t)
        showsubtypes(st, level + 1)
    end
end

function pedigree(t::Type, result=Vector{Type}())
    push!(result, t)
    if supertype(t) != t
        pedigree(supertype(t), result)
    end
    result
end