text:
@testset "input" begin
    a = InteractBase.input()
    @test widgettype(a) == :input
    @test observe(a)[] == ""
    a = InteractBase.input(typ = "number");
    @test observe(a)[] == 0
    s = Observable{Any}(12)
    a = InteractBase.input(s, typ = "number");
    @test observe(a)[] == s[]

    @test widgettype(input(Bool)) == :toggle
    @test widgettype(input(Dates.Date)) == :datepicker
    @test widgettype(input(Dates.Time)) == :timepicker
    @test widgettype(input(Color)) == :colorpicker
    @test widgettype(input(String)) == :textbox
    @test widgettype(input(Int)) == :spinbox

    @test widgettype(widget(true)) == :toggle
    @test widgettype(widget(Dates.Date(Dates.now()))) == :datepicker
    @test widgettype(widget(Dates.Time(Dates.now()))) == :timepicker
    @test widgettype(widget(colorant"red")) == :colorpicker
    @test widgettype(widget("")) == :textbox
    @test widgettype(widget(1)) == :spinbox
    @test widgettype(widget(1:100)) == :slider
    @test widgettype(widget(["a", "b", "c"])) == :togglebuttons
end

@testset "input widgets" begin
    a = filepicker()
    @test widgettype(a) == :filepicker
    @test a["filename"][] == ""
    @test a["path"][] == ""
    a["path"][] = "/home/Jack/documents/test.csv"
    @test a["path"][] == observe(a)[] == "/home/Jack/documents/test.csv"

    a = datepicker(value = Dates.Date(01,01,01))
    b = datepicker(Dates.Date(01,01,01))
    @test observe(a)[] == observe(b)[] == Dates.Date(01,01,01)
    @test widgettype(a) == :datepicker

    a = colorpicker(value = colorant"red")
    b = colorpicker(colorant"red")
    @test observe(a)[] == observe(b)[] == colorant"red"
    @test widgettype(a) == :colorpicker

    a = spinbox(label = "")
    @test widgettype(a) == :spinbox
    @test observe(a)[] == nothing

    a = textbox();
    @test widgettype(a) == :textbox

    @test observe(a)[] == ""
    s = "asd"
    a = textbox(value = s);
    @test observe(a)[] == "asd"

    a = textarea(label = "test");
    @test widgettype(a) == :textarea

    @test observe(a)[] == ""
    s = "asd"
    a = textarea(value = s);
    @test observe(a)[] == "asd"

    a = autocomplete(["aa", "bb", "cc"], value = "a");
    @test widgettype(a) == :autocomplete

    @test observe(a)[] == "a"

    a = button("Press me!", value = 12)
    @test widgettype(a) == :button
    @test observe(a)[] == 12

    a = toggle(label = "Agreed")
    @test widgettype(a) == :toggle

    @test observe(a)[] == false
    s = Observable(true)
    a = toggle(s, label = "Agreed")
    @test observe(a)[] == true

    a = togglecontent(checkbox("Yes, I am sure"), "Are you sure?")
    @test widgettype(a) == :togglecontent

    @test observe(a)[] == false
    s = Observable(true)
    a = togglecontent(checkbox("Yes, I am sure"), "Are you sure?", value = s)
    @test observe(a)[] == true

    v = slider([0, 12, 22], value = 12)
    @test widgettype(v) == :slider

    @test observe(v)[] == 12
    @test v["index"][] == 2
    # v["internalvalue"][] = 3
    # @test observe(v)[] == 22
end

@testset "slider" begin
    @test isfile(InteractBase.nouislider_min_js)
    @test isfile(InteractBase.nouislider_min_css)
    w = Dates.Date("2000-11-11") : Day(1) : Dates.Date("2000-12-12")
    s = InteractBase.rangeslider(w, value = [w[10], w[20]])
    @test observe(s)[] == [w[10], w[20]]
    @test observe(s["index"])[] == [10, 20]
    observe(s["index"])[] = [13, 14]
    sleep(0.1)
    @test observe(s)[] == [w[13], w[14]]

    w = 1:5:500
    s = InteractBase.rangepicker(w, value = [w[10], w[20]])
    @test collect(keys(components(s))) == [:input1, :input2, :slider, :changes]
    @test observe(s)[] == [w[10], w[20]] == observe(s["slider"])[]
    @test observe(s["slider"])[] == [w[10], w[20]]
    observe(s["slider"])[] = [w[13], w[14]]
    sleep(0.1)
    @test observe(s)[] == [w[13], w[14]]
    @test observe(s["input1"])[] == w[13]
    @test observe(s["input2"])[] == w[14]
end

@testset "options" begin
    a = dropdown(["a", "b", "c"])
    @test widgettype(a) == :dropdown

    @test observe(a)[] == "a"
    a = dropdown(OrderedDict("a" => 1, "b" => 2, "c" => 3))
    @test observe(a)[] == 1
    a = dropdown(OrderedDict("a" => 1, "b" => 3, "c" => 4), value = 3)
    @test observe(a)[] == 3
    @test observe(a, "index")[] == 2

    v = [0.1, 0.2, 1.2]
    a = dropdown(OrderedDict("a" => 1, "b" => 3, "c" => v), value = v)
    @test observe(a)[] == v
    @test observe(a, "index")[] == 3

    v = [0.1, 0.2, 1.2]
    a = dropdown(OrderedDict("a" => 1, "b" => 3, "c" => v), value = [3, v], multiple = true)
    @test observe(a)[] == [3, v]
    @test observe(a, "index")[] == [2, 3]

    a = togglebuttons(["a", "b", "c"])
    @test widgettype(a) == :togglebuttons

    @test observe(a)[] == "a"
    a = togglebuttons(OrderedDict("a" => 1, "b" => 2, "c"=>3))
    @test observe(a)[] == 1
    a = togglebuttons(OrderedDict("a" => 1, "b" => 2, "c" => 4), value = 4)
    @test observe(a)[] == 4

    a = radiobuttons(["a", "b", "c"])
    @test widgettype(a) == :radiobuttons

    @test observe(a)[] == "a"
    a = radiobuttons(OrderedDict("a" => 1, "b" => 2, "c" => 3))
    @test observe(a)[] == 1
    a = radiobuttons(OrderedDict("a" => 1, "b" => 2, "c" => 3), value = 3, label = "Test")
    @test observe(a)[] == 3

end

@testset "ijulia" begin
    @test !InteractBase.isijulia()
end

@testset "widget" begin
    s = slider(1:100, value = 12)
    w = InteractBase.Widget{:test}(components(s), scope = Widgets.scope(s), output = Observable(1))
    @test observe(w)[] == 1
    @test widgettype(s) == :slider
    @test widgettype(w) == :test
    @test w["index"][] == 12
    w = Widget(w, output = scope(s)["index"])
    @test observe(w)[] == 12

    w = InteractBase.widget(Observable(1))
    @test w isa Observable
end

@testset "manipulate" begin
    ui = @manipulate for r = 0:.05:1, g = 0:.05:1, b = 0:.05:1
        RGB(r,g,b)
    end
    @test observe(ui)[] == RGB(0.5, 0.5, 0.5)
    observe(ui, :r)[] = 0.1
    sleep(0.1)
    @test observe(ui)[] == RGB(0.1, 0.5, 0.5)

    ui = @manipulate throttle = 1 for r = 0:.05:1, g = 0:.05:1, b = 0:.05:1
        RGB(r,g,b)
    end
    observe(ui, :r)[] = 0.1
    sleep(0.1)
    observe(ui, :r)[] = 0.3
    sleep(0.1)
    observe(ui, :g)[] = 0.1
    sleep(0.1)
    observe(ui, :g)[] = 0.3
    sleep(0.1)
    observe(ui, :b)[] = 0.1
    sleep(0.1)
    observe(ui, :b)[] = 0.3
    sleep(0.1)
    @test observe(ui)[] != RGB(0.3, 0.3, 0.3)
    sleep(1.5)
    @test observe(ui)[] == RGB(0.3, 0.3, 0.3)
end

@testset "output" begin
    @test isfile(InteractBase.katex_min_js)
    @test isfile(InteractBase.katex_min_css)
    l = Observable("\\sum_{i=1}^{\\infty} e^i")
    a = latex(l)
    @test widgettype(a) == :latex
    @test observe(a)[] == l[]
    l[] == "\\sum_{i=1}^{12} e^i"
    @test observe(a)[] == l[]

    @test isfile(joinpath(InteractBase.external_assets_path, "prism.js"))
    @test isfile(joinpath(InteractBase.external_assets_path, "prism.css"))

    l = Observable("1+1+exp(2)")
    a = highlight(l)
    @test widgettype(a) == :highlight
    @test observe(a)[] == l[]
    l[] == "1-1"
    @test observe(a)[] == l[]

    l = Observable("1+1+exp(2)")
    a = widget(Val(:highlight), l)
    @test widgettype(a) == :highlight
    @test observe(a)[] == l[]
    l[] == "1-1"
    @test observe(a)[] == l[]

    a = alert()
    a("Error!")
    @test a["text"] isa Observable
    @test a["text"][] == "Error!"

    a = widget(Val(:alert), "Error 2!")
    a()
    @test a["text"] isa Observable
    @test a["text"][] == "Error 2!"

    a = confirm()
    a("Error!")
    @test a["text"] isa Observable
    @test a["text"][] == "Error!"
    @test observe(a)[] == false
    @test a["function"](1) === nothing

    a = widget(Val(:confirm), "Error 2!")
    a()
    @test a["text"] isa Observable
    @test a["text"][] == "Error 2!"
    @test observe(a)[] == false

    v = Any["A"]
    f = notifications(v)
    sleep(0.1)
    @test observe(f)[] == v
    list = children(f.scope.dom[])[1]

    @test begin
        Widgets.scope(f)[:to_delete][] = 1
        sleep(0.1)
        observe(f)[] == []
    end

    v = OrderedDict("a" => checkbox(), "b" => 12)
    wdg = InteractBase.accordion(v, multiple = true)
    sleep(0.1)
    @test observe(wdg)[] == Int[]
    @test observe(wdg["options"])[] == v
    observe(wdg)[] = [1]
    sleep(0.1)
    @test observe(wdg)[] == [1]
    observe(wdg["options"])[] = OrderedDict("a" => 12)
    sleep(0.1)
    @test observe(wdg)[] == [1]

    v = OrderedDict("a" => checkbox(), "b" => 12)
    wdg = InteractBase.accordion(v, multiple = false)
    sleep(0.1)
    @test observe(wdg)[] == 1
    @test observe(wdg["options"])[] == v
    observe(wdg)[] = 2
    sleep(0.1)
    @test observe(wdg)[] == 2
    observe(wdg["options"])[] = OrderedDict("a" => 12)
    sleep(0.1)
    @test observe(wdg)[] == 2

    a = tabulator(OrderedDict("a" => 1.1, "b" => 1.2, "c" => 1.3))
    @test a[:navbar] isa InteractBase.Widget{:tabs}
    @test a[:navbar][:index][] == 1
    @test observe(a, :navbar)[] == 1
    observe(a)[] = 2
    sleep(0.1)
    @test a[:navbar][:index][] == 2
    @test observe(a, :navbar)[] == 2
    @test observe(a, "key")[] == "b"

    a = tabulator(OrderedDict("a" => 1.1, "b" => 1.2, "c" => 1.3), value = 0)
    @test a[:navbar][:index][] == 0
    @test observe(a, :key)[] == nothing

    v = OrderedDict("a" => checkbox(), "b" => 12)
    wdg = InteractBase.mask(v, multiple = true)
    sleep(0.1)
    @test observe(wdg)[] == Int[]
    @test observe(wdg["options"])[] == v
    observe(wdg)[] = [1]
    sleep(0.1)
    @test observe(wdg)[] == [1]
    observe(wdg["options"])[] = OrderedDict("a" => 12)
    sleep(0.1)
    @test observe(wdg)[] == [1]

    v = OrderedDict("a" => checkbox(), "b" => 12)
    wdg = InteractBase.mask(v; multiple = false)
    sleep(0.1)
    @test observe(wdg)[] == 1
    @test observe(wdg["options"])[] == v
    observe(wdg)[] = 2
    sleep(0.1)
    @test observe(wdg)[] == 2
    observe(wdg["options"])[] = OrderedDict("a" => 12)
    sleep(0.1)
    @test observe(wdg)[] == 2
end

@testset "node" begin
    @test InteractBase.node("a", "b") isa Node
    @test InteractBase.div("a", "b") isa Node
end

@testset "onchange" begin
    value = Observable(50)
    changes = Observable(0)
    s0 = slider(1:100, value = value, changes = changes)
    s1 = onchange(s0)
    onrelease = InteractBase.triggeredby(s0, s0[:changes])
    @test onrelease[] == 50 == s1[]
    s0[] = 12
    sleep(0.1)
    @test onrelease[] == 50 == s1[]
    s0[:changes][] += 1
    sleep(0.1)
    @test onrelease[] == 12 == s1[][]
    @test changes[] == 1
end


text:
"""
Run all tests for simulating FMU's
"""

using Test

include("tests.jl")
include("fmi-cross-check-tests.jl")

@testset "All Tests" begin
    @testset "Simulating simple FMUs" begin
        @testset "Without events" begin
            @test testHelloFMI20World()
        end;
        @testset "With events" begin
            @test testBouncingBall()
        end;
    end;

    runFMICrossTests()
end;

ext: jl

max_stars_repo_name: cuihantao/JuliaFMI

max_issues_repo_name: cuihantao/JuliaFMI

max_forks_repo_name: cuihantao/JuliaFMI


text:
using Thebes, Luxor, Colors, Rotations

function frame(scene, fn)
    background("black")
    eased_n = scene.easingfunction(fn, 0, 1, scene.framerange.stop)
    sethue("gold")
    fontsize(50)
    helloworld()
    eyepoint(400, 50, 200)
    fontface("JuliaMono-Black")
    axes3D(100)
    r = eased_n * 2π
    text3D("Hello Ole",
        Point3D(0, 0, 0),
        halign=:center,
        rotation = RotXYZ(-r, π/2 + r, π/2 + r))
end

oletext = Movie(300, 300, "oletext")
animate(oletext, [Scene(oletext, frame, 1:150, easingfunction=easeinoutcubic)],
    creategif=true,
    pathname="/tmp/oletext.gif")


text:
using GitHub, JLD, Requests

const org = "JuliaSmoothOptimizers"
const docurl = "https://$org.github.io"
const shieldurl = "https://img.shields.io"
const style = "?style=flat" # flat, flat-square or plastic

no_appveyor = ["CUTEst.jl", "HSL.jl", "LinearOperators.jl", "MUMPS.jl",
               "qr_mumps.jl"]

shield(args; label="") = "[![]($shieldurl/$args$style&label=$label)]"
shield(svc, org, pkg; label="") = shield("$svc/$org/$pkg.svg", label=label)
shield(svc, org, pkg, br; label="") = shield(svc, org, "$pkg/$br", label=label)

APV(pkg) = lowercase(replace(pkg, ".", "-"))

travis(pkg, br) = shield("travis", org, pkg, br, label="Travis_$br") * "(https://travis-ci.org/$org/$pkg)"
appveyor(pkg, br) = shield("appveyor/ci", "dpo", APV(pkg), br, label="Appveyor_$br") * "(https://ci.appveyor.com/project/dpo/$(APV(pkg))/branch/$br)"
coverage(pkg, br) = shield("coveralls", org, pkg, br, label="Coveralls_$br") * "(https://coveralls.io/github/$org/$pkg?branch=$br)"

docs(pkg) = shield("badge/docs-latest-3f51b5.svg", label="docs") * "($docurl/$pkg/latest)"

Case(s) = uppercase(s[1:1]) * lowercase(s[2:end])

function generate_packages_table(info)

  N = length(info)
  names = [info[i]["name"] for i = 1:N]
  desc  = [info[i]["desc"] for i = 1:N]

  s = ""
  s *= "This is a brief list of the installed packages and their situation.\n"
  s *= "For help, check the docs of each package,
  [gitter](https://gitter.im/JuliaSmoothOptimizers/JuliaSmoothOptimizers) for
  anyone online, or open an issue on GitHub\n"

  for i = 1:N
    pkg = names[i]
    pkg == "$org.github.io" && continue
    println("Preparing $pkg")

    github = "[![Github Stars](https://img.shields.io/github/stars/$org/$pkg.svg?style=social&label=Star)]($(info[i]["url"]))"
    s *= "## $pkg\n\n"
    s *= "$(desc[i])\n"
    s *= "\n"

    s *= github * " "
    for srv in [travis, coverage]
      for br in ["develop", "master"] ∩ info[i]["branches"]
        s *= srv(pkg, br) * " "
      end
      s *= "\n"
    end
    for srv in [appveyor]
      pkg in no_appveyor && continue
      for br in ["develop", "master"] ∩ info[i]["branches"]
        s *= srv(pkg, br) * " "
      end
      s *= "\n"
    end

    # Docs
    if statuscode(get("$docurl/$pkg/latest")) == 200
      s *= "$(docs(pkg))\n"
    end

    s *= "\n"
  end
  return s
end

function generate_home_page()
  open("src/index.md", "w") do f
    write(f, readstring("README.md"))
  end
end

"""
To be used only when information changes. It's slow because it has to
communicate with GitHub too many times
"""
function save_information()
  println("Downloading repos")
  repos = GitHub.repos(org)[1]
  n = length(repos)
  info = [Dict{String,Any}("name"=>get(repos[i].name),
                           "desc"=>get(repos[i].description),
                           "url"=>get(repos[i].html_url)
                          ) for i = 1:n]
  for i = 1:n
    println("Downloading branches of repo $(info[i]["name"])")
    br = branches(repos[i])
    info[i]["branches"] = map(x->get(x.name), br[1])
  end

  JLD.save("github_info.jld", "info", info)
  return info
end

function get_information()
  return load("github_info.jld")["info"]
end

#info = save_information()
info = get_information()
s = generate_packages_table(info)
open("src/status.md", "w") do f
  write(f, s)
end

generate_home_page()


text:
"""
    SimpleForcing{X, Y, Z, M, F, P}

A callable object that implements a 'simple' forcing function
on a field at the location `X, Y, Z`.
"""
struct SimpleForcing{X, Y, Z, M, P, N, C, F}
                forcing :: F
             parameters :: P
     field_in_signature :: Bool

    function SimpleForcing{X, Y, Z}(field_name, forcing, parameters, field_in_signature) where {X, Y, Z}

        # Deduce the container where `field_name` is found. This is only needed
        # if field_in_signature=true.
        field_container = field_name ∈ (:u, :v, :w) ? :velocities : :tracers

        return new{X, Y, Z, field_in_signature, typeof(parameters), field_name, field_container,
                   typeof(forcing)}(forcing, parameters, field_in_signature)
    end
end

"""
    SimpleForcing(forcing; parameters=nothing, field_in_signature=false)

Construct forcing for a field named `field_name` based on a `forcing` function,
with optional parameters. The keyword arguments determine the expected function
signature of `forcing`. The keyword `field_in_signature=true` specifies that the field
being forced appears in the user-defined function signature.
If `parameters` is anything other than `nothing`, it is assumed to be part of the
function signature of `forcing`.

The four possible signatures of function `forcing` are thus

* `parameters=nothing, field_in_signature=false`:

    `forcing(x, y, z, t)`

This is the function signature for default choices of `parameters` and `field_in_signature`.

* `field_in_signature=false`, specified parameters:

    `forcing(x, y, z, t, parameters)`,

where `parameters` is the object passed as keyword argument. To compile on the GPU
this must be a simple object; typically, a `NamedTuple` of floats and other constants.

* `parameters=nothing, field_in_signature=true`:

    `forcing(x, y, z, t, field)`

where `field` is the value of the field the forcing is applied to at `x, y, z`.

* `field_in_signature=true`, specified parameters:

    `forcing(x, y, z, t, field, parameters)`

Examples
========

* The simplest case: no parameters, additive forcing:

```julia
julia> const a = 2.1

julia> fun_forcing(x, y, z, t) = a * exp(z) * cos(t)

julia> u_forcing = SimpleForcing(fun_forcing)
```

* Parameterized, additive forcing:

```julia
julia> parameterized_forcing(x, y, z, t, p) = p.μ * exp(z / p.λ) * cos(p.ω * t)

julia> v_forcing = SimpleForcing(parameterized_forcing, parameters = (μ=42, λ=0.1, ω=π))
```

* Field-dependent forcing with no parameters:

```julia
julia> growth_in_sunlight(x, y, z, t, P) = exp(z) * P

julia> plankton_forcing = SimpleForcing(growth_in_sunlight, field_in_signature=true)
```

* Field-dependent forcing with parameters. This example relaxes a tracer to some reference
    linear profile.

```julia
julia> tracer_relaxation(x, y, z, t, c, p) = p.μ * exp((z + p.H) / p.λ) * (p.dCdz * z - c)

julia> c_forcing = SimpleForcing(tracer_relaxation, parameters=(μ=1/60, λ=10, H=1000, dCdz=1),
                                 field_in_signature=true)
```
"""
SimpleForcing(forcing; parameters=nothing, field_in_signature=false) =
    SimpleForcing{Cell, Cell, Cell}(:tracer, forcing, parameters, field_in_signature)

# Simple additive forcing without parameters
@inline function (f::SimpleForcing{X, Y, Z, false, <:Nothing})(i, j, k, grid, clock, state) where {X, Y, Z}

    return @inbounds f.forcing(xnode(X, i, grid),
                               ynode(Y, j, grid),
                               znode(Z, k, grid),
                               clock.time)
end

# Simple additive forcing with parameters
@inline function (f::SimpleForcing{X, Y, Z, false})(i, j, k, grid, clock, state) where {X, Y, Z}

    return @inbounds f.forcing(xnode(X, i, grid),
                               ynode(Y, j, grid),
                               znode(Z, k, grid),
                               clock.time,
                               f.parameters)
end

# Simple field-dependent forcing without parameters
@inline function (f::SimpleForcing{X, Y, Z, true, <:Nothing, N, C})(i, j, k, grid, clock, state) where {X, Y, Z, N, C}
    container = getproperty(state, C)
    field = getproperty(container, N)

    return @inbounds f.forcing(xnode(X, i, grid),
                               ynode(Y, j, grid),
                               znode(Z, k, grid),
                               field[i, j, k],
                               clock.time)
end

# Simple field-dependent forcing with parameters
@inline function (f::SimpleForcing{X, Y, Z, true, P, N, C})(i, j, k, grid, clock, state) where {X, Y, Z, P, N, C}
    container = getproperty(state, C)
    field = getproperty(container, N)

    return @inbounds f.forcing(xnode(X, i, grid),
                               ynode(Y, j, grid),
                               znode(Z, k, grid),
                               clock.time,
                               field[i, j, k],
                               f.parameters)
end


text:
using MathProgBase, GLPKMathProgInterface

import Base: rand,
             ∈

export VPolytope,
       vertices_list,
       convex_hull,
       cartesian_product,
       linear_map,
       remove_redundant_vertices

"""
    VPolytope{N<:Real} <: AbstractPolytope{N}

Type that represents a convex polytope in V-representation.

### Fields

- `vertices` -- the list of vertices
"""
struct VPolytope{N<:Real} <: AbstractPolytope{N}
    vertices::Vector{Vector{N}}
end

# constructor for a VPolytope with empty vertices list
VPolytope{N}() where {N<:Real} = VPolytope{N}(Vector{Vector{N}}())

# constructor for a VPolytope with no vertices of type Float64
VPolytope() = VPolytope{Float64}()

# constructor from a polygon in V-representation
function VPolytope(P::VPolygon, share::Bool=false)
    v = vertices_list(P)
    if !share
        v = copy(v)
    end
    return VPolytope(v)
end

# --- LazySet interface functions ---

"""
    dim(P::VPolytope)::Int

Return the dimension of a polytope in V-representation.

### Input

- `P`  -- polytope in V-representation

### Output

The ambient dimension of the polytope in V-representation.
If it is empty, the result is ``-1``.

### Examples

```jldoctest
julia> v = VPolytope();

julia> dim(v) > 0
false

julia> v = VPolytope([ones(3)])
VPolytope{Float64}(Array{Float64,1}[[1.0, 1.0, 1.0]])

julia> dim(v) == 3
true

```
"""
function dim(P::VPolytope)::Int
    return length(P.vertices) == 0 ? -1 : length(P.vertices[1])
end

"""
    σ(d::AbstractVector{N}, P::VPolytope{N}) where {N<:Real}

Return the support vector of a polytope in V-representation in a given
direction.

### Input

- `d` -- direction
- `P` -- polytope in V-representation

### Output

The support vector in the given direction.

### Algorithm

A support vector maximizes the support function.
For a polytope, the support function is always maximized in some vertex.
Hence it is sufficient to check all vertices.
"""
function σ(d::AbstractVector{N}, P::VPolytope{N}) where {N<:Real}
    # base cases
    m = length(P.vertices)
    if m == 0
        error("the support function for an empty polytope is undefined")
    elseif m == 1
        return P.vertices[1]
    end

    # evaluate support function in every vertex
    max_ρ = N(-Inf)
    max_idx = 0
    for (i, vi) in enumerate(P.vertices)
        ρ_i = dot(d, vi)
        if ρ_i > max_ρ
            max_ρ = ρ_i
            max_idx = i
        end
    end
    return P.vertices[max_idx]
end

"""
    ∈(x::AbstractVector{N}, P::VPolytope{N};
      solver=GLPKSolverLP(method=:Simplex))::Bool where {N<:Real}

Check whether a given point is contained in a polytope in vertex representation.

### Input

- `x` -- point/vector
- `P` -- polytope in vertex representation

### Output

`true` iff ``x ∈ P``.

### Algorithm

We check, using linear programming, the definition of a convex polytope that a
point is in the set if and only if it is a convex combination of the vertices.

Let ``\\{v_j\\}`` be the ``m`` vertices of `P`.
Then we solve the following ``m``-dimensional linear program.

```math
\\max 0 \\text{s.t.}
\bigwedge_{i=1}^n \\sum_{j=1}^m λ_j v_j[i] = x[i]
\\sum_{j=1}^m λ_j = 1
\bigwedge_{j=1}^m λ_j ≥ 0
```
"""
function ∈(x::AbstractVector{N}, P::VPolytope{N};
           solver=GLPKSolverLP(method=:Simplex))::Bool where {N<:Real}
    vertices = P.vertices
    m = length(vertices)

    # special cases: 0 or 1 vertex
    if m == 0
        return false
    elseif m == 1
        return x == vertices[1]
    end

    n = length(x)
    @assert n == dim(P) "a vector of length $(length(x)) cannot be " *
        "contained in a polytope of dimension $(dim(P))"

    A = Matrix{N}(undef, n+1, m)
    for j in 1:m
        v_j = vertices[j]
        # ⋀_i Σ_j λ_j v_j[i] = x[i]
        for i in 1:n
            A[i, j] = v_j[i]
        end
        # Σ_j λ_j = 1
        A[n+1, j] = one(N)
    end
    b = [x; one(N)]

    lbounds = zeros(N, m)
    ubounds = Inf
    sense = [fill('=', n); '<']
    obj = zeros(N, m)

    lp = linprog(obj, A, sense, b, lbounds, ubounds, solver)
    if lp.status == :Optimal
        return true
    elseif lp.status == :Infeasible
        return false
    end
    @assert false "LP returned status $(lp.status) unexpectedly"
end

"""
    rand(::Type{VPolytope}; [N]::Type{<:Real}=Float64, [dim]::Int=2,
         [rng]::AbstractRNG=GLOBAL_RNG, [seed]::Union{Int, Nothing}=nothing,
         [num_vertices]::Int=-1)::VPolytope{N}

Create a random polytope in vertex representation.

### Input

- `VPolytope`    -- type for dispatch
- `N`            -- (optional, default: `Float64`) numeric type
- `dim`          -- (optional, default: 2) dimension
- `rng`          -- (optional, default: `GLOBAL_RNG`) random number generator
- `seed`         -- (optional, default: `nothing`) seed for reseeding
- `num_vertices` -- (optional, default: `-1`) upper bound on the number of
                    vertices of the polytope (see comment below)

### Output

A random polytope in vertex representation.

### Algorithm

All numbers are normally distributed with mean 0 and standard deviation 1.

The number of vertices can be controlled with the argument `num_vertices`.
For a negative value we choose a random number in the range `dim:5*dim` (except
if `dim == 1`, in which case we choose in the range `1:2`).
Note that we do not guarantee that the vertices are not redundant.
"""
function rand(::Type{VPolytope};
              N::Type{<:Real}=Float64,
              dim::Int=2,
              rng::AbstractRNG=GLOBAL_RNG,
              seed::Union{Int, Nothing}=nothing,
              num_vertices::Int=-1
             )::VPolytope{N}
    rng = reseed(rng, seed)
    if num_vertices < 0
        num_vertices = (dim == 1) ? rand(1:2) : rand(dim:5*dim)
    end
    vertices = [randn(rng, N, dim) for i in 1:num_vertices]
    return VPolytope(vertices)
end

"""
    linear_map(M::AbstractMatrix{N}, P::VPolytope{N}) where {N<:Real}

Concrete linear map of a polytope in vertex representation.

### Input

- `M` -- matrix
- `P` -- polytope in vertex representation

### Output

A polytope in vertex representation.

### Algorithm

The linear map ``M`` is applied to each vertex of the given set ``P``, obtaining
a polytope in V-representation. The output type is again a `VPolytope`.
"""
function linear_map(M::AbstractMatrix{N}, P::VPolytope{N}) where {N<:Real}
    @assert dim(P) == size(M, 2) "a linear map of size $(size(M)) cannot be applied to a set of dimension $(dim(P))"
    return _linear_map_vrep(M, P)
end

@inline function _linear_map_vrep(M::AbstractMatrix{N}, P::VPolytope{N}) where {N<:Real}
    return broadcast(v -> M * v, vertices_list(P)) |> VPolytope{N}
end

"""
    translate(P::VPolytope{N}, v::AbstractVector{N}) where {N<:Real}

Translate (i.e., shift) a polytope in vertex representation by a given vector.

### Input

- `P` -- polytope in vertex representation
- `v` -- translation vector

### Output

A translated polytope in vertex representation.

### Algorithm

We add the vector to each vertex of the polytope.
"""
function translate(P::VPolytope{N}, v::AbstractVector{N}) where {N<:Real}
    @assert length(v) == dim(P) "cannot translate a $(dim(P))-dimensional " *
                                "set by a $(length(v))-dimensional vector"
    return VPolytope([x + v for x in vertices_list(P)])
end

# --- AbstractPolytope interface functions ---

"""
    vertices_list(P::VPolytope{N})::Vector{Vector{N}} where {N<:Real}

Return the list of vertices of a polytope in V-representation.

### Input

- `P` -- polytope in vertex representation

### Output

List of vertices.
"""
function vertices_list(P::VPolytope{N})::Vector{Vector{N}} where {N<:Real}
    return P.vertices
end

"""
    constraints_list(P::VPolytope{N}) where {N<:Real}

Return the list of constraints defining a polytope in V-representation.

### Input

- `P` -- polytope in V-representation

### Output

The list of constraints of the polytope.

### Algorithm

First the H-representation of ``P`` is computed, then its list of constraints
is returned.
"""
function constraints_list(P::VPolytope{N}) where {N<:Real}
    return constraints_list(tohrep(P))
end

"""
    convex_hull(P1::VPolytope{N}, P2::VPolytope{N};
                [backend]=default_polyhedra_backend(P1, N)) where {N<:Real}

Compute the convex hull of the set union of two polytopes in V-representation.

### Input

- `P1`         -- polytope
- `P2`         -- another polytope
- `backend`    -- (optional, default: `nothing`) the polyhedral
                  computations backend

### Output

The `VPolytope` obtained by the concrete convex hull of `P1` and `P2`.

### Notes

This function takes the union of the vertices of each polytope and then relies
on a concrete convex hull algorithm. For low dimensions, a specialized implementation
for polygons is used. For higher dimensions, `convex_hull` relies on the polyhedral
computations backend that can be specified using the `backend` keyword argument.

For performance reasons, it is suggested to use the `CDDLib.Library()` backend
for the `convex_hull`.
"""
function convex_hull(P1::VPolytope{N}, P2::VPolytope{N};
                     backend=nothing) where {N<:Real}
    vunion = [P1.vertices; P2.vertices]
    convex_hull!(vunion; backend=backend)
    return VPolytope(vunion)
end

"""
    cartesian_product(P1::VPolytope{N}, P2::VPolytope{N};
                      [backend]=default_polyhedra_backend(P1, N)) where {N}

Compute the Cartesian product of two polytopes in V-representation.

### Input

- `P1`         -- polytope
- `P2`         -- another polytope
- `backend`    -- (optional, default: `default_polyhedra_backend(P1, N)`) the polyhedral
                  computations backend, see
                  [Polyhedra's documentation](https://juliapolyhedra.github.io/Polyhedra.jl/latest/installation.html#Getting-Libraries-1)
                  for further information

### Output

The `VPolytope` obtained by the concrete Cartesian product of `P1` and `P2`.
"""
function cartesian_product(P1::VPolytope{N}, P2::VPolytope{N};
                           backend=default_polyhedra_backend(P1, N)) where {N}
    require(:Polyhedra; fun_name="cartesian_product")
    Pcp = vcartesianproduct(polyhedron(P1; backend=backend),
                            polyhedron(P2; backend=backend))
    return VPolytope(Pcp)
end

"""
    remove_redundant_vertices(P::VPolytope{N};
                              [backend]=nothing,
                              [solver]=nothing)::VPolytope{N} where {N<:Real}

Return the polytope obtained by removing the redundant vertices of the given polytope.

### Input

- `P`       -- polytope in vertex representation
- `backend` -- (optional, default: `nothing`) the polyhedral
               computations backend, see `default_polyhedra_backend(P1, N)` or
               [Polyhedra's documentation](https://juliapolyhedra.github.io/Polyhedra.jl/latest/installation.html#Getting-Libraries-1)
               for further information on the available backends
- `solver`  -- (optional, default: `nothing`) the linear programming
               solver used in the backend, if needed; see `default_lp_solver(N)`

### Output

A new polytope such that its vertices are the convex hull of the given polytope.

### Notes

The optimization problem associated to removing redundant vertices is handled
by `Polyhedra`. If the polyhedral computations backend requires an LP solver but
it has not been set, we use `default_lp_solver(N)` to define such solver.
Otherwise, the redundancy removal function of the polyhedral backend is used.
"""
function remove_redundant_vertices(P::VPolytope{N};
                                   backend=nothing,
                                   solver=nothing)::VPolytope{N} where {N<:Real}
    require(:Polyhedra; fun_name="remove_redundant_vertices")
    if backend == nothing
        backend = default_polyhedra_backend(P, N)
    end
    Q = polyhedron(P; backend=backend)
    if Polyhedra.supportssolver(typeof(Q))
        if solver == nothing
            solver = default_lp_solver(N)
        end
        vQ = Polyhedra.vrep(Q)
        Polyhedra.setvrep!(Q, Polyhedra.removevredundancy(vQ, solver))
    else
        Polyhedra.removevredundancy!(Q)
    end
    return VPolytope(Q)
end

"""
    tohrep(P::VPolytope{N};
           [backend]=default_polyhedra_backend(P, N)) where {N<:Real}

Transform a polytope in V-representation to a polytope in H-representation.

### Input

- `P`       -- polytope in vertex representation
- `backend` -- (optional, default: `default_polyhedra_backend(P, N)`) the
               backend for polyhedral computations

### Output

The `HPolytope` which is the constraint representation of the given polytope
in vertex representation.

### Notes

The conversion may not preserve the numeric type (e.g., with `N == Float32`)
depending on the backend.
For further information on the supported backends see
[Polyhedra's documentation](https://juliapolyhedra.github.io/Polyhedra.jl/latest/installation.html#Getting-Libraries-1).
"""
function tohrep(P::VPolytope{N};
                backend=default_polyhedra_backend(P, N)) where {N<:Real}
    vl = P.vertices
    if isempty(vl)
        return EmptySet{N}()
    end
    require(:Polyhedra; fun_name="tohrep")
    return HPolytope(polyhedron(P; backend=backend))
end

"""
    tovrep(P::VPolytope)

Return a vertex representation of the given polytope in vertex
representation (no-op).

### Input

- `P` -- polytope in vertex representation

### Output

The same polytope instance.
"""
function tovrep(P::VPolytope)
    return P
end

# ==========================================
# Lower level methods that use Polyhedra.jl
# ==========================================

function load_polyhedra_vpolytope() # function to be loaded by Requires
return quote
# see the interface file AbstractPolytope.jl for the imports

export vertices_list,
       tohrep,
       tovrep

# VPolytope from a VRep
function VPolytope(P::VRep{N}) where {N}
    vertices = Vector{Vector{N}}()
    for vi in Polyhedra.points(P)
        push!(vertices, vi)
    end
    return VPolytope(vertices)
end

"""
    polyhedron(P::VPolytope{N};
               [backend]=default_polyhedra_backend(P, N)) where {N<:Real}

Return an `VRep` polyhedron from `Polyhedra.jl` given a polytope in V-representation.

### Input

- `P`       -- polytope
- `backend` -- (optional, default: `default_polyhedra_backend(P, N)`) the polyhedral
               computations backend, see [Polyhedra's documentation](https://juliapolyhedra.github.io/Polyhedra.jl/latest/installation.html#Getting-Libraries-1)
               for further information
- `relative_dimension` -- (default, optional: `nothing`) an integer representing
                          the (relative) dimension of the polytope; this argument
                          is mandatory if the polytope is empty

### Output

A `VRep` polyhedron.

### Notes

The *relative dimension* (or just *dimension*) refers to the dimension of the set
relative to itself, independently of the ambient dimension. For example, a point
has (relative) dimension zero, and a line segment has (relative) dimension one.

In this library, `LazySets.dim` always returns the ambient dimension of the set,
such that a line segment in two dimensions has dimension two. However,
`Polyhedra.dim` will assign a dimension equal to one to a line segment
because it uses a different convention.
"""
function polyhedron(P::VPolytope{N};
                    backend=default_polyhedra_backend(P, N),
                    relative_dimension=nothing) where {N<:Real}
    if isempty(P)
        if relative_dimension == nothing
            error("the conversion to a `Polyhedra.polyhedron` requires the (relative) dimension " *
                  "of the `VPolytope` to be known, but it cannot be inferred from an empty set; " *
                  "try passing it in the keyword argument `relative_dimension`")
        end
        return polyhedron(Polyhedra.vrep(P.vertices, d=relative_dimension), backend)
    end
    return polyhedron(Polyhedra.vrep(P.vertices), backend)
end

end # quote
end # function load_polyhedra_vpolytope()


text:
import Libav
using Test

modulepath(name) = realpath(joinpath(dirname(pathof(name)),".."))

# Example
# ffmpeg.exe -y -i foobar.mp4 -f image2pipe -vcodec ppm foobar.ppm
# output:
# 0x50, 0x36, 0x0a                                          P6^J
# 0x31, 0x36, 0x30, 0x20, 0x31, 0x32, 0x30, 0x0a            160 120^J
# 0x32, 0x35, 0x35, 0x0a                                    255^J
# ...
function ppmheader_test()
    m = modulepath(Libav)
    p = joinpath(m, "test/foobar.mp4")
    d = joinpath(m, "test/foobar.ppm")
    Libav.ripvideo(p, d)
    # run(`ffmpeg.exe -y -i $p -f image2pipe -vcodec ppm $d`)
    h = Libav.ppmheader(d)
end
let h = ppmheader_test()
    @info h
end

function rgbintense_test()
    m = modulepath(Libav)
    p = joinpath(m, "test/alexa-hc-red.mp4")
    d = joinpath(m, "test/alexa-hc-red.ppm")
    Libav.ripvideo(p, d)
    rgb, noiseppm = Libav.rgbintense(d, 480, 640, 30, (2,4))
end
let (rgb, noiseppm) = rgbintense_test()
    @info noiseppm
    rm(joinpath(modulepath(Libav), "test/alexa-hc-red.ppm"))
end


text:
using Base: @propagate_inbounds

function compatible_dimensions(arr::AbstractArray, border::Inner)
    true
end

function compatible_dimensions(arr,border)
    ndims(arr) == ndims(border)
end

function convert_border_eltype(inner, border::Fill)
    T = eltype(inner)
    val::T = try
        convert(T, border.value)
    catch err
        msg = "Cannot convert elements of border=$border to eltype(inner)=$T."
        throw(ArgumentError(msg))
    end
    Fill(val, border.lo, border.hi)
end

function convert_border_eltype(inner, border)
    border
end

"""
    BorderArray(inner::AbstractArray, border::AbstractBorder) <: AbstractArray

Construct a thin wrapper around the array `inner`, with given `border`. No data is copied in the constructor, instead border values are computed on the fly in `getindex` calls. Usful for stencil computations. See also [padarray](@ref).

# Examples
```julia
julia> using ImageFiltering

julia> arr = reshape(1:6, (2,3))
2×3 reshape(::UnitRange{Int64}, 2, 3) with eltype Int64:
 1  3  5
 2  4  6

julia> BorderArray(arr, Pad((1,1)))
BorderArray{Int64,2,Base.ReshapedArray{Int64,2,UnitRange{Int64},Tuple{}},Pad{2}} with indices 0:3×0:4:
 1  1  3  5  5
 1  1  3  5  5
 2  2  4  6  6
 2  2  4  6  6

julia> BorderArray(arr, Fill(10, (2,1)))
BorderArray{Int64,2,Base.ReshapedArray{Int64,2,UnitRange{Int64},Tuple{}},Fill{Int64,2}} with indices -1:4×0:4:
 10  10  10  10  10
 10  10  10  10  10
 10   1   3   5  10
 10   2   4   6  10
 10  10  10  10  10
 10  10  10  10  10
```
"""
struct BorderArray{T,N,A,B} <: AbstractArray{T,N}
    inner::A
    border::B
    function BorderArray(inner::AbstractArray{T,N}, border::AbstractBorder) where {T,N}
        if !compatible_dimensions(inner, border)
            throw(error_bad_padding_size(inner, border))
        end
        border = convert_border_eltype(inner, border)
        A = typeof(inner)
        B = typeof(border)
        new{T,N,A,B}(inner, border)
    end
end

# these are adapted from OffsetArrays
function Base.similar(A::BorderArray, ::Type{T}, dims::Dims) where T
    B = similar(A.inner, T, dims)
end
const OffsetAxis = Union{Integer, UnitRange, Base.OneTo, IdentityUnitRange}
function Base.similar(A::BorderArray, ::Type{T}, inds::Tuple{OffsetAxis,Vararg{OffsetAxis}}) where T
    similar(A.inner, T, axes(A))
end

@inline function Base.axes(o::BorderArray)
    _outeraxes(o.inner, o.border)
end

@inline function Base.size(o::BorderArray)
    map(length, axes(o))
end

@inline function _outeraxes(arr, border::Inner)
    axes(arr)
end

@inline function _outeraxes(arr, border)
    map(axes(arr), border.lo, border.hi) do r, lo, hi
        (first(r)-lo):(last(r)+hi)
    end
end

@inline @propagate_inbounds function Base.getindex(o::BorderArray{T,N}, inds::Vararg{Int,N}) where {T,N}
    ci = CartesianIndex(inds)
    o[ci]
end

@inline function Base.getindex(o::BorderArray, ci::CartesianIndex)
    if checkbounds(Bool, o.inner, ci)
        @inbounds o.inner[ci]
    else
        @boundscheck checkbounds(o, ci)
        getindex_outer_inbounds(o.inner, o.border, ci)
    end
end

function getindex_outer_inbounds(arr, b::Fill, index)
    b.value
end

function _inner_index(arr, b::Pad, index::CartesianIndex)
    s = b.style
    inds = if s == :replicate
        map(index.I, axes(arr)) do i, r
            clamp(i, first(r), last(r))
        end
    elseif s == :circular
        map(modrange, index.I, axes(arr))
    elseif s == :reflect
        map(index.I, axes(arr)) do i, r
            if i > last(r)
                i = 2last(r) - i
            elseif i < first(r)
                i = 2first(r) - i
            end
            i
        end
    elseif s == :symmetric
        map(index.I, axes(arr)) do i, r
            if i > last(r)
                i = 2last(r) + 1 - i
            elseif i < first(r)
                i = 2first(r) -1 - i
            end
            i
        end
    else
        error("border style $(s) unrecognized")
    end
    CartesianIndex(inds)
end

function getindex_outer_inbounds(arr, b::Pad, index::CartesianIndex)
    i = _inner_index(arr, b, index)
    @inbounds arr[i]
end

function Base.checkbounds(::Type{Bool}, arr::BorderArray, index::CartesianIndex)
    map(axes(arr), index.I) do r, i
        checkindex(Bool, r, i)
    end |> all
end

function Base.copy!(dst::AbstractArray, src::BorderArray)
    axes(dst) == axes(src) || throw(DimensionMismatch("axes(dst) == axes(src) must hold."))
    _copy!(dst, src, src.border)
end
function Base.copy!(dst::AbstractArray{T,1} where T, src::BorderArray{T,1,A,B} where B where A where T)
    # fix ambiguity
    axes(dst) == axes(src) || throw(DimensionMismatch("axes(dst) == axes(src) must hold."))
    _copy!(dst, src, src.border)
end

function _copy!(dst, src, ::Inner)
    copyto!(dst, src.inner)
end

function _copy!(dst, src, ::Pad)
    inds = padindices(src.inner, src.border)
    copydata!(dst, src.inner, inds)
end

function _copy!(dst, src, ::Fill)
    try
        fill!(dst, src.border.value)
    catch
        error("Unable to fill! an array of element type $(eltype(dst)) with the value $(src.border.value). Supply an appropriate value to `Fill`, such as `zero(eltype(A))`.")
    end
    dst[axes(src.inner)...] = src.inner
    dst
end


text:
include("ACQR.jl")

include("ACQR_simula_test.jl")
include("ACQR_kfilter_test.jl")
include("ACQR_ksmoother_test.jl")
include("ACQR_kfilter_s_test.jl")
include("ACQR_ksmoother_s_test.jl")
include("ACQR_em_starting_test.jl")
include("ACQR_em_test.jl")
include("ACQR_em_s_test.jl")


text:
#demo_tr0_r.jl: tr0 read tests
#-------------------------------------------------------------------------------

using CMDimCircuits
CMDimCircuits.@using_CData()

try
import CppSimData #To access sample data ONLY!!!
catch e
	msg = "Demo requires sampled data stored in the \"CppSimData\" package."
	msg *= "\nTo execute demo, install with the following:"
	msg *= "\n\n    using Pkg; Pkg.add(PackageSpec(path=\"git://github.com/ma-laforge/CppSimData.jl\"))"

	@error msg
	rethrow(e)
end

#Get a demo display:
include(CMDimCircuits.demoplotcfgscript)
#Normally use something like:
#CMDimData.@includepkg EasyPlotInspect

#==Constants
===============================================================================#
LBLAX_AMPV = "Amplitude (V)"
LBLAX_TIME = "Time (s)"
LBLAX_RAT_100 = "Ratio (%)"
color1 = cons(:a, line = set(color=:red))
color2 = cons(:a, line = set(color=:blue))

#==Input data
===============================================================================#
testfile = "test.tr0"
filepath = joinpath(CppSimData.rootpath, "core/data", testfile)

@info("Reading $filepath...")
EDAData.open_tr0(filepath) do data
	global vin, div_val, out
	@show data
	@show names(data)
	out = read(data, "out")
	vin = read(data, "vin")
	div_val = read(data, "div_val")
end

#==Generate plot
===============================================================================#
plot = cons(:plot, nstrips=2,
	ystrip1 = set(axislabel=LBLAX_AMPV, striplabel="VCO"),
	ystrip2 = set(axislabel=LBLAX_RAT_100, striplabel="Divide Ratio"),
	xaxis = set(label=LBLAX_TIME),
)
push!(plot,
	cons(:wfrm, out, color2, label="Vout", strip=1),
	cons(:wfrm, vin, color1, label="Vcontrol", strip=1),
	cons(:wfrm, div_val, color1, strip=2),
)

pcoll = push!(cons(:plot_collection, title="EDAData Tests: tr0 Format"), plot)
	pcoll.ncolumns = 1

#==Display results in pcoll
===============================================================================#
EasyPlot.displaygui(pcoll)

#==Return pcoll to user (call evalfile(...))
===============================================================================#
pcoll


text:
immutable Gamma <: ContinuousUnivariateDistribution
    shape::Float64
    scale::Float64
    function Gamma(sh::Real, sc::Real)
        if sh > 0.0 && sc > 0.0
            new(float64(sh), float64(sc))
        else
            error("Both shape and scale must be positive")
        end
    end
end

Gamma(sh::Real) = Gamma(sh, 1.0)
Gamma() = Gamma(1.0, 1.0) # Standard exponential distribution

@_jl_dist_2p Gamma gamma

function entropy(d::Gamma)
    x = (1.0 - d.shape) * digamma(d.shape)
    x += lgamma(d.shape) + log(d.scale) + d.shape
    return x
end

insupport(d::Gamma, x::Number) = isreal(x) && isfinite(x) && 0.0 <= x

kurtosis(d::Gamma) = 6.0 / d.shape

mean(d::Gamma) = d.shape * d.scale

function mgf(d::Gamma, t::Real)
    k, theta = d.shape, d.scale
    return (1.0 - t * theta)^(-k)
end

function cf(d::Gamma, t::Real)
    k, theta = d.shape, d.scale
    return (1.0 - im * t * theta)^(-k)
end

function modes(d::Gamma)
    if d.shape >= 1.0
        d.scale * (d.shape - 1.0)
    end
end

# rand()
#
#  A simple method for generating gamma variables - Marsaglia and Tsang (2000)
#  http://www.cparity.com/projects/AcmClassification/samples/358414.pdf
#  Page 369
#  basic simulation loop for pre-computed d and c
function randg2(d::Float64, c::Float64)
    while true
        x = v = 0.0
        while v <= 0.0
            x = randn()
            v = 1.0 + c * x
        end
        v = v^3
        U = rand()
        x2 = x^2
        if U < 1.0 - 0.331 * x2^2 ||
           log(U) < 0.5 * x2 + d * (1.0 - v + log(v))
            return d * v
        end
    end
end

function rand(d::Gamma)
    dpar = (d.shape <= 1.0 ? d.shape + 1.0 : d.shape) - 1.0 / 3.0
    return d.scale *
           randg2(dpar, 1.0 / sqrt(9.0 * dpar)) *
           (d.shape > 1.0 ? 1.0 : rand()^(1.0 / d.shape))
end

function rand!(d::Gamma, A::Array{Float64})
    dpar = (d.shape <= 1.0 ? d.shape + 1.0 : d.shape) - 1.0 / 3.0
    cpar = 1.0 / sqrt(9.0 * dpar)
    for i in 1:length(A)
        A[i] = randg2(dpar, cpar)
    end
    if d.shape <= 1.0
        ainv = 1.0 / d.shape
        for i in 1:length(A)
            A[i] *= rand()^ainv
        end
    end
    return d.scale * A
end

skewness(d::Gamma) = 2.0 / sqrt(d.shape)

var(d::Gamma) = d.shape * d.scale * d.scale


text:
struct SingleVector{T} <: AbstractVector{T}
    data::T
end

Base.axes(a::SingleVector) = (Base.OneTo(1),)
Base.size(a::SingleVector) = (1,)
Base.IndexStyle(::Type{<:SingleVector}) = IndexLinear()
@inline function Base.getindex(a::SingleVector, i::Integer)
    @boundscheck if i != 1
        throw(BoundsError(a, i))
    end
    return a.data
end

Base.getindex(a::SingleVector) = a.data


text:
using Base.Test
using Laplacians

# generate chimeric graphs, and test run many routines on them

function isTree(gr::SparseMatrixCSC)
    isConnected(gr) && (nnz(gr) == 2*(gr.n-1))
end

function testSolvers(a)

    n = a.n
    excess = zeros(n); excess[1] = excess[n] = 1e-5;
    la = lap(a);
    sdd = la + spdiagm(excess);
    b = rand(n); b = b - mean(b);

    for solver in SDDSolvers
        f = solver(sdd, tol=1e-6, maxtime=0.001);
        x = f(b);
    end

    for solver in LapSolvers
        f = solver(a, tol=1e-6, maxtime=0.001);
        x = f(b);
    end

end

n = 1000
tol = 1e-11

for i in 1:100
    println("wtedChimera($n, $i)")
    gr = wtedChimera(n,i)
    @test gr.n == n
    @test isConnected(gr)
    @test isTree(kruskal(gr))
    @test isTree(randishKruskal(gr))
    @test isTree(randishPrim(gr))
    @test abs(sum(kruskal(gr)) - sum(prim(gr))) < tol
    @test sum(kruskal(gr,kind=:min)) <= sum(kruskal(gr,kind=:max))

    testSolvers(gr)
end

# every single function/method should have at least one test!

# assume each vertex pair has at most ONE unique edge
function isSubsetOfEdges(subset::SparseMatrixCSC, superset::SparseMatrixCSC, marginOfError)
  nSubsetVertices = subset.n

  for vInd in 1:nSubsetVertices
    for eSubInd in subset.colptr[vInd]:(subset.colptr[vInd+1]-1)
      for eSupInd in superset.colptr[vInd]:(superset.colptr[vInd+1]-1)
        if superset.rowval[eSupInd] == subset.rowval[eSubInd]
          if abs(superset.nzval[eSupInd] - subset.nzval[eSubInd]) > marginOfError
            println("edge: v1: ", vInd, " v2: ", subset.rowval[eSubInd], " edgeWeight: ", subset.nzval[eSubInd])
            println("\t=/= edge: v1: ", vInd, " v2: ", superset.rowval[eSupInd], " edgeWeight: ", superset.nzval[eSupInd])
            return false
          end

          break
        end
      end
    end
  end #for

  return true

end #isSubsetOfEdges

# Test Function

# uses four different types of graphs:
# dim will increase by "iterationsize" each time for each number of iterations
"""
testAKPW runs tests on akpw to make sure it runs correctly. It takes a few parameters:

startingDim is the dimension of the first set of graphs it will test (the number of nodes is dim squared)
iterationsize is the amount that dim will grow each iteration
numIterations is the number of iterations
kindVar, randomClustersVar, metisClusteringVar, exponentialX, and shuffleClustersVar each pertain to the equivalent options for akpw
marginOfError is the allowable margin of error for running edgeSubset tests (if the difference in edge weights between the tree and graph is above this margin, it will not pass the test)

In each iteration, 4 graphs (one grid, one product, one necklace and one chimera) each with
(startingDim+((numIterations-1)*iterationsize))^2 nodes are produced, and akpw is run on each. The resulting tree
is tested using two functions: isTree and isSubsetOfEdges. isTree tests to see if it is, in fact, a tree.
isSubsetOfEdges makes sure all of the edges in the tree are also edges in the original graph (this works with
  akpw! as well, if you want to test that function instead, by using the marginOfError option).

for example:

  testAKPW(4, 2, 3) should 3 sets of four (grid, product, necklace and chimera) graphs: with 16, 36, and 64 nodes.
  None should produce any errors
"""
function testAKPW(startingDim, iterationSize, numIterations; kindVar=:max, randomClustersVar=false, metisClusteringVar=false, shuffleClustersVar=true, exponentialX=true, marginOfError=.0000001)

  if startingDim < 4
    println("please choose a higher startingDim (min: 4)")
    return
  end

  dim = startingDim

  for i in 1:numIterations

    a = grid2(dim)
    n = size(a)[1]
    (ai,aj,av) = findnz(triu(a))
    gridGraph = sparse(ai,aj,rand(size(av)),n,n)
    gridGraph = gridGraph + gridGraph';

    a = productGraph(generalizedRing(dim,[1 5]), grownGraphD(dim,3));
    n = size(a)[1]
    (ai,aj,av) = findnz(triu(a))
    productGraph1 = sparse(ai,aj,rand(size(av)),n,n)
    productGraph1 = productGraph1 + productGraph1'

    a = generalizedNecklace(generalizedRing(dim,[1 5]), grownGraphD(dim,3),5);
    n = size(a)[1]
    (ai,aj,av) = findnz(triu(a))
    necklaceGraph = sparse(ai,aj,rand(size(av)),n,n)
    necklaceGraph = necklaceGraph + necklaceGraph'

    chimGraph = wtedChimera(dim^2)

    rows, columns, edgeWeights = findnz(gridGraph)
    akpwTreeGrid = akpw(gridGraph, kind=kindVar, metisClustering = metisClusteringVar, randomClusters = randomClustersVar, shuffleClusters = shuffleClustersVar, exponentialX = exponentialX)
    oldGrid = sparse(rows, columns, edgeWeights)

    rows, columns, edgeWeights = findnz(productGraph1)
    akpwTreeProduct = akpw(productGraph1, kind=kindVar, metisClustering = metisClusteringVar, randomClusters = randomClustersVar, shuffleClusters = shuffleClustersVar, exponentialX = exponentialX)
    oldProductGraph = sparse(rows, columns, edgeWeights)

    rows, columns, edgeWeights = findnz(necklaceGraph)
    akpwTreeNecklace = akpw(necklaceGraph, kind=kindVar, metisClustering = metisClusteringVar, randomClusters = randomClustersVar, shuffleClusters = shuffleClustersVar, exponentialX = exponentialX)
    oldNecklace = sparse(rows, columns, edgeWeights)

    rows, columns, edgeWeights = findnz(chimGraph)
    akpwTreeChim = akpw(chimGraph, kind=kindVar, metisClustering = metisClusteringVar, randomClusters = randomClustersVar, shuffleClusters = shuffleClustersVar, exponentialX = exponentialX)
    oldChimera = sparse(rows, columns, edgeWeights)

    if !isTree(akpwTreeGrid)
      write(STDERR, "Grid is not a tree\n")
    end

    if !isSubsetOfEdges(akpwTreeGrid, oldGrid, marginOfError)
      write(STDERR, "akpwTreeGrid is not a subset of the edges in gridGraph\n")
    end

    if !isTree(akpwTreeProduct)
      write(STDERR, "productGraph is not a tree\n")
    end

    if !isSubsetOfEdges(akpwTreeProduct, oldProductGraph, marginOfError)
      write(STDERR, "akpwTreeProduct is not a subset of the edges in productGraph\n")
    end

    if !isTree(akpwTreeNecklace)
      write(STDERR, "necklace is not a tree\n")
    end

    if !isSubsetOfEdges(akpwTreeNecklace, oldNecklace, marginOfError)
      write(STDERR, "akpwTreeNecklace is not a subset of the edges in necklaceGraph\n")
    end

    if !isTree(akpwTreeChim)
      write(STDERR, "chimera is not a tree\n")
    end

    if !isSubsetOfEdges(akpwTreeChim, oldChimera, marginOfError)
      write(STDERR, "akpwTreeChim is not a subset of the edges in chimGraph\n")
    end

    dim += iterationSize

  end#for

  println("Tests PASS!")

end #testAKPW

include("testPCG.jl")


text:
"""
    Haar{T,TM<:Manifolds.AbstractGroupManifold} <: Manifolds.MPointDistribution{TM}

Haar (invariant, i.e. uniform) measure on group manifold ``M``.

# Constructor

    Haar(M::AbstractGroupManifold, p)

Constructs a Haar distribution on the group manifold ``M`` with points represented like the
provided point ``p``, used for random sampling when implemented for compact groups.
"""
struct Haar{P,M<:Manifolds.AbstractGroupManifold} <: Manifolds.MPointDistribution{M}
    manifold::M
    point::P # point stored for random sampling
end

function Distributions.logpdf(d::Haar, p)
    T = real(
        Base.promote_eltype(Manifolds.number_eltype(p), Manifolds.number_eltype(d.point))
    )
    return zero(T)
end

Base.eltype(::Type{<:Haar{T}}) where {T} = T

function Distributions.convolve(
    d1::D1, d2::D2
) where {P1,P2,G<:Manifolds.AbstractGroupManifold,D1<:Haar{P1,G},D2<:Haar{P2,G}}
    M = d1.manifold
    # compose to promote points
    p = Manifolds.compose(M, d1.point, d2.point)
    return Haar(M, p)
end

Distributions.insupport(d::Haar, p) = Manifolds.is_manifold_point(d.manifold, p)

inversion(d::Haar) = d

## random rotation matrices

function Distributions._rand!(
    rng::AbstractRNG, ::Haar{P,M}, q::AbstractMatrix
) where {P,N,M<:SpecialOrthogonal{N}}
    n = size(q, 1)
    @assert n == N
    randn!(rng, q)
    F = qr(q)
    idx = diagind(N, N)
    q .= F.Q .* sign.(getindex.(Ref(F.R), idx'))
    if det(q) < 0 # q is orthogonal but not special; swap first two columns to make it special
        for i in 1:N
            @inbounds q[i, 1], q[i, 2] = q[i, 2], q[i, 1]
        end
    end
    return q
end

function Base.rand(
    rng::AbstractRNG, s::Random.SamplerTrivial{<:Haar{P,M}}
) where {P,M<:SpecialOrthogonal}
    d = s.self
    q = Manifolds.allocate_result(d.manifold, rand, d.point)
    return Random.rand!(rng, d, q)
end


text:
# Tests of BATSRUS.jl

using Batsrus, Test, SHA, LazyArtifacts
using Batsrus.UnitfulBatsrus, Unitful

function filecmp(path1::AbstractString, path2::AbstractString)
   stat1, stat2 = stat(path1), stat(path2)
   if !(isfile(stat1) && isfile(stat2)) || filesize(stat1) != filesize(stat2)
      return false
   end
   stat1 == stat2 && return true # same file
   open(path1, "r") do file1
      open(path2, "r") do file2
         buf1 = Vector{UInt8}(undef, 32768)
         buf2 = similar(buf1)
         while !eof(file1) && !eof(file2)
            n1 = readbytes!(file1, buf1)
            n2 = readbytes!(file2, buf2)
            n1 != n2 && return false
            0 != Base._memcmp(buf1, buf2, n1) && return false
         end
         return eof(file1) == eof(file2)
      end
   end
end

@testset "Batsrus.jl" begin
   datapath = artifact"testdata"
   @testset "Reading 1D ascii" begin
      filename = "1d__raw_2_t25.60000_n00000258.out"
      data = readdata(filename, dir=datapath, verbose=true)
      @test startswith(repr(data), "filename : 1d")
      @test Batsrus.setunits(data.head, "NORMALIZED")
      @test isa(data.head, NamedTuple)
      @test extrema(data.x) == (-127.5, 127.5)
      @test extrema(data.w) == (-0.79960780498, 1.9394335293)
   end

   @testset "Reading 2D structured binary" begin
      filename = "z=0_raw_1_t25.60000_n00000258.out"
      data = readdata(filename, dir=datapath)
      @test data.head.time == 25.6f0
      @test extrema(data.x) == (-127.5f0, 127.5f0)
      @test extrema(data.w) == (-0.79985905f0, 1.9399388f0)
   end

   @testset "Reading 2D unstructured binary" begin
      #filename = "z=0_raw_1_t25.60000_n00000258.out"
      #data = readdata(filename)
   end

   @testset "Reading 3D structured binary" begin
      filename = "3d_raw.out"
      data = readdata(filename, dir=datapath)
      plotrange = [-50.0, 50.0, -0.5, 0.5]
      X, Z, p = cutdata(data, "p"; dir="y", sequence=1, plotrange)
      @test p[1] ≈ 0.560976f0
      @test p[2] ≈ 0.53704995f0
      vars = getvars(data, ["p"])
      @test size(vars["p"]) == (8,8,8)
   end

   @testset "Log" begin
      logfilename = joinpath(datapath, "log_n000001.log")
      head, data = readlogdata(logfilename)
      @test isa(head, NamedTuple)
      @test extrema(data) == (-0.105, 258.0)
   end

   @testset "VTK" begin
      filename = joinpath(datapath, "3d_bin.dat")
      head, data, connectivity = readtecdata(filename)
      @test maximum(connectivity) ≤ head[:nNode] # check if it's read correctly
      convertTECtoVTU(head, data, connectivity)
      sha_str = bytes2hex(open(sha1, "out.vtu"))
      @test sha_str == "5b04747666542d802357dec183177f757754a254"
      rm("out.vtu")

      filetag = joinpath(datapath, "3d_mhd_amr/3d__mhd_1_t00000000_n00000000")
      batl = Batl(readhead(filetag*".info"), readtree(filetag)...)
      # local block index check
      @test Batsrus.find_grid_block(batl, [1.0, 0.0, 0.0]) == 4
      @test Batsrus.find_grid_block(batl, [100.0, 0.0, 0.0]) == -100
      connectivity = getConnectivity(batl)
      sha_str = bytes2hex(sha256(string(connectivity)))
      @test sha_str == "c6c5a65a46d86a9ba4096228c1516f89275e45e295cd305eb70c281a770ede74"
   end

   @testset "Plotting" begin
      using PyPlot
      ENV["MPLBACKEND"]="agg" # no GUI
      @testset "1D ascii" begin
         filename = "1d__raw_2_t25.60000_n00000258.out"
         data = readdata(filename, dir=datapath, verbose=false)
         plotdata(data, "p", plotmode="line")
         line = gca().lines[1]
         @test line.get_xdata() ≈ data.x
         @test line.get_ydata() ≈ data.w[:,10]
      end

      @testset "2D structured binary" begin
         filename = "z=0_raw_1_t25.60000_n00000258.out"
         data = readdata(filename, dir=datapath)
         plotdata(data, "p bx;by", plotmode="contbar streamover")
         @test isa(gca(), PyPlot.PyObject)
         contourf(data, "p")
         @test isa(gca(), PyPlot.PyObject)
         @test_throws ErrorException contourf(data, "rho", innermask=true)
         contour(data, "rho")
         @test isa(gca(), PyPlot.PyObject)
         tricontourf(data, "rho")
         @test isa(gca(), PyPlot.PyObject)
         streamplot(data, "bx;by")
         @test isa(gca(), PyPlot.PyObject)
         plt.close()
         fig = plt.figure()
         ax = fig.add_subplot(111, projection="3d")
         plot_surface(data, "rho")
         @test isa(gca(), PyPlot.PyObject)
         plt.close()
      end

      @testset "2D AMR Cartesian" begin
         filename = "bx0_mhd_6_t00000100_n00000352.out"
         data = readdata(filename, dir=datapath)
         plotdata(data, "P", plotmode="contbar")
         ax = gca()
         @test isa(ax, PyPlot.PyObject)
      end
   end

   @testset "Units" begin
      @test 1.0bu"R" > 2.0bu"Rg"
      filename = "y=0_var_1_t00000000_n00000000.out"
      data = readdata(filename, dir=datapath)
      varunit = getunit(data, "Rho")
      @test varunit == bu"amucc"
      varunit = getunit(data, "Ux")
      @test varunit == u"km/s"
      varunit = getunit(data, "Bx")
      @test varunit == u"nT"
      varunit = getunit(data, "P")
      @test varunit == u"nPa"
      varunit = getunit(data, "jx")
      @test dimension(varunit) == dimension(Unitful.A/Unitful.m^2)
      varunit = getunit(data, "ex")
      @test varunit == u"mV/m"
   end

end


text:
using ThingSpeak
using Base.Test

# write your own tests here
using ThingSpeak
function loadsecret(name="TEST_SECRET")
  if name in keys(ENV)
    println("Loading $name from environment variable")
    return ENV[name]
  else
    println("Loading $name from file")
    return open(name) do f
      strip(readline(f))
    end
  end
end

ThingSpeak.setdefaultchannel(loadsecret("THINGSPEAK_API_KEY")) #this is technically a secret but its a dedicated channel for testing
response = tsupdate({3=>3.5},status="ran runtests.jl for ThingSpeak.jl")
response = tsupdate(1.1,2.2)

tssuccess(response)


text:
abstract type AbstractColor end
abstract type AbstractRGB <: AbstractColor end
abstract type AbstractHex <: AbstractColor end

# Plots.plot(color::AbstractColor; x=0, y=0, ms=100) = Plots.scatter(
#     [x], [y], c=Hex(color).string, ms=ms, xaxis=nothing, yaxis=nothing, xticks=[], yticks=[]
# )

# ---------------------------------------------------------------------------- #
#                                      RGB                                     #
# ---------------------------------------------------------------------------- #

"""
    RGB(r, g, b)

Represents a color in RGB space
"""
mutable struct RGB <: AbstractRGB
    r:: Int64
    g:: Int64
    b:: Int64
end

"""
Alternative constructor for `RGB` passing the color values
∈ [0, 1] range as floats.
"""
RGB(r::Float64, g::Float64, b::Float64) = RGB(toint(r * 255), toint(g * 255), toint(b * 255))

RGB(rgb::Union{Tuple, AbstractVector}) = RGB(rgb...)

# ---------------------------------- methods --------------------------------- #
Base.iterate(rgb::RGB) = (rgb.r, 1)
function Base.iterate(rgb::RGB, state)
    if state==1
        return (rgb.g, 2)
    elseif state==2
        return (rgb.b, 3)
    else
        return nothing
    end
end

Base.length(rgb::RGB) = 3

Base.Vector(rgb::RGB) = [rgb.r, rgb.g, rgb.b]

# ---------------------------------------------------------------------------- #
#                                      HEX                                     #
# ---------------------------------------------------------------------------- #
"""
    Hex("#ffffff")

Base type for hex style colors
"""
mutable struct Hex <: AbstractHex
    code::String
    string::String
    function Hex(s::String)
        @assert s[1] == '#' "hex color should start with '#'"
        @assert length(s) == 7 "hex color string should have length == 7 (1+6)"
        return new(s[2:end], s)
    end
end

# ---------------------------------- methods --------------------------------- #
Base.iterate(hex::Hex) = (hex.code[1], 2)

function Base.iterate(hex::Hex, state)
    state < 7 ? (hex.code[state], state+1) : nothing
end

Base.length(hex::Hex) = 6


text:
abstract type HloOp{Opcode}; end
function opcode(op::HloOp{Opcode}) where {Opcode}
    if Opcode == :reduceWindow
        return "reduce-window"
    end
    String(Opcode)
end

function Base.convert(::Type{Type{<:XRTArray}}, s::Shape)
    XRTArray{convert(Type, XlaType(s.element_type)),
        tuple(s.dimensions...),
        length(s.dimensions)}
end

struct HloTuple <: HloOp{:tuple}
end
fill_fields!(proto::HloInstructionProto, tup::HloTuple) = nothing

struct GenericHloOp{Opcode} <: HloOp{Opcode}
    T::Type
    shape::Tuple
end

struct HloParameter <: HloOp{:parameter}
    T::Type
    shape::Tuple
    id::Int64
end
fill_fields!(proto::HloInstructionProto, p::HloParameter) = proto.parameter_number = p.id

fill_fields!(proto::HloInstructionProto, gen::GenericHloOp) = nothing

struct HloConstant <: HloOp{:constant}
    T::Type
    shape::Tuple
    value::LiteralProto
end
HloConstant(x::XLAScalar) = HloConstant(typeof(x), (), convert(LiteralProto, x))
HloConstant(a::Array{<:XLAScalar}) = HloConstant(eltype(a), size(a), convert(LiteralProto, a))
fill_fields!(proto::HloInstructionProto, c::HloConstant) = proto.literal = c.value

# Equivalent to DotDimensionNumbers, but in a format the compiler has an
# easier time understanding
struct DimNums{N1, N2, N3, N4}
    lhs_contracting_dimensions::NTuple{N1, Int64}
    rhs_contracting_dimensions::NTuple{N2, Int64}
    lhs_batch_dimensions::NTuple{N3, Int64}
    rhs_batch_dimensions::NTuple{N4, Int64}
end

function Base.convert(::Type{DotDimensionNumbers}, dds::DimNums)
    DotDimensionNumbers(
        lhs_contracting_dimensions = collect(dds.lhs_contracting_dimensions),
        rhs_contracting_dimensions = collect(dds.rhs_contracting_dimensions),
        lhs_batch_dimensions = collect(dds.lhs_batch_dimensions),
        rhs_batch_dimensions = collect(dds.rhs_batch_dimensions),
    )
end

struct HloDot <: HloOp{:dot}
    dims::DimNums
end
function fill_fields!(proto::HloInstructionProto, d::HloDot)
    proto.dot_dimension_numbers = convert(DotDimensionNumbers, d.dims)

    # Ask for increased precision
    #proto.precision_config = xla.PrecisionConfig(operand_precision = [2 for _ in proto.operand_ids])
end

# This opcode isn't really real - there's very little point to a map
# unless we're able to compile it
struct HloMap{T} <: HloOp{:map}
end
fill_fields!(proto::HloInstructionProto, s::HloMap) = nothing

struct WindowDims
    size::Int64
    stride::Int64
    padding_low::Int64
    padding_high::Int64
    window_dilation::Int64
    base_dilation::Int64
    window_reversal::Bool
end

struct ConvDimNums{N1, N2, N3}
    input_batch_dimension::Int64
    input_feature_dimension::Int64
    input_spatial_dimensions::NTuple{N1, Int64}
    kernel_input_feature_dimension::Int64
    kernel_output_feature_dimension::Int64
    kernel_spatial_dimensions::NTuple{N2, Int64}
    output_batch_dimension::Int64
    output_feature_dimension::Int64
    output_spatial_dimensions::NTuple{N3, Int64}
end

struct HloConv <: HloOp{:convolution}
    window::NTuple{N, WindowDims} where N
    dims::ConvDimNums
end

xla.WindowDimension(w::WindowDims) = xla.WindowDimension(
    size = w.size,
    stride = w.stride,
    padding_low = w.padding_low,
    padding_high = w.padding_high,
    window_dilation = w.window_dilation,
    base_dilation = w.base_dilation,
    window_reversal = w.window_reversal
)

xla.ConvolutionDimensionNumbers(cdn::ConvDimNums) = xla.ConvolutionDimensionNumbers(
    input_batch_dimension = cdn.input_batch_dimension,
    input_feature_dimension = cdn.input_feature_dimension,
    input_spatial_dimensions = collect(cdn.input_spatial_dimensions),
    kernel_input_feature_dimension = cdn.kernel_input_feature_dimension,
    kernel_output_feature_dimension = cdn.kernel_output_feature_dimension,
    kernel_spatial_dimensions = collect(cdn.kernel_spatial_dimensions),
    output_batch_dimension = cdn.output_batch_dimension,
    output_feature_dimension = cdn.output_feature_dimension,
    output_spatial_dimensions = collect(cdn.output_spatial_dimensions),
)

function fill_fields!(proto::HloInstructionProto, d::HloConv)
    window = xla.Window(dimensions = map(xla.WindowDimension, collect(d.window)))
    proto.window = window
    proto.convolution_dimension_numbers = xla.ConvolutionDimensionNumbers(d.dims)

    # Ask for increased precision
    #proto.precision_config = xla.PrecisionConfig(operand_precision = [2 for _ in proto.operand_ids])
end

struct HloReduceWindow{fT, N} <: HloOp{Symbol("reduce-window")}
    window::NTuple{N, WindowDims}
end
(::Type{HloReduceWindow{fT}})(a::NTuple{N, WindowDims}) where {fT, N} =
    HloReduceWindow{fT, N}(a)
function fill_fields!(proto::HloInstructionProto, d::HloReduceWindow)
    window = xla.Window(dimensions = map(xla.WindowDimension, collect(d.window)))
    proto.window = window
end

struct HloReduce{fT, N} <: HloOp{:reduce}
    dims::NTuple{N, Int64}
end
(::Type{HloReduce{fT}})(dims::NTuple{N, Int64}) where {fT, N} =
    HloReduce{fT, N}(dims)
function fill_fields!(proto::HloInstructionProto, d::HloReduce)
    proto.dimensions = collect(Int64, d.dims)
end

struct HloReshape{N} <: HloOp{:reshape}
    result_shape::NTuple{N, Int}
end
fill_fields!(proto::HloInstructionProto, r::HloReshape) = nothing

struct HloBroadcast{N1, N2} <: HloOp{:broadcast}
    dim_mappings::NTuple{N1, Int}
    result_shape::NTuple{N2, Int}
end
function fill_fields!(proto::HloInstructionProto, r::HloBroadcast)
    mappings = collect(Int64, r.dim_mappings)
    proto.dimensions = mappings
end

struct HloTranspose{N} <: HloOp{:transpose}
    permutation::NTuple{N, Int}
end
function fill_fields!(proto::HloInstructionProto, r::HloTranspose)
    proto.dimensions = collect(Int64, r.permutation)
end

struct HloRev{N} <: HloOp{Symbol("reverse")}
    dims::NTuple{N, Int}
end
function fill_fields!(proto::HloInstructionProto, r::HloRev)
    proto.dimensions = collect(r.dims)
end

struct HloSelectAndScatter{T, S, N} <: HloOp{Symbol("select-and-scatter")}
    window::NTuple{N, WindowDims}
end
(::Type{HloSelectAndScatter{T, S}})(a::NTuple{N, WindowDims}) where {T, S, N} =
    HloSelectAndScatter{T, S, N}(a)

function fill_fields!(proto::HloInstructionProto, r::HloSelectAndScatter)
    window = xla.Window(dimensions = map(xla.WindowDimension, collect(r.window)))
    proto.window = window
end

struct Argument
    shape::Shape
    id::Int64
end

struct HloRng <: HloOp{:rng}
    T::Type
    shape::NTuple{N, Int} where N
    kind::Int
end
function fill_fields!(proto::HloInstructionProto, r::HloRng)
    proto.distribution = Int32(r.kind)
end

struct HloGetTupleElement <: HloOp{Symbol("get-tuple-element")}
    idx::Int
end
function fill_fields!(proto::HloInstructionProto, r::HloGetTupleElement)
    proto.tuple_index = r.idx
end

let global_id = 0
    global make_id
    make_id() = (global_id += 1; global_id)
end

struct HloDynamicUpdateSlice <: HloOp{Symbol("dynamic-update-slice")}
end
fill_fields!(proto::HloInstructionProto, r::HloDynamicUpdateSlice) = nothing

struct HloDynamicSlice{N} <: HloOp{Symbol("dynamic-slice")}
    sizes::NTuple{N, Int64}
end
fill_fields!(proto::HloInstructionProto, r::HloDynamicSlice) = proto.dynamic_slice_sizes = collect(r.sizes)

const SliceDimensions = HloInstructionProto_SliceDimensions
struct SliceDim
    start::Int64
    limit::Int64
    stride::Int64
end
Base.convert(::Type{SliceDimensions}, s::SliceDim) =
    SliceDimensions(start=s.start, limit=s.limit, stride=s.stride)

struct HloSlice{N} <: HloOp{:slice}
    dims::NTuple{N, SliceDim}
end
function fill_fields!(proto::HloInstructionProto, r::HloSlice)
    proto.slice_dimensions = map(x->convert(SliceDimensions, x), collect(r.dims))
end

struct GatherDims{N1, N2, N3}
    offset_dims::NTuple{N1, Int64}
    collapsed_slice_dims::NTuple{N2, Int64}
    start_index_map::NTuple{N3, Int64}
    index_vector_dim::Int64
end

function Base.convert(::Type{GatherDimensionNumbers}, gds::GatherDims)
    GatherDimensionNumbers(
        offset_dims = collect(gds.offset_dims),
        collapsed_slice_dims = collect(gds.collapsed_slice_dims),
        start_index_map = collect(gds.start_index_map),
        index_vector_dim = gds.index_vector_dim
    )
end

struct HloGather{N} <: HloOp{:gather}
    dims::GatherDims
    slice_sizes::NTuple{N, Int64}
end
function fill_fields!(proto::HloInstructionProto, r::HloGather)
    proto.gather_dimension_numbers = convert(GatherDimensionNumbers, r.dims)
    proto.gather_slice_sizes = collect(r.slice_sizes)
end

struct HloConcatenate <: HloOp{Symbol("concatenate")}
    dim::Int64
end
fill_fields!(proto::HloInstructionProto, r::HloConcatenate) = proto.dimensions = Int64[r.dim]

struct HloInfeed{T} <: HloOp{Symbol("infeed")}
    infeed_shape::Type{T}
end
function fill_fields!(proto::HloInstructionProto, r::HloInfeed)
    # TODO: What does this field do? TF doesn't seem to use it
    proto.infeed_config = UInt8[]
end

struct HloOutfeed <: HloOp{Symbol("outfeed")}
end
function fill_fields!(proto::HloInstructionProto, r::HloOutfeed)
    # TODO: What does this field do? TF doesn't seem to use it
    proto.outfeed_config = UInt8[]
end

struct HloAfterAll <: HloOp{Symbol("after-all")}
end
fill_fields!(proto::HloInstructionProto, r::HloAfterAll) = nothing

struct ScatterDimNums{N1, N2, N3}
    update_window_dims::NTuple{N1, Int64}
    inserted_window_dims::NTuple{N2, Int64}
    scatter_dims_to_operand_dims::NTuple{N3, Int64}
    index_vector_dim::Int64
end

function Base.convert(::Type{ScatterDimensionNumbers}, sds::ScatterDimNums)
    ScatterDimensionNumbers(
        update_window_dims = collect(sds.update_window_dims),
        inserted_window_dims = collect(sds.inserted_window_dims),
        scatter_dims_to_operand_dims = collect(sds.scatter_dims_to_operand_dims),
        index_vector_dim = sds.index_vector_dim
    )
end

struct HloScatter{fT} <: HloOp{Symbol("scatter")}
    dims::ScatterDimNums
end
function fill_fields!(proto::HloInstructionProto, r::HloScatter)
    proto.scatter_dimension_numbers = convert(ScatterDimensionNumbers, r.dims)
end

struct HloIota <: HloOp{Symbol("iota")}
    result_type::Type
    result_shape::NTuple{N, Int64} where N
    iota_dimension::Int64
end
function fill_fields!(proto::HloInstructionProto, r::HloIota)
    proto.dimensions = Int[r.iota_dimension]
end

const PaddingConfigDimension = PaddingConfig_PaddingConfigDimension
struct PaddingConfigDim
    edge_padding_low::Int64
    edge_padding_high::Int64
    interior_padding::Int64
end
Base.convert(::Type{PaddingConfigDimension}, p::PaddingConfigDim) =
    PaddingConfigDimension(
        edge_padding_low = p.edge_padding_low,
        edge_padding_high = p.edge_padding_high,
        interior_padding = p.interior_padding)

struct HloPad <: HloOp{:pad}
    padding::NTuple{N, PaddingConfigDim} where {N}
end
function fill_fields!(proto::HloInstructionProto, r::HloPad)
    proto.padding_config = PaddingConfig(dimensions=[convert(PaddingConfigDimension, p) for p in r.padding])
end

struct RepGroup{N}
    replica_ids::NTuple{N, Int64}
end
Base.convert(::Type{ReplicaGroup}, x::RepGroup) = ReplicaGroup(replica_ids=collect(x.replica_ids))

struct HloCrossReplicaSum{fT} <: HloOp{Symbol("cross-replica-sum")}
    replica_groups::NTuple{N, RepGroup} where N
    all_reduce_id::Int
    barrier::String
end
function fill_fields!(proto::HloInstructionProto, r::HloCrossReplicaSum)
    proto.replica_groups = ReplicaGroup[convert(ReplicaGroup, x) for x in r.replica_groups]
    proto.all_reduce_id = r.all_reduce_id
    proto.cross_replica_sum_barrier = r.barrier
end

struct HloSort{fT} <: HloOp{:sort}
    dimension::Int
    is_stable::Bool
end

function fill_fields!(proto::HloInstructionProto, r::HloSort)
    proto.dimensions = [r.dimension]
    #proto.is_stable = r.is_stable
end

function HloInstructionProto(comp::HloComputationProto, opcode::String; id=make_id(), name=nothing)
    proto = HloInstructionProto(
        opcode=opcode,
        id=id,
        name=something(name, "comp$(comp.id)_$(opcode)$id")
    )
    push!(comp.instructions, proto)
    proto
end

function shape_from_operands(op::HloOp, operands...)
    if isa(op, HloGetTupleElement)
        xshape = operands[1].shape.tuple_shapes[op.idx+1]
    elseif isa(op, Union{HloTuple, HloSort})
        xshape = Shape(
            element_type = xla.PrimitiveType.TUPLE,
            tuple_shapes = collect(op.shape for op in operands)
        )
    elseif isa(op, Union{HloMap, HloReduce, HloReduceWindow, HloScatter, HloCrossReplicaSum})
        T, shape = shape_infer(op,
            typeof(op).parameters[1],
            map(operands) do op
                convert(Type{<:XRTArray}, op.shape)
            end...
        )
        xshape = Shape(T, shape)
    elseif isa(op, Union{HloSelectAndScatter})
        T, shape = shape_infer(op,
            typeof(op).parameters[1],
            typeof(op).parameters[2],
            map(operands) do op
                convert(Type{<:XRTArray}, op.shape)
            end...
        )
        xshape = Shape(T, shape)
    elseif isa(op, GenericHloOp)
        xshape = Shape(op.T, op.shape)
    else
        T, shape = shape_infer(op,
            map(operands) do op
                convert(Type{<:XRTArray}, op.shape)
            end...
        )
        xshape = Shape(T, shape)
    end
    return xshape
end

function HloInstructionProto(op::HloOp, @nospecialize(operands::Union{Argument, HloInstructionProto}...); id=make_id(), name="$(opcode(op))$id",
        shape=shape_from_operands(op, operands...), kwargs...)
    proto = HloInstructionProto(
        opcode = opcode(op),
        shape = shape,
        id = id,
        name = name,
        operand_ids = collect(map(x->x.id, operands)),
        kwargs...
    )
    fill_fields!(proto, op)
    proto
end

function HloInstructionProto(comp::HloComputationProto, op::HloOp, @nospecialize(args...); id=make_id(), name=nothing, shape=shape_from_operands(op, args...))
    proto = HloInstructionProto(op, args...; id=id, name=something(name, "comp$(comp.id)_$(opcode(op))$id"), shape=shape)
    push!(comp.instructions, proto)
    proto
end


text:
# Custom pretty-printing

Base.show(io::IO, ::MIME"text/plain", ::ScalarQuantity) = print(io, "Scalar")
Base.show(io::IO, ::MIME"text/plain", ::VectorQuantity) = print(io, "Vector")

function grid_description(grid)
    mins = minimum(grid)
    maxs = maximum(grid)

    join(join.(zip(mins, maxs), " … "), " × " )
end

function Base.show(io::IO, m::MIME"text/plain", f::AbstractPICDataStructure)
    show(io, m, scalarness(typeof(f)))
    data = unwrapdata(f)
    grid = getdomain(f)
    print(io, " with data:\n")
    ctx = IOContext(io, :limit=>true, :compact=>true, :displaysize => (10,50))
    Base.print_array(ctx, data)
    print(io, "\nand $(parameterless_type(typeof(grid))) grid ")
    print(io, grid_description(grid))
end


text:
struct RateParameter
end

@snn_kw mutable struct Rate{VFT=Vector{Float32}}
    param::RateParameter = RateParameter()
    N::Int32 = 100
    x::VFT = 0.5randn(N)
    r::VFT = tanh.(x)
    g::VFT = zeros(N)
    I::VFT = zeros(N)
    records::Dict = Dict()
end

"""
[Rate Neuron](https://neuronaldynamics.epfl.ch/online/Ch15.S3.html)
"""
Rate

function integrate!(p::Rate, param::RateParameter, dt::Float32)
    @unpack N, x, r, g, I = p
    @inbounds for i = 1:N
        x[i] += dt * (-x[i] + g[i] + I[i])
        r[i] = tanh(x[i]) #max(0, x[i])
    end
end


text:
# Gradient Boosted Learner
module GBBaseLearner

export GBBL,
       build_base_func,
       learner_fit,
       learner_predict

using GradientBoost.GB
using GradientBoost.LossFunctions
using GradientBoost.Util

# Gradient boosted base learner algorithm.
mutable struct GBBL{F <: AbstractFloat} <: GBAlgorithm
  loss_function::LossFunction
  sampling_rate::F
  learning_rate::F
  num_iterations::Int
  learner

  function GBBL{F}(learner; loss_function=LeastSquares(),
    sampling_rate=0.8, learning_rate=0.1,
    num_iterations=100) where F <: AbstractFloat

    new(loss_function, sampling_rate, learning_rate, num_iterations, learner)
  end
end

function GB.build_base_func(
    gb::GBBL,
    instances,
    labels,
    prev_func_pred,
    psuedo)

  # Train learner
  lf = gb.loss_function
  learner = gb.learner
  model = learner_fit(lf, learner, instances, psuedo)
  psuedo_pred = learner_predict(lf, learner, model, instances)
  model_const =
    fit_best_constant(lf, labels, psuedo, psuedo_pred, prev_func_pred)

  # Produce function that delegates prediction to model
  return (instances) ->
    model_const .* learner_predict(lf, learner, model, instances)
end

# Fits base learner.
# The learner must be instantiated within this function.
#
# @param lf Loss function (typically, this is not used).
# @param learner Base learner.
# @param instances Instances.
# @param labels Labels.
# @return Model.
function learner_fit(lf::LossFunction, learner, instances, labels)
  error("This function must be implemented by $(learner) for $(lf)")
end

# Predicts on base learner.
#
# @param lf Loss function (typically, this is not used).
# @param learner Base learner.
# @param model Model produced by base learner.
# @param instances Instances.
# @return Predictions.
function learner_predict(lf::LossFunction, learner, model, instances)
  error("This function must be implemented by $(learner) for $(lf)")
end

# Loss function fits
function fit_best_constant(lf::LeastSquares,
  labels, psuedo, psuedo_pred, prev_func_pred)

  # No refitting required
  1.0
end

function fit_best_constant(lf::LeastAbsoluteDeviation,
  labels, psuedo, psuedo_pred, prev_func_pred)

  weights = abs(psuedo_pred)
  values = labels .- prev_func_pred

  for i = 1:length(labels)
    if weights[i] != 0.0
      values[i] /= psuedo_pred[i]
    end
  end

  weighted_median(weights, values)
end
function fit_best_constant(lf::BinomialDeviance,
  labels, psuedo, psuedo_pred, prev_func_pred)

  # TODO(svs14): Add fit_best_constant (BinomialDeviance) for base learner.
  error("$(typeof(lf)) is not implemented for GBBaseLearner.")
end

end # module


text:
include("../shared_source/boilerplate.jl")

import ClimateMachine.BalanceLaws:
    vars_state,
    init_state_prognostic!,
    init_state_auxiliary!,
    compute_gradient_argument!,
    compute_gradient_flux!,
    flux_first_order!,
    flux_second_order!,
    source!,
    wavespeed,
    boundary_conditions,
    boundary_state!
import ClimateMachine.DGMethods: DGModel
import ClimateMachine.NumericalFluxes: numerical_flux_first_order!

"""
    TwoDimensionalCompressibleNavierStokesEquations <: BalanceLaw

A `BalanceLaw` for shallow water modeling.

write out the equations here

# Usage

    TwoDimensionalCompressibleNavierStokesEquations()

"""
abstract type AbstractFluid2D <: AbstractFluid end
struct Fluid2D <: AbstractFluid2D end
struct TwoDimensionalCompressibleNavierStokesEquations{
    I,
    D,
    A,
    T,
    C,
    F,
    BC,
    FT,
} <: AbstractFluid2D
    initial_value_problem::I
    domain::D
    advection::A
    turbulence::T
    coriolis::C
    forcing::F
    boundary_conditions::BC
    g::FT
    c::FT
    function TwoDimensionalCompressibleNavierStokesEquations{FT}(
        initial_value_problem::I,
        domain::D,
        advection::A,
        turbulence::T,
        coriolis::C,
        forcing::F,
        boundary_conditions::BC;
        g = FT(10), # m/s²
        c = FT(0),  #m/s
    ) where {FT <: AbstractFloat, I, D, A, T, C, F, BC}
        return new{I, D, A, T, C, F, BC, FT}(
            initial_value_problem,
            domain,
            advection,
            turbulence,
            coriolis,
            forcing,
            boundary_conditions,
            g,
            c,
        )
    end
end
CNSE2D = TwoDimensionalCompressibleNavierStokesEquations

function vars_state(m::CNSE2D, ::Prognostic, T)
    @vars begin
        ρ::T
        ρu::SVector{2, T}
        ρθ::T
    end
end

function init_state_prognostic!(m::CNSE2D, state::Vars, aux::Vars, localgeo, t)
    cnse_init_state!(m, state, aux, localgeo, t)
end

# default initial state if IVP == nothing
function cnse_init_state!(model::CNSE2D, state, aux, localgeo, t)
    ρ = 1

    state.ρ = ρ
    state.ρu = ρ * @SVector [-0, -0]
    state.ρθ = ρ

    return nothing
end

# user defined initial state
function cnse_init_state!(
    model::CNSE2D{<:InitialValueProblem},
    state,
    aux,
    localgeo,
    t,
)
    x = aux.x
    y = aux.y
    z = aux.z

    params = model.initial_value_problem.params
    ic = model.initial_value_problem.initial_conditions

    state.ρ = ic.ρ(params, x, y, z)
    state.ρu = ic.ρu(params, x, y, z)
    state.ρθ = ic.ρθ(params, x, y, z)

    return nothing
end

function vars_state(m::CNSE2D, ::Auxiliary, T)
    @vars begin
        x::T
        y::T
        z::T
    end
end

function init_state_auxiliary!(
    model::CNSE2D,
    state_auxiliary::MPIStateArray,
    grid,
    direction,
)
    init_state_auxiliary!(
        model,
        (model, aux, tmp, geom) -> cnse_init_aux!(model, aux, geom),
        state_auxiliary,
        grid,
        direction,
    )
end

function cnse_init_aux!(::CNSE2D, aux, geom)
    @inbounds begin
        aux.x = geom.coord[1]
        aux.y = geom.coord[2]
        aux.z = geom.coord[3]
    end

    return nothing
end

function vars_state(m::CNSE2D, ::Gradient, T)
    @vars begin
        ∇u::SVector{2, T}
        ∇θ::T
    end
end

function compute_gradient_argument!(
    model::CNSE2D,
    grad::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
)
    compute_gradient_argument!(model.turbulence, grad, state, aux, t)
end

compute_gradient_argument!(::LinearDrag, _...) = nothing

@inline function compute_gradient_argument!(
    ::ConstantViscosity,
    grad::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
)
    ρ = state.ρ
    ρu = state.ρu
    ρθ = state.ρθ

    u = ρu / ρ
    θ = ρθ / ρ

    grad.∇u = u
    grad.∇θ = θ

    return nothing
end

function vars_state(m::CNSE2D, ::GradientFlux, T)
    @vars begin
        ν∇u::SMatrix{3, 2, T, 6}
        κ∇θ::SVector{3, T}
    end
end

function compute_gradient_flux!(
    model::CNSE2D,
    gradflux::Vars,
    grad::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)
    compute_gradient_flux!(
        model,
        model.turbulence,
        gradflux,
        grad,
        state,
        aux,
        t,
    )
end

compute_gradient_flux!(::CNSE2D, ::LinearDrag, _...) = nothing

@inline function compute_gradient_flux!(
    ::CNSE2D,
    turb::ConstantViscosity,
    gradflux::Vars,
    grad::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)
    ν = Diagonal(@SVector [turb.ν, turb.ν, -0])
    κ = Diagonal(@SVector [turb.κ, turb.κ, -0])

    gradflux.ν∇u = -ν * grad.∇u
    gradflux.κ∇θ = -κ * grad.∇θ

    return nothing
end

@inline function flux_first_order!(
    model::CNSE2D,
    flux::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
    direction,
)

    ρ = state.ρ
    ρu = @SVector [state.ρu[1], state.ρu[2], -0]
    ρθ = state.ρθ

    ρₜ = flux.ρ
    ρuₜ = flux.ρu
    θₜ = flux.ρθ

    g = model.g

    Iʰ = @SMatrix [
        1 -0
        -0 1
        -0 -0
    ]

    flux.ρ += ρu
    flux.ρu += g * ρ^2 * Iʰ / 2

    advective_flux!(model, model.advection, flux, state, aux, t)

    return nothing
end

advective_flux!(::CNSE2D, ::Nothing, _...) = nothing

@inline function advective_flux!(
    ::CNSE2D,
    ::NonLinearAdvectionTerm,
    flux::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)
    ρ = state.ρ
    ρu = state.ρu
    ρv = @SVector [state.ρu[1], state.ρu[2], -0]
    ρθ = state.ρθ

    flux.ρu += ρv ⊗ ρu / ρ
    flux.ρθ += ρv * ρθ / ρ

    return nothing
end

function flux_second_order!(
    model::CNSE2D,
    flux::Grad,
    state::Vars,
    gradflux::Vars,
    ::Vars,
    aux::Vars,
    t::Real,
)
    flux_second_order!(model, model.turbulence, flux, state, gradflux, aux, t)
end

flux_second_order!(::CNSE2D, ::LinearDrag, _...) = nothing

@inline function flux_second_order!(
    ::CNSE2D,
    ::ConstantViscosity,
    flux::Grad,
    state::Vars,
    gradflux::Vars,
    aux::Vars,
    t::Real,
)
    flux.ρu += gradflux.ν∇u
    flux.ρθ += gradflux.κ∇θ

    return nothing
end

@inline function source!(
    model::CNSE2D,
    source::Vars,
    state::Vars,
    gradflux::Vars,
    aux::Vars,
    t::Real,
    direction,
)
    coriolis_force!(model, model.coriolis, source, state, aux, t)
    forcing_term!(model, model.forcing, source, state, aux, t)
    linear_drag!(model, model.turbulence, source, state, aux, t)

    return nothing
end

coriolis_force!(::CNSE2D, ::Nothing, _...) = nothing

@inline function coriolis_force!(
    model::CNSE2D,
    coriolis::fPlaneCoriolis,
    source,
    state,
    aux,
    t,
)
    ρu = @SVector [state.ρu[1], state.ρu[2], -0]

    # f × u
    f = [-0, -0, coriolis_parameter(model, coriolis, aux.coords)]
    id = @SVector [1, 2]
    fxρu = (f × ρu)[id]

    source.ρu -= fxρu

    return nothing
end

forcing_term!(::CNSE2D, ::Nothing, _...) = nothing

@inline function forcing_term!(
    model::CNSE2D,
    forcing::KinematicStress,
    source,
    state,
    aux,
    t,
)
    source.ρu += kinematic_stress(model, forcing, aux.coords)

    return nothing
end

linear_drag!(::CNSE2D, ::ConstantViscosity, _...) = nothing

@inline function linear_drag!(::CNSE2D, turb::LinearDrag, source, state, aux, t)
    source.ρu -= turb.λ * state.ρu

    return nothing
end

@inline wavespeed(m::CNSE2D, _...) = m.c

roe_average(ρ⁻, ρ⁺, var⁻, var⁺) =
    (sqrt(ρ⁻) * var⁻ + sqrt(ρ⁺) * var⁺) / (sqrt(ρ⁻) + sqrt(ρ⁺))

function numerical_flux_first_order!(
    ::RoeNumericalFlux,
    model::CNSE2D,
    fluxᵀn::Vars{S},
    n⁻::SVector,
    state⁻::Vars{S},
    aux⁻::Vars{A},
    state⁺::Vars{S},
    aux⁺::Vars{A},
    t,
    direction,
) where {S, A}
    numerical_flux_first_order!(
        CentralNumericalFluxFirstOrder(),
        model,
        fluxᵀn,
        n⁻,
        state⁻,
        aux⁻,
        state⁺,
        aux⁺,
        t,
        direction,
    )

    FT = eltype(fluxᵀn)

    # constants and normal vectors
    g = model.g
    @inbounds nˣ = n⁻[1]
    @inbounds nʸ = n⁻[2]

    # get minus side states
    ρ⁻ = state⁻.ρ
    @inbounds ρu⁻ = state⁻.ρu[1]
    @inbounds ρv⁻ = state⁻.ρu[2]
    ρθ⁻ = state⁻.ρθ

    u⁻ = ρu⁻ / ρ⁻
    v⁻ = ρv⁻ / ρ⁻
    θ⁻ = ρθ⁻ / ρ⁻

    # get plus side states
    ρ⁺ = state⁺.ρ
    @inbounds ρu⁺ = state⁺.ρu[1]
    @inbounds ρv⁺ = state⁺.ρu[2]
    ρθ⁺ = state⁺.ρθ

    u⁺ = ρu⁺ / ρ⁺
    v⁺ = ρv⁺ / ρ⁺
    θ⁺ = ρθ⁺ / ρ⁺

    # averages for roe fluxes
    ρ = (ρ⁺ + ρ⁻) / 2
    ρu = (ρu⁺ + ρu⁻) / 2
    ρv = (ρv⁺ + ρv⁻) / 2
    ρθ = (ρθ⁺ + ρθ⁻) / 2

    u = roe_average(ρ⁻, ρ⁺, u⁻, u⁺)
    v = roe_average(ρ⁻, ρ⁺, v⁻, v⁺)
    θ = roe_average(ρ⁻, ρ⁺, θ⁻, θ⁺)

    # normal and tangent velocities
    uₙ = nˣ * u + nʸ * v
    uₚ = nˣ * v - nʸ * u

    # differences for difference vector
    Δρ = ρ⁺ - ρ⁻
    Δρu = ρu⁺ - ρu⁻
    Δρv = ρv⁺ - ρv⁻
    Δρθ = ρθ⁺ - ρθ⁻

    Δφ = @SVector [Δρ, Δρu, Δρv, Δρθ]

    """
    # jacobian
    ∂F∂φ = [
        0 nˣ nʸ 0
        (nˣ * c^2 - u * uₙ) (uₙ + nˣ * u) (nʸ * u) 0
        (nʸ * c^2 - v * uₙ) (nˣ * v) (uₙ + nʸ * v) 0
        (-θ * uₙ) (nˣ * θ) (nʸ * θ) uₙ
    ]

    # eigen decomposition
    λ, R = eigen(∂F∂φ)
    """

    # eigen values matrix
    c = sqrt(g * ρ)
    λ = @SVector [uₙ, uₙ + c, uₙ - c, uₙ]
    Λ = Diagonal(abs.(λ))

    # eigenvector matrix
    R = @SMatrix [
        0 1 1 0
        -nʸ (u+nˣ * c) (u-nˣ * c) 0
        nˣ (v+nʸ * c) (v-nʸ * c) 0
        0 θ θ 1
    ]

    # inverse of eigenvector matrix
    R⁻¹ = @SMatrix [
        -uₚ -nʸ nˣ 0
        (c - uₙ)/(2c) nˣ/(2c) nʸ/(2c) 0
        (c + uₙ)/(2c) -nˣ/(2c) -nʸ/(2c) 0
        -θ 0 0 1
    ]

    # @test norm(R⁻¹ * R - I) ≈ 0

    # actually calculate flux
    # parent(fluxᵀn) .-= R * Λ * (R \ Δφ) / 2
    parent(fluxᵀn) .-= R * Λ * R⁻¹ * Δφ / 2

    return nothing
end

boundary_conditions(model::CNSE2D) = model.boundary_conditions

"""
    boundary_state!(nf, ::CNSE2D, args...)

applies boundary conditions for the hyperbolic fluxes
dispatches to a function in CNSEBoundaryConditions
"""
@inline function boundary_state!(nf, bc, model::CNSE2D, args...)
    return _cnse_boundary_state!(nf, bc, model, args...)
end

"""
    CNSE_boundary_state!(nf, bc::FluidBC, ::CNSE2D)

splits boundary condition application into velocity
"""
@inline function cnse_boundary_state!(nf, bc::FluidBC, m::CNSE2D, args...)
    return cnse_boundary_state!(nf, bc.momentum, m, m.turbulence, args...)
    return cnse_boundary_state!(nf, bc.temperature, m, args...)
end

include("bc_momentum.jl")
include("bc_tracer.jl")

"""
STUFF FOR ANDRE'S WRAPPERS
"""

function get_boundary_conditions(
    model::SpatialModel{BL},
) where {BL <: AbstractFluid2D}
    bcs = model.boundary_conditions

    west_east = (check_bc(bcs, :west), check_bc(bcs, :east))
    south_north = (check_bc(bcs, :south), check_bc(bcs, :north))

    return (west_east..., south_north...)
end

function DGModel(
    model::SpatialModel{BL};
    initial_conditions = nothing,
) where {BL <: AbstractFluid2D}
    params = model.parameters
    physics = model.physics

    Lˣ, Lʸ = length(model.grid.domain)
    bcs = get_boundary_conditions(model)
    FT = eltype(model.grid.numerical.vgeo)

    if !isnothing(initial_conditions)
        initial_conditions = InitialValueProblem(params, initial_conditions)
    end

    balance_law = CNSE2D{FT}(
        initial_conditions,
        (Lˣ, Lʸ),
        physics.advection,
        physics.dissipation,
        physics.coriolis,
        nothing,
        bcs,
        c = params.c,
        g = params.g,
    )

    numerical_flux_first_order = model.numerics.flux # should be a function

    rhs = DGModel(
        balance_law,
        model.grid.numerical,
        numerical_flux_first_order,
        CentralNumericalFluxSecondOrder(),
        CentralNumericalFluxGradient(),
    )

    return rhs
end


text:
__precompile__(false)

module Contextual

using TinyCassette
using MacroTools

export @context, @contextualized, @contextual, @withctx

include("shared.jl")
include("inline_contextual.jl")
include("outline_contextual.jl")

end # module


text:
using .LibClang.CEnum: name_value_pairs

cxname2clname(x::AbstractString) = "CL" * last(split(x, '_'; limit=2))
cxname2clname(x::Symbol) = cxname2clname(string(x))

"Remove FirstXX and LastXX from enumerations."
function remove_ranges(pairs)
    filter(pairs) do (sym, val)
        !occursin(r"first|last"i, string(sym))
    end
end

# each CXCursorKind enum gets a specific Julia type wrapping
# so that we can dispatch directly on node kinds.
abstract type CLCursor end

const CXCursorMap = Dict()

for (sym, val) in name_value_pairs(CXCursorKind)
    clsym = Symbol(cxname2clname(sym))
    @eval begin
        struct $clsym <: CLCursor
            cursor::CXCursor
        end
    end
end

for (sym, val) in remove_ranges(name_value_pairs(CXCursorKind))
    clsym = Symbol(cxname2clname(sym))
    @eval CXCursorMap[$sym] = $clsym
end

CLCursor(c::CXCursor) = CXCursorMap[c.kind](c)
Base.convert(::Type{T}, x::CXCursor) where {T<:CLCursor} = CLCursor(x)
Base.cconvert(::Type{CXCursor}, x::CLCursor) = x
Base.unsafe_convert(::Type{CXCursor}, x::CLCursor) = x.cursor

Base.show(io::IO, x::CXCursor) = print(io, "$(kind(x)): $(spelling(x))")
Base.show(io::IO, x::CLCursor) = print(io, "CLCursor (", typeof(x), ") ", name(x))

# each CXTypeKind enum gets a specific Julia type wrapping
# so that we can dispatch directly on node kinds.
abstract type CLType end

const CLTypeMap = Dict()

for (sym, val) in name_value_pairs(CXTypeKind)
    clsym = Symbol(cxname2clname(sym))
    @eval begin
        struct $clsym <: CLType
            type::CXType
        end
    end
end

for (sym, val) in remove_ranges(name_value_pairs(CXTypeKind))
    clsym = Symbol(cxname2clname(sym))
    @eval CLTypeMap[$sym] = $clsym
end

CLType(c::CXType) = CLTypeMap[c.kind](c)
Base.convert(::Type{T}, x::CXType) where {T<:CLType} = CLType(x)
Base.cconvert(::Type{CXType}, x::CLType) = x
Base.unsafe_convert(::Type{CXType}, x::CLType) = x.type

Base.show(io::IO, x::CXType) = print(io, "$(kind(x)): $(spelling(x))")
Base.show(io::IO, x::CLType) = print(io, "CLType (", typeof(x), ") ")

# CLToken
abstract type CLToken end

const CLTokenMap = Dict()

for (sym, val) in name_value_pairs(CXTokenKind)
    clsym = Symbol(last(split(string(sym), '_'; limit=2)))
    @eval begin
        struct $clsym <: CLToken
            token::CXToken
            kind::CXTokenKind
            text::String
        end
        CLTokenMap[$sym] = $clsym
    end
end

for (sym, val) in remove_ranges(name_value_pairs(CXTokenKind))
    clsym = Symbol(last(split(string(sym), '_'; limit=2)))
    @eval CLTokenMap[$sym] = $clsym
end

Base.cconvert(::Type{CXToken}, x::CLToken) = x
Base.unsafe_convert(::Type{CXToken}, x::CLToken) = x.token

Base.show(io::IO, x::CXToken) = print(io, "$(kind(x))")
Base.show(io::IO, x::CLToken) = print(io, typeof(x), "(\"", x.text, "\")")

# CLToken
abstract type CLComment end

const CLCommentMap = Dict()

for (sym, val) in name_value_pairs(CXCommentKind)
    clsym = Symbol(last(split(string(sym), '_'; limit=2)))
    @eval begin
        struct $clsym <: CLComment
            comment::CXComment
        end
    end
end

for (sym, val) in remove_ranges(name_value_pairs(CXCommentKind))
    clsym = Symbol(last(split(string(sym), '_'; limit=2)))
    @eval CLCommentMap[$sym] = $clsym
end

CLComment(c::CXComment) = CLCommentMap[kind(c)](c)
Base.convert(::Type{T}, x::CXComment) where {T<:CLComment} = T(x)
Base.cconvert(::Type{CXComment}, x::CLComment) = x
Base.unsafe_convert(::Type{CXComment}, x::CLComment) = x.comment

Base.show(io::IO, x::CXComment) = print(io, "$(kind(x))")
Base.show(io::IO, x::CLComment) = print(io, "CLComment (", typeof(x), ")")

getTranslationUnit(c::CLComment) = c.TranslationUnit


text:
using Documenter, Junet

makedocs()