text:
# Define a circle by radius r
circle(x, y, r) = ellipse_by_eccentricity(x, y, a = r, e = 0)

# Define an ellipse with: semi-major axis a, eccentricity e
ellipse_by_eccentricity(x, y; a, e) = ellipse_by_axis(x, y, a = a, b = √(a^2 * (1 - e^2)))

function ellipse_by_axis(x, y; a, b, num=500)
    # Define an ellipse with:
    #   semi-major axis a
    #   semi-minor axis b
    θ = LinRange(0, 2π, num)
    return @. (x + a * cos(θ), y + b * sin(θ))
end

@recipe function f(sols::AbstractArray{<:Trajectory})
    sols, sols[1].frame
end

@recipe function f(sols::AbstractArray{<:Trajectory}, frame::Abstract_ReferenceFrame)
    values = get(plotattributes, :values, nothing)
    plot_cbar = !isnothing(values) && length(values) == length(sols)

    user_xlim = pop!(plotattributes, :xlims, nothing)
    user_ylim = pop!(plotattributes, :ylims, nothing)

    legend --> false # XXX: Legend seems to mess up the plot badly!

    # Work out the maximum extent of the orbit
    xlim = (Inf, -Inf)
    ylim = (Inf, -Inf)

    for (i, sol) in enumerate(sols)
        current_xlim, current_ylim = get_margin_lims(convert_to_frame(sol, frame), plotattributes)
        if isnothing(user_xlim)
            xlim = (min(xlim[1], current_xlim[1]), max(xlim[2], current_xlim[2]))
        else
            xlim = user_xlim
        end
        if isnothing(user_ylim)
            ylim = (min(ylim[1], current_ylim[1]), max(ylim[2], current_ylim[2]))
        else
            ylim = user_ylim
        end

        @series begin
            label := false
            if plot_cbar
                seriescolor --> cgrad(:thermal)
                arrow --> false
                linewidth --> 3
                label --> false
                colorbar --> :left
                line_z --> values[i]
            else
                seriescolor --> :blue
                arrow --> true
                linewidth --> 1.5
                label --> "Trajectory"
            end
            xlims := xlim
            ylims := ylim
            if i < length(sols)
                nomodel := true
            end
            sol, frame
        end
    end
end

# ------------------
# NEW STYLE PLOTTING
# ------------------

@recipe function f(sol::Trajectory)
    trace_vars = get(plotattributes, :trace, get(plotattributes, :trace_stability, false))
    if trace_vars === false
        # Plot the frame
        if !get(plotattributes, :nomodel, false)
            @series begin
                seriesalpha := 1.0
                (sol.model, sol.frame)
            end
            framestyle --> :zerolines
        else
            framestyle --> :none
        end

        # Plot the trajectory
        vars --> (1, 2)  # (x, y)
        denseplot --> get(plotattributes, :denseplot, true)

        xlim, ylim = get_margin_lims(sol, plotattributes)
        xlims --> xlim
        ylims --> ylim

        # Formatting
        arrow --> true
        dpi --> 150
        size --> (400, 500)
        legend --> false
        xaxis --> (rotation=45)
        aspect_ratio --> 1
        ticks --> false
    end

    sol.sol
end

@recipe function f(sol::Union{DiffEqBase.ODESolution{T}, OrdinaryDiffEq.ODECompositeSolution{T}}) where {T <: ForwardDiff.Dual}
    # TODO: Fix the type dispatch on this, since it's doing piracy.

    trace_vars = get(plotattributes, :trace, false)
    trace_stability = get(plotattributes, :trace_stability, false)
    denseplot = get(plotattributes, :denseplot, true)
    plotdensity = get(plotattributes, :plotdensity, 1000)

    tspan = ForwardDiff.value.(denseplot ? range(sol.t[begin], sol.t[end], length=plotdensity) : sol.t)
    tspan_norm = @. (tspan - tspan[begin]) / (tspan[end] - tspan[begin])
    u_vals = sol.(tspan)

    if trace_vars
        STMs = hcat([reshape(v, length(v)) for v in extract_STMs(u_vals)]...)'
        @series begin
            label --> ""
            legend --> false
            tspan_norm, STMs
        end
    elseif trace_stability
        stability_indices = norm.(extract_stability(u_vals))'
        sorted_indices = hcat(map(sort, eachslice(stability_indices, dims=1))...)[4:6,:]'
        max_eigenvalues = map(maximum, eachslice(sorted_indices, dims=1))
        @series begin
            label --> ""
            legend --> false
            tspan_norm, max_eigenvalues
        end
    else
        vars = get(plotattributes, :vars, (1,2))
        ([[u[v].value for u in u_vals] for v in vars]...,)
    end
end

@recipe function f(model::Abstract_DynamicalModel, frame::Abstract_ReferenceFrame)
    nothing
end

@recipe function f(model::Abstract_DynamicalModel, frame::SynodicFrame)
    # User arguments
    nolabels = get(plotattributes, :nolabels, get(plotattributes, :nomodel, false))
    plot_libration = get(plotattributes, :libration_points, true)
    circ_props = R3BPSystemProperties(primary_body(model), secondary_body(model))

    vars = get(plotattributes, :vars, (1,2))

    @series begin
        seriestype := :shape
        seriescolor := get(plotattributes, :primary_color, :blue)
        linecolor := :black
        linewidth := 0.5
        label := nolabels ? "" : titlecase(String(primary_body(model)))
        line_z := nothing

        primary_pos = (-circ_props.μ, 0., 0.)
        ellipse_by_axis(primary_pos[vars[1]], primary_pos[vars[2]];
                        a = circ_props.R1[vars[1]] / circ_props.L,
                        b = circ_props.R1[vars[2]] / circ_props.L)
    end

    @series begin
        seriestype := :shape
        seriescolor := get(plotattributes, :secondary_color, :brown)
        linecolor := :black
        linewidth := 0.5
        label := nolabels ? "" : titlecase(String(secondary_body(model)))
        line_z := nothing

        secondary_pos = (1 - circ_props.μ, 0., 0.)
        ellipse_by_axis(secondary_pos[vars[1]], secondary_pos[vars[2]];
                        a = circ_props.R2[vars[1]] / circ_props.L,
                        b = circ_props.R2[vars[2]] / circ_props.L)
    end

    @series begin
        seriestype := :hline
        seriescolor := :black
        linestyle := :dash
        linewidth := 0.75
        line_z := nothing
        label := ""
        [0]
    end

    if get(plotattributes, :origin_secondary, true)
        @series begin
            seriestype := :vline
            seriescolor := :black
            linestyle := :dash
            linewidth := 0.75
            label := nolabels ? "" : "Origin of $(titlecase(String(secondary_body(model))))"
            line_z := nothing
            if vars[1] == 1
                [1 - circ_props.μ]
            else
                [0.]
            end
        end
    end

    @series begin
        # Find the libration points
        L = libration_points(circ_props)
        if !isnothing(L) && plot_libration
            seriestype := :scatter
            markercolor := :yellow
            markerstrokewidth := 1
            markerstrokecolor := :red
            markershape := :diamond
            markersize := 3
            label := nolabels ? "" : "Libration points"
            line_z := nothing
            [l[vars[1]] for l in L], [l[vars[2]] for l in L]
        else
            [], []
        end
    end

end

@recipe function f(sol::Trajectory, frame::Abstract_ReferenceFrame)
    convert_to_frame(sol, frame)
end

@recipe function f(state::State)
    seriestype := :scatter
    ([state.prob.u0[1]], [state.prob.u0[2]])
end

@recipe function f(state::State, frame::Abstract_ReferenceFrame)
    @series begin
        state.model, frame
    end
    convert_to_frame(state, frame)
end

function get_margin_lims(sol::Trajectory, plotattributes)
    margins = get(plotattributes, :padding, 0.10)
    a, b = get(plotattributes, :vars, (1, 2))

    # Work out the maximum extent of the orbit
    x, y = (ForwardDiff.value.(sol.sol[a,:]), ForwardDiff.value.(sol.sol[b,:]))
    xlim = (minimum(x), maximum(x))
    ylim = (minimum(y), maximum(y))

    # Add some margins to the figure
    xdiff = margins * (xlim[2] - xlim[1])
    ydiff = margins * (ylim[2] - ylim[1])

    xlims = (xlim[1] - xdiff, xlim[2] + xdiff)
    ylims = (ylim[1] - ydiff, ylim[2] + ydiff)

    return xlims, ylims
end


text:
### A Pluto.jl notebook ###
# v0.19.0

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ 1954018c-4a07-4ca5-9011-d340ccd9996f
begin
    import Pkg
    Pkg.activate() # disable Pluto's package management
end

# ╔═╡ 1b8554d8-fd89-11ea-3760-d960e6f55ebd
begin
    using SpecialFunctions
    x = 0:0.1:15
    j0, j1, j2, j3 = [besselj.(order, x) for order in 0:3]
end

# ╔═╡ 05ec2818-fd89-11ea-3da8-23ea406ba0da
begin
    ENV["GRDISPLAY"] = "pluto"
    using GR
end

# ╔═╡ 6a7c464e-fd8b-11ea-042f-cfbcaddf2d0e
begin
    using PlutoUI
    @bind xmax Slider(5:15)
end

# ╔═╡ 25f44cb2-fd89-11ea-1bfb-ad13a07d717e
plot(x, j0, x, j1, x, j2, x, j3,
     labels=["J_0", "J_1", "J_2", "J_3"], location=11, xlim=(0, xmax))

# ╔═╡ Cell order:
# ╠═1954018c-4a07-4ca5-9011-d340ccd9996f
# ╠═1b8554d8-fd89-11ea-3760-d960e6f55ebd
# ╠═05ec2818-fd89-11ea-3da8-23ea406ba0da
# ╠═6a7c464e-fd8b-11ea-042f-cfbcaddf2d0e
# ╠═25f44cb2-fd89-11ea-1bfb-ad13a07d717e

ext: jl


text:
################################################################################
# implementation of the GenML interface

inputs{F, LD}(::Type{MultilayerPerceptron{F, LD}}) = first(LD)
outputs{F, LD}(::Type{MultilayerPerceptron{F, LD}}) = last(LD)

@generated function parameters{F, LD}(::Type{MultilayerPerceptron{F, LD}})
  totalparams = 0
  for idx = 1:(length(LD) - 1)
    #just use a dummy function here, since it doesn't affect the parameter count.
    totalparams += parameters(FullyConnectedLayer{F, LD[idx], LD[idx + 1], TF.sigmoid})
  end
  :($totalparams)
end

@generated function flatten!{F, LD}(v::Vector{F}, mlp::MultilayerPerceptron{F, LD}, offset::Integer = 0)
  l = length(LD)
  #initialize various parameters.
  current_offset = 0
  code = :()
  #manually unroll this.
  for idx = 1:(length(LD) - 1)
    code = quote
      $code
      flatten!(v, mlp.layers[$idx], $current_offset)
    end
    current_offset += parameters(FullyConnectedLayer{F, LD[idx], LD[idx + 1], TF.sigmoid})
  end
  code
end

@generated function unflatten!{F, LD}(mlp::MultilayerPerceptron{F, LD}, storage::Vector{F}, offset::Integer = 0)
  l = length(LD)
  #initialize various parameters.
  current_offset = 0
  code = :( )
  #manually unroll this.
  for idx = 1:(length(LD) - 1)
    code = quote
      $code
      unflatten!(mlp.layers[$idx], storage, $current_offset)
    end
    current_offset += parameters(FullyConnectedLayer{F, LD[idx], LD[idx + 1], TF.sigmoid})
  end
  code
end

max_stars_repo_name: interplanetary-robot/GenML

max_issues_repo_name: interplanetary-robot/GenML

max_forks_repo_name: interplanetary-robot/GenML


text:
""" Graphviz support for Catlab's graph types.
"""
module GraphvizGraphs
export parse_graphviz, to_graphviz

using StaticArrays: StaticVector, SVector

using ...Present, ...Theories
using ...Graphs, ...CategoricalAlgebra.Subobjects
import ..Graphviz

# Property graphs
#################

""" Parse Graphviz output in JSON format.

Returns a property graph with graph layout and other metadata. Each node has a
position and size.

All units are in points. Note that Graphviz has 72 points per inch.
"""
function parse_graphviz(doc::AbstractDict)::AbstractPropertyGraph
  graph = doc["directed"] ? PropertyGraph{Any}() : SymmetricPropertyGraph{Any}()
  nsubgraphs = doc["_subgraph_cnt"] # Subgraphs are ignored.

  # Graph-level layout: bounds and padding.
  # It seems, but is not documented, that the first two numbers in the Graphviz
  # bounding box are always zero.
  set_gprops!(graph,
    bounds = SVector{2}(parse_vector(doc["bb"])[3:4]),
    pad = 72*parse_point(get(doc, "pad", "0,0")),
    rankdir = get(doc, "rankdir", "TB"),
  )

  # Add vertex for each Graphviz node.
  node_keys = ("id", "name", "comment", "label", "shape", "style")
  for node in doc["objects"][nsubgraphs+1:end]
    props = Dict{Symbol,Any}(
      Symbol(k) => node[k] for k in node_keys if haskey(node, k))
    props[:position] = parse_point(node["pos"])
    props[:size] = 72*SVector(parse(Float64, node["width"]),
                              parse(Float64, node["height"]))
    add_vertex!(graph, props)
  end

  # Add edge for each Graphviz edge.
  edge_keys = ("id", "comment", "label", "xlabel", "headlabel", "taillabel",
               "headport", "tailport")
  for edge in doc["edges"]
    if get(edge, "style", nothing) == "invis"
      # Omit invisible edges, which are used to tweak the layout in Graphviz.
      continue
    end
    props = Dict{Symbol,Any}(
      Symbol(k) => edge[k] for k in edge_keys if haskey(edge, k))
    props[:spline] = parse_spline(edge["pos"])
    src = edge["tail"] - nsubgraphs + 1
    tgt = edge["head"] - nsubgraphs + 1
    add_edge!(graph, src, tgt, props)
  end

  graph
end

""" Parse an array of floats in Graphviz's comma-separated format.
"""
parse_vector(s::AbstractString) = [ parse(Float64, x) for x in split(s, ",") ]

""" Parse Graphviz point.

http://www.graphviz.org/doc/info/attrs.html#k:point
"""
parse_point(s::AbstractString) = SVector{2}(parse_vector(s))

""" Parse Graphviz spline.

In Graphviz, a "spline" is a cubic B-spline of overlapping cubic Bezier curves.
It consists of 3n+1 points, where n is the number of Bezier curves.

http://www.graphviz.org/doc/info/attrs.html#k:splineType
https://web.archive.org/web/20170418034924/http://www.graphviz.org/content/how-convert-b-spline-bezier
"""
function parse_spline(spline::AbstractString)
  points = StaticVector{2,Float64}[]
  start, stop = nothing, nothing
  for s in split(spline, " ")
    if startswith(s, "s,"); start = parse_point(s[3:end])
    elseif startswith(s, "e,"); stop = parse_point(s[3:end])
    else push!(points, parse_point(s)) end
  end
  # Prefer explicit start or end points to the spline start and end points.
  # Thus, endpoints may pass into the node shape but should not fall short.
  if !isnothing(start); points[0] = start end
  if !isnothing(stop); points[end] = stop end
  points
end

""" Convert a property graph to a Graphviz graph.

This method is usually more convenient than direct AST manipulation for creating
simple Graphviz graphs. For more advanced features, like nested subgraphs, you
must use the Graphviz AST.
"""
function to_graphviz(g::AbstractPropertyGraph)::Graphviz.Graph
  gv_name(v::Int) = "n$v"
  gv_path(e::Int) = [gv_name(src(g,e)), gv_name(tgt(g,e))]

  # Add Graphviz node for each vertex in property graph.
  stmts = Graphviz.Statement[]
  for v in vertices(g)
    push!(stmts, Graphviz.Node(gv_name(v), vprops(g, v)))
  end

  # Add Graphviz edge for each edge in property graph.
  is_directed = !(g isa SymmetricPropertyGraph)
  for e in edges(g)
    # In undirected case, only include one edge from each pair.
    if is_directed || (e <= inv(g,e))
      push!(stmts, Graphviz.Edge(gv_path(e), eprops(g, e)))
    end
  end

  attrs = gprops(g)
  Graphviz.Graph(
    name = get(attrs, :name, "G"),
    directed = is_directed,
    prog = get(attrs, :prog, is_directed ? "dot" : "neato"),
    stmts = stmts,
    graph_attrs = Graphviz.as_attributes(get(attrs, :graph, Dict())),
    node_attrs = Graphviz.as_attributes(get(attrs, :node, Dict())),
    edge_attrs = Graphviz.as_attributes(get(attrs, :edge, Dict())),
  )
end

# Graphs
########

# FIXME: Should be an abstract type after StructACSets refactor.
const IsGraph = Union{AbstractGraph, AbstractSymmetricGraph,
                      AbstractReflexiveGraph, AbstractSymmetricReflexiveGraph}

""" Convert a graph to a Graphviz graph.

A simple default style is applied. For more control over the visual appearance,
first convert the graph to a property graph, define the Graphviz attributes as
needed, and then convert to a Graphviz graph.
"""
function to_graphviz(g::IsGraph; kw...)
  to_graphviz(to_graphviz_property_graph(g; kw...))
end

function to_graphviz_property_graph(g::AbstractGraph;
    prog::AbstractString="dot", graph_attrs::AbstractDict=Dict(),
    node_attrs::AbstractDict=Dict(), edge_attrs::AbstractDict=Dict(),
    node_labels::Bool=false, edge_labels::Bool=false)
  node_labeler(v) = Dict(:label => node_labels ? string(v) : "")
  edge_labeler(e) = if edge_labels
    Dict(:label => string(e))
  else
    Dict{Symbol,String}()
  end
  PropertyGraph{Any}(g, node_labeler, edge_labeler;
    prog = prog,
    graph = merge!(Dict(:rankdir => "LR"), graph_attrs),
    node = merge!(default_node_attrs(node_labels), node_attrs),
    edge = merge!(Dict(:arrowsize => "0.5"), edge_attrs),
  )
end

default_node_attrs(show_labels::Bool) = Dict(
  :shape => show_labels ? "circle" : "point",
  :width => "0.05",
  :height => "0.05",
  :margin => "0",
)

# Symmetric graphs
##################

function to_graphviz_property_graph(g::AbstractSymmetricGraph;
    prog::AbstractString="neato", graph_attrs::AbstractDict=Dict(),
    node_attrs::AbstractDict=Dict(), edge_attrs::AbstractDict=Dict(),
    node_labels::Bool=false, edge_labels::Bool=false)
  node_labeler(v) = Dict(:label => node_labels ? string(v) : "")
  edge_labeler(e) = if edge_labels
    e′ = inv(g,e)
    Dict(:label => "($(min(e,e′)),$(max(e,e′)))")
  else
    Dict{Symbol,String}()
  end
  SymmetricPropertyGraph{Any}(g, node_labeler, edge_labeler;
    prog = prog,
    graph = graph_attrs,
    node = merge!(default_node_attrs(node_labels), node_attrs),
    edge = merge!(Dict(:len => "0.5"), edge_attrs),
  )
end

# Reflexive graphs
##################

function to_graphviz_property_graph(g::AbstractReflexiveGraph;
    prog::AbstractString="dot", graph_attrs::AbstractDict=Dict(),
    node_attrs::AbstractDict=Dict(), edge_attrs::AbstractDict=Dict(),
    node_labels::Bool=false, edge_labels::Bool=false,
    show_reflexive::Bool=false)
  pg = PropertyGraph{Any}(; prog = prog,
    graph = graph_attrs,
    node = merge!(default_node_attrs(node_labels), node_attrs),
    edge = merge!(Dict(:arrowsize => "0.5"), edge_attrs),
  )
  for v in vertices(g)
    add_vertex!(pg, label=node_labels ? string(v) : "")
  end
  reflexive_edges = Set(refl(g))
  for e in edges(g)
    is_reflexive = e ∈ reflexive_edges
    if show_reflexive || !is_reflexive
      e′ = add_edge!(pg, src(g,e), tgt(g,e))
      if is_reflexive; set_eprop!(pg, e′, :style, "dashed") end
      if edge_labels; set_eprop!(pg, e′, :label, string(e)) end
    end
  end
  pg
end

# Symmetric reflexive graphs
############################

function to_graphviz_property_graph(g::AbstractSymmetricReflexiveGraph;
    prog::AbstractString="neato", graph_attrs::AbstractDict=Dict(),
    node_attrs::AbstractDict=Dict(), edge_attrs::AbstractDict=Dict(),
    node_labels::Bool=false, edge_labels::Bool=false,
    show_reflexive::Bool=false)
  pg = SymmetricPropertyGraph{Any}(; prog = prog,
    graph = graph_attrs,
    node = merge!(default_node_attrs(node_labels), node_attrs),
    edge = merge!(Dict(:len => "0.5"), edge_attrs),
  )
  for v in vertices(g)
    add_vertex!(pg, label=node_labels ? string(v) : "")
  end
  reflexive_edges = Set(refl(g))
  for e in edges(g)
    is_reflexive = e ∈ reflexive_edges
    if is_reflexive && show_reflexive
      e′ = add_edge!(pg, src(g,e), tgt(g,e), style="dashed")
      if edge_labels; set_eprop!(pg, e′, :label, string(e)) end
    elseif !is_reflexive && e <= inv(g,e)
      e′ = add_edge!(pg, src(g,e), tgt(g,e))
      if edge_labels; set_eprop!(pg, e′, :label, "($e,$(inv(g,e)))") end
    end
  end
  pg
end

# Half-edge graphs
##################

to_graphviz(g::AbstractHalfEdgeGraph; kw...) =
  to_graphviz(to_graphviz_property_graph(g; kw...))

function to_graphviz_property_graph(g::AbstractHalfEdgeGraph;
    prog::AbstractString="neato", graph_attrs::AbstractDict=Dict(),
    node_attrs::AbstractDict=Dict(), edge_attrs::AbstractDict=Dict(),
    node_labels::Bool=false, edge_labels::Bool=false)
  pg = SymmetricPropertyGraph{Any}(; prog = prog,
    graph = graph_attrs,
    node = merge!(default_node_attrs(node_labels), node_attrs),
    edge = merge!(Dict(:len => "0.5"), edge_attrs),
  )
  for v in vertices(g)
    add_vertex!(pg, label=node_labels ? string(v) : "")
  end
  for e in half_edges(g)
    if e == inv(g,e)
      # Dangling half-edge: add an invisible vertex.
      v = add_vertex!(pg, style="invis", shape="none", label="")
      e′ = add_edge!(pg, vertex(g,e), v)
      if edge_labels; set_eprop!(pg, e′, :label, string(e)) end
    elseif e < inv(g,e)
      # Pair of distict half-edges: add a standard edge.
      e′ = add_edge!(pg, vertex(g,e), vertex(g,inv(g,e)))
      if edge_labels; set_eprop!(pg, e′, :label, "($e,$(inv(g,e)))") end
    end
  end
  pg
end

# Subgraphs
###########

to_graphviz(subgraph::Subobject{<:IsGraph}; kw...) =
  to_graphviz(to_graphviz_property_graph(subgraph; kw...))

function to_graphviz_property_graph(
    subgraph::Subobject{<:Union{AbstractGraph,AbstractSymmetricGraph}};
    subgraph_node_attrs::AbstractDict=Dict(:color => "cornflowerblue"),
    subgraph_edge_attrs::AbstractDict=Dict(:color => "cornflowerblue"), kw...)
  pg = to_graphviz_property_graph(ob(subgraph); kw...)
  f = hom(subgraph)
  for v in vertices(dom(f))
    set_vprops!(pg, f[:V](v), subgraph_node_attrs)
  end
  for e in edges(dom(f))
    set_eprops!(pg, f[:E](e), subgraph_edge_attrs)
  end
  pg
end

# Schemas
#########

to_graphviz(pres::Presentation{Schema}) =
  to_graphviz(to_graphviz_property_graph(pres))

function to_graphviz_property_graph(pres::Presentation{Schema})
  obs,homs,attr_types,attrs = generators.(Ref(pres), [:Ob,:Hom,:AttrType,:Attr])
  g = PropertyGraph{Any}()

  ob_vertices = add_vertices!(g,length(obs))
  for (v, ob) in zip(ob_vertices, obs)
    set_vprop!(g,v,:label,string(nameof(ob)))
    set_vprop!(g,v,:shape,"plain")
    set_vprop!(g,v,:margin,"2")
  end

  attr_vertices = add_vertices!(g,length(attr_types))
  for (v, attr_type) in zip(attr_vertices, attr_types)
    set_vprop!(g,v,:label,string(nameof(attr_type)))
  end

  hom_edges = add_edges!(g,
    ob_vertices[generator_index.(Ref(pres), nameof.(dom.(homs)))],
    ob_vertices[generator_index.(Ref(pres), nameof.(codom.(homs)))])
  for (e, hom) in zip(hom_edges, homs)
    set_eprop!(g,e,:label,string(nameof(hom)))
    set_eprop!(g,e,:len,"2")
  end

  attr_edges = add_edges!(g,
    ob_vertices[generator_index.(Ref(pres), nameof.(dom.(attrs)))],
    attr_vertices[generator_index.(Ref(pres), nameof.(codom.(attrs)))])
  for (e, attr) in zip(attr_edges, attrs)
    set_eprop!(g,e,:label,string(nameof(attr)))
    set_eprop!(g,e,:len,"2")
  end

  set_gprop!(g,:graph,Dict(:rankdir => "LR"))
  set_gprop!(g,:prog,"neato")
  g
end

end


text:
#
# Print
#

Base.show(io::IO, c::Contract) = _show_tree(io, c)
Base.show(io::IO, o::Observable) = _show_tree(io, o)
Base.show(io::IO, o::BufferedObservable) = print(io, "BufferedObservable(#$(repr(UInt64(pointer_from_objref(o)))))")

function Base.show(io::IO, sc::SpotCurrency; color::Bool=true)
    if color
        printstyled(io, "SpotCurrency(", sc.currency, ")", color=:yellow)
    else
        print(io, "SpotCurrency(", sc.currency, ")")
    end
end

function Base.show(io::IO, disc::DiscountFactor; color::Bool=true)
    if color
        printstyled(io, "DiscountFactor(", disc.sym, ", ", disc.maturity, ")", color=:yellow)
    else
        print(io, "DiscountFactor(", disc.sym, ", ", disc.maturity, ")")
    end
end

for T in (:FixedCashRiskFactor, :FixedNonCashRiskFactor)
    @eval begin
        function Base.show(io::IO, i::$T; color::Bool=true)
            if color
                printstyled(io, $T,"(", i.rf, " = ", i.val, ")", color=:yellow)
            else
                print(io, $T,"(", i.rf, " = ", i.val, ")")
            end
        end
    end
end

_show_root(io::IO, x) = print(io, x)
_children(x) = ()

_show_root(io::IO, c::Contract) = printstyled(io, string(typeof(c).name.name), color=:magenta)
_show_root(io::IO, evn::FixedIncomeEvent) = printstyled(io, string(typeof(evn).name.name), "{:", event_symbol(evn), "}", color=:magenta)
_children(c::Contract) = (getfield(c,i) for i = 1:nfields(c))

_show_root(io::IO, o::Observable) = printstyled(io, string(typeof(o).name.name), color=:cyan)
_children(o::Observable) = (getfield(o,i) for i = 1:nfields(o))

_show_root(io::IO, o::Konst) = printstyled(io, string(o.val), color=:yellow)
_children(o::Konst) = ()

#=
function _show_root(io::IO, o::LiftObs)
    printstyled(io, "{", color=:cyan)
    print(io,o.f)
    printstyled(io, "}", color=:cyan)
end
_children(o::LiftObs) = o.a
=#

function _show_tree(io::IO, x, indent_root="", indent_leaf="")
    print(io, indent_root)
    _show_root(io, x)
    print(io, '\n')
    cs = _children(x)
    for (i,c) in enumerate(cs)
        if i < length(cs)
            _show_tree(io, c, indent_leaf*" ├─",indent_leaf*" │ ")
        else
            _show_tree(io, c, indent_leaf*" └─",indent_leaf*"   ")
        end
    end
end

# Bash colors
const COLOR_RED = 196
const COLOR_PURPLE = 165
const COLOR_BLUE = 20
const COLOR_GREY = 237
const COLOR_GREEN = 70

function pretty_formula(x::LiftObs)
    str_f = nameof(x.f)
    str_vp = pretty_formula(x.o)
    "$str_f($str_vp)"
end

function pretty_formula(x::LiftObs2)
    str_vp1 = pretty_formula(x.o1)
    str_vp2 = pretty_formula(x.o2)
    str_f = nameof(x.f)
    #str_new_line = ((x.f in [+, -]) ? "\n" : "")
    str_new_line = ""

    if x.f in [+, -, /, *, ^]
        result = "($str_vp1 $str_f$str_new_line $str_vp2)"
    else
        result = "$str_f($str_vp1, $str_vp2)"
    end
    result
end

pretty_formula(x::Konst) = "$(x.val)"
function pretty_formula(x::Konst{T}) where {T<:AbstractFloat}
    @sprintf("%15.6f", x.val)
end

# Risk Factors
function pretty_formula(x::SpotCurrency)
    io = IOBuffer()
    printstyled(io, x.currency, color=COLOR_PURPLE)
    print(io, "@")
    printstyled(io, "t", color=COLOR_BLUE)
    String(take!(io))
end

function pretty_formula(x::Stock)
    io = IOBuffer()
    printstyled(io, x.ticker, color=COLOR_GREY)
    print(io, "@")
    printstyled(io, "t", color=COLOR_BLUE)
    String(take!(io))
end

function pretty_formula(x::DiscountFactor)
    io = IOBuffer()
    printstyled(io, x.sym, color=COLOR_GREEN)
    print(io, "@")
    printstyled(io, "t", color=COLOR_BLUE)
    print(io, "→")
    printstyled(io, x.maturity, color=COLOR_BLUE)
    String(take!(io))
end

function pretty_formula(x::DiscountFactorForward)
    io = IOBuffer()
    printstyled(io, x.sym, color=COLOR_GREEN)
    print(io, "@")
    printstyled(io, x.start_date, color=COLOR_BLUE)
    print(io, "→")
    printstyled(io, x.end_date, color=COLOR_BLUE)
    String(take!(io))
end


text:
function dot_simd(a::AbstractVector, b::AbstractVector)
    s = zero(eltype(a))
    @fastmath @inbounds @simd for i ∈ eachindex(a)
        s += a[i]' * b[i]
    end
    s
end
function cdot_mat(ca::AbstractVector{Complex{T}}, cb::AbstractVector{Complex{T}}) where {T}
    a = reinterpret(reshape, T, ca)
    b = reinterpret(reshape, T, cb)
    re = zero(T); im = zero(T)
    @turbo for i ∈ axes(a,2)
        re += a[1,i] * b[1,i] + a[2,i] * b[2,i]
        im += a[1,i] * b[2,i] - a[2,i] * b[1,i]
    end
    Complex(re, im)
end
function cdot_swizzle(ca::AbstractVector{Complex{T}}, cb::AbstractVector{Complex{T}}) where {T}
  a = reinterpret(T, ca)
  b = reinterpret(T, cb)
  reim = Vec(zero(T),zero(T))
  @turbo for i ∈ eachindex(a)
    reim = vfmsubadd(vmovsldup(a[i]), b[i], vfmsubadd(vmovshdup(a[i]), vpermilps177(b[i]), reim))
  end
  Complex(reim(1), reim(2))
end
function cdot_affine(ca::AbstractVector{Complex{T}}, cb::AbstractVector{Complex{T}}) where {T}
    a = reinterpret(T, ca);
    b = reinterpret(T, cb);
    re = zero(T); im = zero(T)
    # with a multiplier, we go from `i = 1 -> 2i = 2` to `i = 0 -> 2i = 0
    # 2(i+1-1) = 2i + 2 - 2, so....
    @turbo for i ∈ 1:length(a)>>>1
        re += a[2i-1] * b[2i-1] + a[2i] * b[2i  ]
        im += a[2i-1] * b[2i  ] - a[2i] * b[2i-1]
    end
    Complex(re, im)
end
function cdot_stride(ca::AbstractVector{Complex{T}}, cb::AbstractVector{Complex{T}}) where {T}
    a = reinterpret(T, ca);
    b = reinterpret(T, cb);
    re = zero(T); im = zero(T)
    @turbo for i ∈ 1:2:length(a)
        re += a[i] * b[i  ] + a[i+1] * b[i+1]
        im += a[i] * b[i+1] - a[i+1] * b[i  ]
    end
    Complex(re, im)
end
function qdot_simd(x::AbstractVector{NTuple{4,T}}, y::AbstractVector{NTuple{4,T}}) where {T}
    a = zero(T)
    b = zero(T)
    c = zero(T)
    d = zero(T)
    @fastmath @inbounds @simd for i ∈ eachindex(x)
        a₁, b₁, c₁, d₁ = x[i]
        a₂, b₂, c₂, d₂ = y[i]
        a += a₁*a₂ + b₁*b₂ + c₁*c₂ + d₁*d₂
        b += a₁*b₂ - b₁*a₂ - c₁*d₂ + d₁*c₂
        c += a₁*c₂ + b₁*d₂ - c₁*a₂ - d₁*b₂
        d += a₁*d₂ - b₁*c₂ + c₁*b₂ - d₁*a₂
    end
    (a,b,c,d)
end
function qdot_mat(x::AbstractMatrix, y::AbstractMatrix)
    a = zero(eltype(x))
    b = zero(eltype(x))
    c = zero(eltype(x))
    d = zero(eltype(x))
    @turbo for i ∈ axes(x,2)
        a₁ = x[1,i]
        b₁ = x[2,i]
        c₁ = x[3,i]
        d₁ = x[4,i]
        a₂ = y[1,i]
        b₂ = y[2,i]
        c₂ = y[3,i]
        d₂ = y[4,i]
        a += a₁*a₂ + b₁*b₂ + c₁*c₂ + d₁*d₂
        b += a₁*b₂ - b₁*a₂ - c₁*d₂ + d₁*c₂
        c += a₁*c₂ + b₁*d₂ - c₁*a₂ - d₁*b₂
        d += a₁*d₂ - b₁*c₂ + c₁*b₂ - d₁*a₂
    end
    (a,b,c,d)
end
function qdot_affine(x::AbstractVector, y::AbstractVector)
    a = zero(eltype(x))
    b = zero(eltype(x))
    c = zero(eltype(x))
    d = zero(eltype(x))
    @turbo for i ∈ 1:length(x)>>2
        a₁ = x[4i-3]
        b₁ = x[4i-2]
        c₁ = x[4i-1]
        d₁ = x[4i]
        a₂ = y[4i-3]
        b₂ = y[4i-2]
        c₂ = y[4i-1]
        d₂ = y[4i]
        a += a₁*a₂ + b₁*b₂ + c₁*c₂ + d₁*d₂
        b += a₁*b₂ - b₁*a₂ - c₁*d₂ + d₁*c₂
        c += a₁*c₂ + b₁*d₂ - c₁*a₂ - d₁*b₂
        d += a₁*d₂ - b₁*c₂ + c₁*b₂ - d₁*a₂
    end
    (a,b,c,d)
end
function qdot_stride(x::AbstractVector, y::AbstractVector)
    a = zero(eltype(x))
    b = zero(eltype(x))
    c = zero(eltype(x))
    d = zero(eltype(x))
    @turbo for i ∈ 1:4:length(x)
        a₁ = x[i]
        b₁ = x[i+1]
        c₁ = x[i+2]
        d₁ = x[i+3]
        a₂ = y[i]
        b₂ = y[i+1]
        c₂ = y[i+2]
        d₂ = y[i+3]
        a += a₁*a₂ + b₁*b₂ + c₁*c₂ + d₁*d₂
        b += a₁*b₂ - b₁*a₂ - c₁*d₂ + d₁*c₂
        c += a₁*c₂ + b₁*d₂ - c₁*a₂ - d₁*b₂
        d += a₁*d₂ - b₁*c₂ + c₁*b₂ - d₁*a₂
    end
    (a,b,c,d)
end
function cmatmul_array!(Cc::AbstractMatrix{Complex{T}}, Ac::AbstractMatrix{Complex{T}}, Bc::AbstractMatrix{Complex{T}}) where {T}
  C = reinterpret(reshape, Float64, Cc);
  A = reinterpret(reshape, Float64, Ac);
  B = reinterpret(reshape, Float64, Bc);
  @turbo for n ∈ indices((C,B),3), m ∈ indices((C,A),2)
    Cre = zero(T)
    Cim = zero(T)
    for k ∈ indices((A,B),(3,2))
      Cre += A[1,m,k] * B[1,k,n] - A[2,m,k] * B[2,k,n]
      Cim += A[1,m,k] * B[2,k,n] + A[2,m,k] * B[1,k,n]
    end
    C[1,m,n] = Cre
    C[2,m,n] = Cim
  end
  return Cc
end
function cmatmul_array_v2!(Cc::AbstractMatrix{Complex{T}}, Ac::AbstractMatrix{Complex{T}}, Bc::AbstractMatrix{Complex{T}}) where {T}
  C = reinterpret(Float64, Cc);
  A = reinterpret(Float64, Ac);
  B = reinterpret(reshape, Float64, Bc);
  @turbo vectorize=2 for n ∈ indices((C,B),(2,3)), m ∈ indices((C,A),1)
    Cmn = zero(T)
    for k ∈ indices((A,B),(2,2))
      Amk = A[m,k]
      Aperm = vpermilps177(Amk)
      Cmn = vfmaddsub(Amk, B[1,k,n], vfmaddsub(Aperm, B[2,k,n], Cmn))
    end
    C[m,n] = Cmn
  end
  return Cc
end

function issue209(M, G, J, H, B, ϕ)
    # tmp = similar(ϕ, G-1, (2*J+1)*(H + 1));
    tmp = view(fill(eltype(ϕ)(123456789), G+15, (2*J+1)*(H + 1)+16), 9:G+7, 9:(2*J+1)*(H + 1)+8);
    Bf = reinterpret(reshape, Float64, B);
    ϕf = reinterpret(reshape, Float64, ϕ);
    tmpf = reinterpret(reshape, Float64, tmp);
    jmax = 2*J + 1
    # B is being indexed via ptr offsetting
    # thus B's initial `gesp`ing must set it up for this
    # currently, it isn't because  `jj` and `hh` are loop induct vars
    for mm = 1:M
        m_idx = M + 2 - mm
        @turbo for hh = 1:H+1
            h_idx = (hh - 1)*jmax
            for jj = 1:jmax, gg = 1:G-1
                tmpf[1, gg, jj + h_idx] = ϕf[1, jj, gg+1, hh, m_idx] +
                                          Bf[1, jj, gg,   hh, m_idx]
                tmpf[2, gg, jj + h_idx] = ϕf[2, jj, gg+1, hh, m_idx] +
                                          Bf[2, jj, gg,   hh, m_idx]
            end
        end
    end
    parent(tmp)
end
function issue209_noavx(M, G, J, H, B, ϕ)
    tmp = view(fill(eltype(ϕ)(123456789), G+15, (2*J+1)*(H + 1)+16), 9:G+7, 9:(2*J+1)*(H + 1)+8);
    Bf = reinterpret(reshape, Float64, B);
    ϕf = reinterpret(reshape, Float64, ϕ);
    tmpf = reinterpret(reshape, Float64, tmp);
    jmax = 2*J + 1
    for mm = 1:M
        m_idx = M + 2 - mm
        for hh = 1:H+1
            h_idx = (hh - 1)*jmax
            for jj = 1:jmax, gg = 1:G-1
                tmpf[1, gg, jj + h_idx] = ϕf[1, jj, gg+1, hh, m_idx] +
                                          Bf[1, jj, gg,   hh, m_idx]
                tmpf[2, gg, jj + h_idx] = ϕf[2, jj, gg+1, hh, m_idx] +
                                          Bf[2, jj, gg,   hh, m_idx]
            end
        end
    end
    parent(tmp)
end
using LoopVectorization

function sumdim2_turbo!(r1, r2)
  @turbo thread=true for j = indices((r1,r2),(3,4)), i ∈ indices((r1,r2),(2,3))
    r1[1,i,j] = r2[1,1,i,j] + r2[1,2,i,j]
    r1[2,i,j] = r2[2,1,i,j] - r2[2,2,i,j]
    r1[3,i,j] = r2[3,1,i,j] * r2[3,2,i,j]
    r1[4,i,j] = r2[4,1,i,j] / r2[4,2,i,j]
  end
  r1
end
function sumdim2!(r1, r2)
  @inbounds @fastmath for j = indices((r1,r2),(3,4)), i ∈ indices((r1,r2),(2,3))
    r1[1,i,j] = r2[1,1,i,j] + r2[1,2,i,j]
    r1[2,i,j] = r2[2,1,i,j] - r2[2,2,i,j]
    r1[3,i,j] = r2[3,1,i,j] * r2[3,2,i,j]
    r1[4,i,j] = r2[4,1,i,j] / r2[4,2,i,j]
  end
  r1
end

# Issue 287
 function my_gemm_noturbo!(out, s::Matrix{UInt8}, V)
           Vcols = size(V, 2)
           srows = size(s, 1)
           scols = size(s, 2)
           k = srows >> 2
           rem = srows & 3
           @inbounds @fastmath for c in 1:Vcols
               for j in 1:scols
                   for l in 1:k
                       block = s[l, j]
                       for p in 1:4
                           Aij = (block >> (2 * (p - 1))) & 3
                           out[4*(l - 1) + p, c] += ((Aij >= 2) + (Aij == 3)) * V[j, c]
                       end
                   end
               end
           end
           # TODO handle rem
       end
function my_gemm_unroll(out, s::Matrix{UInt8}, V)
  Vcols = size(V, 2)
  srows = size(s, 1)
  scols = size(s, 2)
  k = srows >> 2
  rem = srows & 3
  @avx for c in 1:Vcols
    for j in 1:scols
      for l in 1:k
        block = s[l, j]
        for p in 1:4
          Aij = (block >> (2 * (p - 1))) & 3
          out[4*(l - 1) + p, c] += ((Aij >= 2) + (Aij == 3)) * V[j, c]
        end
      end
    end
  end
  # TODO handle rem
end
function my_gemm_manual_unroll(out, s::Matrix{UInt8}, V)
  Vcols = size(V, 2)
  srows = size(s, 1)
  scols = size(s, 2)
  k = srows >> 2
  rem = srows & 3
  @avx for c in 1:Vcols
    for j in 1:scols
      for l in 1:k
        block = s[l, j]
        # unrolled loop
        thisiszero = 0
        p = 1
        Aij = (block >> (2 * (p - 1))) & 3
        out[4*(l - 1) + p, c + thisiszero] += ((Aij >= 2) + (Aij == 3)) * V[j, c]
        p = 2
        Aij = (block >> (2 * (p - 1))) & 3
        out[4*(l - 1) + p, c] += ((Aij >= 2) + (Aij == 3)) * V[j, c]
        p = 3
        Aij = (block >> (2 * (p - 1))) & 3
        out[4*(l - 1) + p + thisiszero, c] += ((Aij >= 2) + (Aij == 3)) * V[j, c]
        p = 4
        Aij = (block >> (2 * (p - 1))) & 3
        out[4*(l - 1) + p, c] += ((Aij >= 2) + (Aij == 3)) * V[j, c]
      end
    end
  end
  # TODO handle rem
end
function my_gemm_nexpr_unroll(out, s::Matrix{UInt8}, V)
  Vcols = size(V, 2)
  srows = size(s, 1)
  scols = size(s, 2)
  k = srows >> 2
  rem = srows & 3
  @turbo for c in 1:Vcols
    for j in 1:scols
      for l in 1:k
        block = s[l, j]
        # unrolled loop
        Base.Cartesian.@nexprs 4 p -> begin
          Aij = (block >> (2 * (p - 1))) & 3
          out[4*(l - 1) + p, c] += ((Aij >= 2) + (Aij == 3)) * V[j, c]
        end
      end
    end
  end
  # TODO handle rem
end

function readraw_turbo!(img, raw)
  npack = length(raw) ÷ 3
  @turbo for i in 0:npack-1
    img[1+4i] = raw[2+3i] << 4
    img[2+4i] = raw[1+3i]
    img[3+4i] = raw[2+3i]
    img[4+4i] = raw[3+3i]
  end
  img
end
function readraw!(img, raw)
  npack = length(raw) ÷ 3
  @inbounds @simd for i in 0:npack-1
    img[1+4i] = raw[2+3i] << 4
    img[2+4i] = raw[1+3i]
    img[3+4i] = raw[2+3i]
    img[4+4i] = raw[3+3i]
  end
  img
end

@testset "shuffles load/stores" begin
  @show @__LINE__
  for i ∈ 1:128
    ac = rand(Complex{Float64}, i);
    bc = rand(Complex{Float64}, i);
    dsimd = dot_simd(ac, bc)
    if VERSION ≥ v"1.6.0-rc1"
      @test dsimd ≈ cdot_mat(ac, bc)
    end
    @test dsimd ≈ cdot_affine(ac, bc) ≈ cdot_stride(ac, bc)

    xq = [ntuple(_ -> rand(), Val(4)) for _ ∈ 1:i];
    yq = [ntuple(_ -> rand(), Val(4)) for _ ∈ 1:i];
    xqv = reinterpret(Float64, xq);
    yqv = reinterpret(Float64, yq);
    qsimd = Base.vect(qdot_simd(xq, yq)...);
    if VERSION ≥ v"1.6.0-rc1"
      xqm = reinterpret(reshape, Float64, xq);
      yqm = reinterpret(reshape, Float64, yq);
      @test qsimd ≈ Base.vect(qdot_mat(xqm, yqm)...)
    end
    @test qsimd ≈ Base.vect(qdot_affine(xqv, yqv)...) ≈ Base.vect(qdot_stride(xqv, yqv)...)

    if VERSION ≥ v"1.6.0-rc1"
      for j ∈ max(1,i-5):i+5, k ∈ max(1,i-5,i+5)
        Ac = rand(Complex{Float64}, j, i);
        Bc = rand(Complex{Float64}, i, k);
        Cc1 = Ac*Bc;
        Cc2 = similar(Cc1);
        Cc3 = similar(Cc1);
        @test Cc1 ≈ cmatmul_array!(Cc2, Ac, Bc)
        Cc2 .= NaN;
        @test Cc1 ≈ cmatmul_array_v2!(Cc2, Ac, Bc)
      end
    end
  end
  @show @__LINE__
  if VERSION ≥ v"1.6.0-rc1"
    M = 10
    G = 50
    J = 50
    H = 30

    # B = rand(Complex{Float64}, 2*J+1, G-1, H+1, M+1);
    # ϕ = rand(Complex{Float64}, 2*J+1, G+1, H+1, M+1);
    rbc = let rb = 1.0:((2*J+17) * (G+15) * (H+17) * (M+17)), rbr = reverse(rb)
      Complex{Float64}[rb[i] + im * rbr[i] for i ∈ eachindex(rb)];
    end
    B = view(reshape(rbc, (2*J+17, G+15, H+17, M+17)), 9:2*J+9, 9:G+9, 9:H+9, 9:M+9) .= rand.() .+ rand.().*im;
    ϕ = view(fill(1e5+1e7im, 2*J+17, G+17, H+17, M+17), 9:2*J+9, 9:G+9, 9:H+9, 9:M+9) .= rand.() .+ rand.().*im;
    @test issue209(M, G, J, H, B, ϕ) ≈ issue209_noavx(M, G, J, H, B, ϕ)
  end

  s = Array{Float64}(undef, 4, 128, 128);
  s2 = rand(4, 2, 128, 128);
  @test sumdim2_turbo!(s, s2) ≈ sumdim2!(similar(s), s2)

  # issue 287
  out_test = zeros(100, 10);
  out_test1 = zeros(100, 10);
  s = rand(UInt8, 25, 100);
  V = rand(100, 10);
  my_gemm_noturbo!(out_test, s, V);
  my_gemm_unroll(out_test1, s, V);
  @test out_test ≈ out_test1
  my_gemm_manual_unroll(fill!(out_test1, 0), s, V);
  @test out_test ≈ out_test1
  my_gemm_nexpr_unroll(fill!(out_test1, 0), s, V);
  @test out_test ≈ out_test1

  w = 2048; raw = rand(UInt8, (3w*w) ÷ 4);
  img1 = Matrix{UInt8}(undef, w, w);
  img2 = Matrix{UInt8}(undef, w, w);
  @test readraw!(img1, raw) == readraw_turbo!(img2, raw)
end


text:
# Julia wrapper for header: /usr/include/libavutil/xtea.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_xtea_alloc,
    av_xtea_init,
    av_xtea_le_init,
    av_xtea_crypt,
    av_xtea_le_crypt

function av_xtea_alloc()
    ccall((:av_xtea_alloc, libavutil), Ptr{AVXTEA}, ())
end

function av_xtea_init(ctx, key::NTuple{16, UInt8})
    ccall((:av_xtea_init, libavutil), Cvoid, (Ptr{AVXTEA}, NTuple{16, UInt8}), ctx, key)
end

function av_xtea_le_init(ctx, key::NTuple{16, UInt8})
    ccall((:av_xtea_le_init, libavutil), Cvoid, (Ptr{AVXTEA}, NTuple{16, UInt8}), ctx, key)
end

function av_xtea_crypt(ctx, dst, src, count::Integer, iv, decrypt::Integer)
    ccall((:av_xtea_crypt, libavutil), Cvoid, (Ptr{AVXTEA}, Ptr{UInt8}, Ptr{UInt8}, Cint, Ptr{UInt8}, Cint), ctx, dst, src, count, iv, decrypt)
end

function av_xtea_le_crypt(ctx, dst, src, count::Integer, iv, decrypt::Integer)
    ccall((:av_xtea_le_crypt, libavutil), Cvoid, (Ptr{AVXTEA}, Ptr{UInt8}, Ptr{UInt8}, Cint, Ptr{UInt8}, Cint), ctx, dst, src, count, iv, decrypt)
end


text:
@testset "minimal_inf_test" begin
    @test inf(emptyinterval()) == +Inf
    @test inf(interval(-Inf,+Inf)) == -Inf
    @test inf(interval(1.0,2.0)) == 1.0
    @test inf(interval(-3.0,-2.0)) == -3.0
    @test inf(interval(-Inf,2.0)) == -Inf
    @test inf(interval(-Inf,0.0)) == -Inf
    @test inf(interval(-Inf,-0.0)) == -Inf
    @test inf(interval(-2.0,Inf)) == -2.0
    @test inf(interval(0.0,Inf)) == -0.0
    @test inf(interval(-0.0,Inf)) == -0.0
    @test inf(interval(-0.0,0.0)) == -0.0
    @test inf(interval(0.0,-0.0)) == -0.0
    @test inf(interval(0.0,0.0)) == -0.0
    @test inf(interval(-0.0,-0.0)) == -0.0
end

@testset "minimal_sup_test" begin
    @test sup(emptyinterval()) == -Inf
    @test sup(interval(-Inf,+Inf)) == +Inf
    @test sup(interval(1.0,2.0)) == 2.0
    @test sup(interval(-3.0,-2.0)) == -2.0
    @test sup(interval(-Inf,2.0)) == 2.0
    @test sup(interval(-Inf,0.0)) == 0.0
    @test sup(interval(-Inf,-0.0)) == 0.0
    @test sup(interval(-2.0,Inf)) == Inf
    @test sup(interval(0.0,Inf)) == Inf
    @test sup(interval(-0.0,Inf)) == Inf
    @test sup(interval(-0.0,0.0)) == 0.0
    @test sup(interval(0.0,-0.0)) == 0.0
    @test sup(interval(0.0,0.0)) == 0.0
    @test sup(interval(-0.0,-0.0)) == 0.0
end

@testset "minimal_mid_test" begin
    @test isnan(mid(emptyinterval()))
    @test mid(interval(-Inf,+Inf)) == 0.0
    @test mid(interval(-0x1.FFFFFFFFFFFFFp1023,+0x1.FFFFFFFFFFFFFp1023)) == 0.0
    @test mid(interval(0.0,2.0)) == 1.0
    @test mid(interval(2.0,2.0)) == 2.0
    @test mid(interval(-2.0,2.0)) == 0.0
    @test mid(interval(0.0,Inf)) == 0x1.FFFFFFFFFFFFFp1023
    @test mid(interval(-Inf,1.2)) == -0x1.FFFFFFFFFFFFFp1023
    @test mid(interval(-0x0.0000000000002P-1022,0x0.0000000000001P-1022)) == 0.0
    @test mid(interval(-0x0.0000000000001P-1022,0x0.0000000000002P-1022)) == 0.0
    @test mid(interval(0x1.FFFFFFFFFFFFFP+1022,0x1.FFFFFFFFFFFFFP+1023)) == 0x1.7FFFFFFFFFFFFP+1023
    @test mid(interval(0x0.0000000000001P-1022,0x0.0000000000003P-1022)) == 0x0.0000000000002P-1022
end

@testset "minimal_rad_test" begin
    @test radius(interval(0.0,2.0)) == 1.0
    @test radius(interval(2.0,2.0)) == 0.0
    @test isnan(radius(emptyinterval()))
    @test radius(interval(-Inf,+Inf)) == Inf
    @test radius(interval(0.0,Inf)) == Inf
    @test radius(interval(-Inf, 1.2)) == Inf
    @test radius(interval(-0x0.0000000000002P-1022,0x0.0000000000001P-1022)) == 0x0.0000000000002P-1022
    @test radius(interval(0x0.0000000000001P-1022,0x0.0000000000002P-1022)) == 0x0.0000000000001P-1022
    @test radius(interval(0x1P+0,0x1.0000000000003P+0)) == 0x1P-51
end

@testset "minimal_mid_rad_test" begin
    @test isnan(midpoint_radius(emptyinterval())[1]) && isnan(midpoint_radius(emptyinterval())[2])
    @test midpoint_radius(interval(-Inf,Inf))[1] == 0.0 && midpoint_radius(interval(-Inf,Inf))[2] == Inf
    @test midpoint_radius(interval(-0x1.FFFFFFFFFFFFFP+1023,0x1.FFFFFFFFFFFFFP+1023))[1] == 0.0 && midpoint_radius(interval(-0x1.FFFFFFFFFFFFFP+1023,0x1.FFFFFFFFFFFFFP+1023))[2] == 0x1.FFFFFFFFFFFFFP+1023
    @test midpoint_radius(interval(0.0,2.0))[1] == 1.0 && midpoint_radius(interval(0.0,2.0))[2] == 1.0
    @test midpoint_radius(interval(2.0,2.0))[1] == 2.0 && midpoint_radius(interval(2.0,2.0))[2] == 0.0
    @test midpoint_radius(interval(-2.0,2.0))[1] == 0.0 && midpoint_radius(interval(-2.0,2.0))[2] == 2.0
    @test midpoint_radius(interval(0.0,Inf))[1] == 0x1.FFFFFFFFFFFFFP+1023 && midpoint_radius(interval(0.0,Inf))[2] == Inf
    @test midpoint_radius(interval(-Inf, 1.2))[1] == -0x1.FFFFFFFFFFFFFP+1023 && midpoint_radius(interval(-Inf, 1.2))[2] == Inf
    @test midpoint_radius(interval(-0x0.0000000000002P-1022,0x0.0000000000001P-1022))[1] == 0.0 && midpoint_radius(interval(-0x0.0000000000002P-1022,0x0.0000000000001P-1022))[2] == 0x0.0000000000002P-1022
    @test midpoint_radius(interval(-0x0.0000000000001P-1022,0x0.0000000000002P-1022))[1] == 0.0 && midpoint_radius(interval(-0x0.0000000000001P-1022,0x0.0000000000002P-1022))[2] == 0x0.0000000000002P-1022
    @test midpoint_radius(interval(0x1.FFFFFFFFFFFFFP+1022,0x1.FFFFFFFFFFFFFP+1023))[1] == 0x1.7FFFFFFFFFFFFP+1023 && midpoint_radius(interval(0x1.FFFFFFFFFFFFFP+1022,0x1.FFFFFFFFFFFFFP+1023))[2] == 0x1.0p+1022
    @test midpoint_radius(interval(0x0.0000000000001P-1022,0x0.0000000000003P-1022))[1] == 0x0.0000000000002P-1022 && midpoint_radius(interval(0x0.0000000000001P-1022,0x0.0000000000003P-1022))[2] == 0x0.0000000000001P-1022
end

@testset "minimal_wid_test" begin
    @test diam(interval(2.0,2.0)) == 0.0
    @test diam(interval(1.0,2.0)) == 1.0
    @test diam(interval(1.0,Inf)) == Inf
    @test diam(interval(-Inf,2.0)) == Inf
    @test diam(interval(-Inf,+Inf)) == Inf
    @test isnan(diam(emptyinterval()))
    @test diam(interval(0x1P+0,0x1.0000000000001P+0)) == 0x1P-52
    @test diam(interval(0x1P-1022,0x1.0000000000001P-1022)) == 0x0.0000000000001P-1022
end

@testset "minimal_mag_test" begin
    @test mag(interval(1.0,2.0)) == 2.0
    @test mag(interval(-4.0,2.0)) == 4.0
    @test mag(interval(-Inf,2.0)) == Inf
    @test mag(interval(1.0,Inf)) == Inf
    @test mag(interval(-Inf,+Inf)) == Inf
    @test isnan(mag(emptyinterval()))
    @test mag(interval(-0.0,0.0)) == 0.0
    @test mag(interval(-0.0,-0.0)) == 0.0
end

@testset "minimal_mig_test" begin
    @test mig(interval(1.0,2.0)) == 1.0
    @test mig(interval(-4.0,2.0)) == 0.0
    @test mig(interval(-4.0,-2.0)) == 2.0
    @test mig(interval(-Inf,2.0)) == 0.0
    @test mig(interval(-Inf,-2.0)) == 2.0
    @test mig(interval(-1.0,Inf)) == 0.0
    @test mig(interval(1.0,Inf)) == 1.0
    @test mig(interval(-Inf,+Inf)) == 0.0
    @test isnan(mig(emptyinterval()))
    @test mig(interval(-0.0,0.0)) == 0.0
    @test mig(interval(-0.0,-0.0)) == 0.0
end


text:
cora_init() = register(DataDep(
    "Cora",
    """
    The full Cora citation network dataset from the
    `"Deep Gaussian Embedding of Graphs: Unsupervised Inductive Learning via
    Ranking" <https://arxiv.org/abs/1707.03815>`_ paper.
    Nodes represent documents and edges represent citation links.
    """,
    "https://github.com/abojchevski/graph2gauss/raw/master/data/cora.npz",
    "62e054f93be00a3dedb15b7ac15a2a07168ceab68b40bf95f54d2289d024c6bc";
    post_fetch_method=preprocess_cora,
))

function preprocess_cora(local_path)
    py"""
    import numpy as np
    import scipy.sparse as sp
    data = np.load($local_path, allow_pickle=True)
    A = sp.csr_matrix((data['adj_data'], data['adj_indices'], data['adj_indptr']), shape=data['adj_shape'])
    X = sp.csr_matrix((data['attr_data'], data['attr_indices'], data['attr_indptr']), shape=data['attr_shape'])
    """

    graph = SparseMatrixCSC(Array(py"A.toarray()"))
    X = SparseMatrixCSC(Array(py"X.toarray()"))
    y = py"data['labels']"
    raw = Dict(:graph=>graph, :all_X=>X, :all_y=>y)

    all_X = sparse(X')
    all_y = Matrix{UInt16}(y')
    sg = to_simpledigraph(graph)
    meta = (graph=(num_V=nv(sg), num_E=ne(sg)),
            all=(features_dim=size(all_X), labels_dim=size(all_y))
            )

    graphfile = replace(local_path, "cora.npz"=>"cora.graph.jld2")
    # trainfile = replace(local_path, "cora.npz"=>"cora.train.jld2")
    # testfile = replace(local_path, "cora.npz"=>"cora.test.jld2")
    allfile = replace(local_path, "cora.npz"=>"cora.all.jld2")
    rawfile = replace(local_path, "cora.npz"=>"cora.raw.jld2")
    metadatafile = replace(local_path, "cora.npz"=>"cora.metadata.jld2")

    @save graphfile sg
    # @save trainfile train_X train_y
    # @save testfile test_X test_y
    @save allfile all_X all_y
    @save rawfile raw
    @save metadatafile meta
end

struct Cora <: Dataset
end

function graphdata(cora::Cora)
    file = datadep"Cora/cora.graph.jld2"
    @load file sg
    sg
end

# function traindata(cora::Cora)
#     file = datadep"Cora/cora.train.jld2"
#     @load file graph train_X train_y
#     graph, train_X, train_y
# end

# function testdata(dataset::Symbol)
#     file = datadep"Cora/cora.test.jld2"
#     @load file graph test_X test_y
#     graph, test_X, test_y
# end

function alldata(cora::Cora)
    file = datadep"Cora/cora.all.jld2"
    @load file all_X all_y
    all_X, all_y
end

function rawdata(cora::Cora)
    file = datadep"Cora/cora.raw.jld2"
    @load file raw
    raw
end

function metadata(cora::Cora)
    file = datadep"Cora/cora.metadata.jld2"
    @load file meta
    meta
end


text:
N = 2
A = randn(N,N)
x = randn(N)
b2 = zeros(N)
b = zeros(N)
for k in 1:N
    local_flux_3 = x[k]
    for n in 1:N
        b[n] += A[k, n] * local_flux_3
    end
    b2[k] += b[k]
end
A' * x - b  # correct
A' * x - b2 # incorrect
##

## Figuring out what the different things mean:
if dim == 1
    Np = (N + 1)
    Nfp = 1
    nface = 2
elseif dim == 2
    Np = (N + 1) * (N + 1)
    Nfp = (N + 1)
    nface = 4
elseif dim == 3
    Np = (N + 1) * (N + 1) * (N + 1)
    Nfp = (N + 1) * (N + 1)
    nface = 6
end
# no 19
for f in [1,2,3,4]
    for e in [1]
        for n in collect(1:(N+1)^(dims-1))
        normal_vector = SVector(
            grid.sgeo[grid.n1id, n, f, e],
            grid.sgeo[grid.n2id, n, f, e],
            grid.sgeo[grid.n3id, n, f, e],
        )
        # Get surface mass, volume mass inverse
        sM, vMI = grid.sgeo[grid.sMid, n, f, e], grid.sgeo[_vMI, n, f, e]
        id⁻, id⁺ = grid.vmap⁻[n, f, e], grid.vmap⁺[n, f, e]
        e⁺ = ((id⁺ - 1) ÷ Np) + 1
        e⁻ = ((id⁻ - 1) ÷ Np) + 1
        # not sure I understand vid⁻ and vid⁺
        vid⁻, vid⁺ = ((id⁻ - 1) % Np) + 1, ((id⁺ - 1) % Np) + 1
        local_flux  = normal_vector[1] * (flux[vid⁻, e⁻, 1] + flux[vid⁺, e⁺, 1])/2
        local_flux += normal_vector[2] * (flux[vid⁻, e⁻, 2] + flux[vid⁺, e⁺, 2])/2
        local_flux += normal_vector[3] * (flux[vid⁻, e⁻, 3] + flux[vid⁺, e⁺, 3])/2
        println("----------------------")
        println("For face = ", f, ", element e = ", e, ", gl point = ", n)
        println("The (x,y,z) coordinates are ")
        println((x[id⁻], y[id⁻], z[id⁻]))
        println("The normal vector is ")
        println(normal_vector)
        println("vmap⁻=", id⁻)
        println("vid⁻=", vid⁻)
        println("local flux=", local_flux)
        println("----------------------")
        x[vid⁻, e]
        normal_vector
        end
    end
end

##
v_flux_divergence[grid.vmap⁻]
s_flux_divergence[grid.vmap⁻]
scatter(x[interior], y[interior])
##
Φ = reshape(flux, (N+1, N+1, size(flux)[2], size(flux)[3]))
divΦ = reshape(v_flux_divergence, (N+1, N+1, size(x)[2]))
∮divΦ = reshape(s_flux_divergence, (N+1, N+1, size(x)[2]))
## grid.sgeo[grid.sMid, :, :, :]
Φ = reshape(flux, (N+1, N+1, N+1, size(flux)[2], size(flux)[3]))
divΦ = reshape(v_flux_divergence, (N+1, N+1, N+1, size(x)[2]))
∮divΦ = reshape(s_flux_divergence, (N+1, N+1, N+1, size(x)[2]))

##
Mr = reshape(M, (N+1,N+1,N+1))
Mir = reshape(Mi, (N+1,N+1,N+1))
i = 1
j = 1
# in the single element case
for i in 1:N+1
    for j in 1:N+1
        println("------------------")
        println("for i = ", i , " and j = ", j)
        println("The local flux 3 should be ")
        println(Mr[i,j,:] .* Φ[i,j,:,1,3])
        println("The value should be")
        println(Mir[i,j,:] .* (grid.D[1]' * ( Mr[i,j,:] .* Φ[i,j,:,1,3] )))
        println("The volume term is ")
        println(divΦ[i,j,:])
        println("The surface term is ")
        println(∮divΦ[i,j,:])
        println("------------------")
    end
end
## more details
for i in 1:N+1
    for j in 1:N+1
        for k in 1:N+1
            println("------------------")
            println("for i = ", i , " and j = ", j, " and k = ", k)
            println("The local flux 3 should be ")
            println(Mr[i,j,k] .* Φ[i,j,k,1,3])
            println("The value should be")
            tmp = Mir[i,j,k] .* (grid.D[1]' * ( Mr[i,j,:] .* Φ[i,j,:,1,3] ))
            println(tmp[k])
            println("The volume term is ")
            println(divΦ[i,j,k])
            println("The surface term is ")
            println(∮divΦ[i,j,k])
            println("------------------")
        end
    end
end

##
i = 1
j = 1
tmp = zeros(N+1)
s_D = copy(grid.D[1])
flux = zeros(N+1)
for k in 1:N+1
    for n in 1:N+1
        tmp[n] += Mir[i,j,k] * (s_D[k,n] * ( Mr[i,j,n] * Φ[i,j,k,1,3] ))
    end
    flux[k] += tmp[k]
end
println("correct computation")
println(tmp)
println("current computation")
println(flux)
println("should be")
println(∮divΦ[i,j,:])
#=
println("computation from kernel")
println(divΦ[i,j,:])
=#

##
## Figuring out what the different things mean:
#=
if dim == 1
    Np = (N + 1)
    Nfp = 1
    nface = 2
elseif dim == 2
    Np = (N + 1) * (N + 1)
    Nfp = (N + 1)
    nface = 4
elseif dim == 3
    Np = (N + 1) * (N + 1) * (N + 1)
    Nfp = (N + 1) * (N + 1)
    nface = 6
end
# no 19
for f in [1,2,3,4]
    for e in [1]
        for n in collect(1:(N+1)^(dims-1))
        normal_vector = SVector(
            grid.sgeo[grid.n1id, n, f, e],
            grid.sgeo[grid.n2id, n, f, e],
            grid.sgeo[grid.n3id, n, f, e],
        )
        # Get surface mass, volume mass inverse
        sM, vMI = grid.sgeo[grid.sMid, n, f, e], grid.sgeo[_vMI, n, f, e]
        id⁻, id⁺ = grid.vmap⁻[n, f, e], grid.vmap⁺[n, f, e]
        e⁺ = ((id⁺ - 1) ÷ Np) + 1
        e⁻ = ((id⁻ - 1) ÷ Np) + 1
        # not sure I understand vid⁻ and vid⁺
        vid⁻, vid⁺ = ((id⁻ - 1) % Np) + 1, ((id⁺ - 1) % Np) + 1
        local_flux  = normal_vector[1] * (flux[vid⁻, e⁻, 1] + flux[vid⁺, e⁺, 1])/2
        local_flux += normal_vector[2] * (flux[vid⁻, e⁻, 2] + flux[vid⁺, e⁺, 2])/2
        local_flux += normal_vector[3] * (flux[vid⁻, e⁻, 3] + flux[vid⁺, e⁺, 3])/2
        println("----------------------")
        println("For face = ", f, ", element e = ", e, ", gl point = ", n)
        println("The (x,y,z) coordinates are ")
        println((x[id⁻], y[id⁻], z[id⁻]))
        println("The normal vector is ")
        println(normal_vector)
        println("vmap⁻=", id⁻)
        println("vid⁻=", vid⁻)
        println("local flux=", local_flux)
        println("----------------------")
        x[vid⁻, e]
        normal_vector
        end
    end
end
=#
#=
##
v_flux_divergence[grid.vmap⁻]
s_flux_divergence[grid.vmap⁻]
scatter(x[interior], y[interior])
##
Φ = reshape(flux, (N+1, N+1, size(flux)[2], size(flux)[3]))
divΦ = reshape(v_flux_divergence, (N+1, N+1, size(x)[2]))
∮divΦ = reshape(s_flux_divergence, (N+1, N+1, size(x)[2]))
## grid.sgeo[grid.sMid, :, :, :]
Φ = reshape(flux, (N+1, N+1, N+1, size(flux)[2], size(flux)[3]))
divΦ = reshape(v_flux_divergence, (N+1, N+1, N+1, size(x)[2]))
∮divΦ = reshape(s_flux_divergence, (N+1, N+1, N+1, size(x)[2]))

##
Mr = reshape(M, (N+1,N+1,N+1))
Mir = reshape(Mi, (N+1,N+1,N+1))
i = 1
j = 1
# in the single element case
for i in 1:N+1
    for j in 1:N+1
        println("------------------")
        println("for i = ", i , " and j = ", j)
        println("The local flux 3 should be ")
        println(Mr[i,j,:] .* Φ[i,j,:,1,3])
        println("The value should be")
        println(Mir[i,j,:] .* (grid.D[1]' * ( Mr[i,j,:] .* Φ[i,j,:,1,3] )))
        println("The volume term is ")
        println(divΦ[i,j,:])
        println("The surface term is ")
        println(∮divΦ[i,j,:])
        println("------------------")
    end
end
## more details
for i in 1:N+1
    for j in 1:N+1
        for k in 1:N+1
            println("------------------")
            println("for i = ", i , " and j = ", j, " and k = ", k)
            println("The local flux 3 should be ")
            println(Mr[i,j,k] .* Φ[i,j,k,1,3])
            println("The value should be")
            tmp = Mir[i,j,k] .* (grid.D[1]' * ( Mr[i,j,:] .* Φ[i,j,:,1,3] ))
            println(tmp[k])
            println("The volume term is ")
            println(divΦ[i,j,k])
            println("The surface term is ")
            println(∮divΦ[i,j,k])
            println("------------------")
        end
    end
end

##
# this is another way to check the answer
Φ = reshape(flux, (N+1, N+1, N+1, size(flux)[2], size(flux)[3]))
∂ˣΦ = zeros(N+1, N+1, N+1, size(flux)[2])
metricterm = grid.vgeo[:, 1, :, :]
for j in 1:(N+1)
    for k in 1:(N+1)
        for e in 1:size(flux)[2]
            ∂ˣΦ[:, j, k, e] = grid.D[1] * Φ[:, j, k, e, 1] * metricterm[1]
        end
    end
end
check = reshape(∂ˣΦ, ((N+1)^dim, size(flux)[2]))
check - analytic_flux_divergence
##
# The check
v_flux_divergence - s_flux_divergence
=#

max_stars_repo_name: CliMA/ImperoHooks

max_issues_repo_name: CliMA/ImperoHooks

max_forks_repo_name: CliMA/ImperoHooks


text:
abstract type ReadCounter end
abstract type BiasModel end

function pushto!(t::IntervalTrees.IntervalBTree{K, V, B}, key::I, value) where {K, V, B, I <: AbstractInterval{K}}
    return _push!(t.root, key, value)
end

function _push!(t::IntervalTrees.InternalNode{K, V, B}, key::AbstractInterval{K}, value) where {K, V, B}
    i = IntervalTrees.findidx(t, key)
    if 1 <= length(t) - 1 && key >= t.keys[i]
        return _push!(t.children[i+1], key, value)
    else
        return _push!(t.children[i], key, value)
    end
end

function _push!(t::IntervalTrees.LeafNode{K, V, B}, key::AbstractInterval{K}, value) where {K, V, B}
    i = IntervalTrees.findidx(t, key)
    if 1 <= i <= length(t) &&
        first(t.entries[i]) == first(key) && last(t.entries[i]) == last(key)
        push!( t.entries[i].value, value )
        return true
    end
    return false
end

function pushzero!( collection::Dict{K,V}, key, value ) where {K, V <: ReadCounter}
   if !haskey( collection, key )
      collection[key] = zero(V)
   end
   push!( collection[key], value )
end

function pushzero!( collection::IntervalMap{K,V}, key, value ) where {K <: Integer, V <: ReadCounter}
   if !haskey( collection, (key.first,key.last) )
      collection[key] = zero(V)
   end
   pushto!( collection, key, value )
end

function Base.get( cnt::R ) where R <: ReadCounter
   if !cnt.isadjusted
      error("Cannot get() read count from model thats not already adjusted!")
   else
      return cnt.count
   end
end

function Base.push!( cnt::R, value::F ) where {R <: ReadCounter, F <: AbstractFloat}
   if cnt.isadjusted
      cnt.count += value
   else
      cnt.count += value
      cnt.isadjusted = true
      #error("Cannot push! a float value to PrimerBiasCounter unless isadjusted is true!")
   end
   cnt
end

Base.identity( cnt::R ) where R <: ReadCounter = get( cnt )
Base.identity( cnt::IntervalValue{I,R} ) where R <: ReadCounter = get( cnt.value )

struct DefaultBiasMod <: BiasModel end

mutable struct DefaultCounter <: ReadCounter
   count::Float64
   isadjusted::Bool

   DefaultCounter() = new( 1.0, true )
   DefaultCounter( i::Float64 ) = new( i, true )
end

const DEFAULTCOUNTER_ZERO = DefaultCounter(0.0)
const DEFAULTCOUNTER_ONE  = DefaultCounter(1.0)

Base.zero( ::Type{DefaultCounter} ) = DefaultCounter(0.0) #DEFAULTCOUNTER_ZERO
Base.one( ::Type{DefaultCounter} )  = DefaultCounter(1.0) #DEFAULTCOUNTER_ONE
Base.push!( cnt::DefaultCounter, value::Float64 ) = (cnt.count += value)

primer_adjust!( cnt::DefaultCounter, m::B ) where B <: BiasModel = (cnt.isadjusted = true)
gc_adjust!( cnt::DefaultCounter, m::B ) where B <: BiasModel = (cnt.isadjusted = true)

count!( mod::DefaultBiasMod, _ ) = 1.0
count!( mod::DefaultBiasMod, x, y ) = 1.0

# 5' Sequence bias model
# Implementation of Hansen et al, NAR 2010

# hash by identity for speed
Base.hash( v::UInt16 ) = v

function add_to_back!( back::Vector{Float64}, exp::ExpectedGC, weight::Float64 )
   for i in 1:length(exp)
      back[i] += exp[i] * weight
   end
   mod
end

primer_normalize!( mod::DefaultBiasMod )   = 1.0
gc_normalize!( mod::DefaultBiasMod, x, y ) = 1.0

struct JointBiasTemp
   primer::Vector{UInt16}
   gc::Vector{Int32}
end

struct JointBiasMod <: BiasModel
   fore::Vector{Float64}
   back::Vector{Float64}
   size::Int64
   backoffset::Int64
   backnpos::Int64
   foreoffset::Int64
   forenpos::Int64

   gcfore::Vector{Float64}
   gcback::Vector{Float64}
   gcsize::Int64
   gcwidth::Float64
   gcoffset::Int64
   gcincrement::Int64

   temp::JointBiasTemp

   function JointBiasMod( ksize::Int=6, backoff::Int=23, backn::Int=10, foreoff::Int=1, foren::Int=2,
                          gcoffset::Int=5, gcincrement::Int=10; gc_param = GCBiasParams(50, 0.05) )
      fore = zeros( 4^ksize )
      back = zeros( 4^ksize )
      primer_temp = zeros(UInt16, foren)

      gcfore = zeros(Float64, Int(ceil(1.0 / gc_param.width)+1))
      gcback = zeros(Float64, Int(ceil(1.0 / gc_param.width)+1))
      gctemp = zeros(Int32, 100)

      temp = JointBiasTemp( primer_temp, gctemp )

      new( fore, back, ksize, backoff, backn, foreoff, foren,
           gcfore, gcback, gc_param.length, gc_param.width, gcoffset, gcincrement,
           temp )
   end
end

function primer_normalize!( mod::JointBiasMod )
   foresum = sum(mod.fore)
   backsum = sum(mod.back)
   @inbounds for i in 1:length(mod.fore)
      @fastmath mod.fore[i] = mod.fore[i] / foresum
      @fastmath mod.back[i] = mod.back[i] / backsum
   end
end

function primer_count!( mod::JointBiasMod, seq::BioSequence{A} ) where A <: BioSequences.Alphabet
   # 5 prime bias
   idx = 1
   for i in mod.foreoffset:(mod.foreoffset+mod.forenpos-1)
      hept = kmer_index_trailing(UInt16, seq[i:(i+mod.size-1)])
      mod.fore[hept+1] += 1.0
      mod.temp.primer[idx] = hept
      idx += 1
   end
   for i in mod.backoffset:(mod.backoffset+mod.backnpos-1)
      mod.back[kmer_index_trailing(UInt16, seq[i:(i+mod.size-1)])+1] += 1.0
   end

   return mod.temp
end

function gc_count!( mod::JointBiasMod, seq::BioSequence{A} ) where A <: BioSequences.Alphabet
   # GC fragment bias
   idx = 1
   for i in mod.gcoffset:mod.gcincrement:(length(seq)-mod.gcsize)
      gc = gc_content(seq[i:(i+mod.gcsize-1)])
      bin = Int(div( gc, mod.gcwidth )+1)
      mod.gcfore[bin] += 1.0
      mod.temp.gc[idx] = bin
      idx += 1
   end

   return mod.temp
end

function count!( mod::JointBiasMod, seq::BioSequence{A} ) where A <: BioSequences.Alphabet
   # 5 prime bias
   primer_count!( mod, seq )
   # GC fragment bias
   fill!(mod.temp.gc, zero(Int32))
   gc_count!( mod, seq )

   return mod.temp
end

function count!( mod::JointBiasMod, fwd::BioSequence{A}, rev::BioSequence{A} ) where A <: BioSequences.Alphabet
   # 5 prime bias
   primer_count!( mod, fwd )
   # GC fragment bias
   fill!(mod.temp.gc, zero(Int32))
   gc_count!( mod, fwd )
   gc_count!( mod, rev )

   return mod.temp
end

@fastmath value!( mod::JointBiasMod, hept::UInt16 ) = mod.back[hept+1] / mod.fore[hept+1]
@fastmath value!( mod::JointBiasMod, bin::Int ) = mod.fore[bin] > 0.0 ? mod.back[bin] / mod.fore[bin] : 1.0

mutable struct JointBiasCounter <: ReadCounter
   count::Float64
   map::Dict{UInt16,Int32}
   gc::Vector{Int32}
   isadjusted::Bool

   function JointBiasCounter( gc_param = GCBiasParams(50, 0.05) )
      gc = zeros(Int32, Int(ceil(1.0 / gc_param.width)+1))
      new(0.0, Dict{UInt16,Int32}(), gc, false)
   end
end

Base.zero(::Type{JointBiasCounter}) = JointBiasCounter()

function Base.push!( cnt::JointBiasCounter, temp::JointBiasTemp )
   for i in 1:length(temp.primer)
      push!( cnt, temp.primer[i] )
   end
   for i in 1:length(temp.gc)
      if temp.gc[i] > 0
         push!( cnt, temp.gc[i] )
      else
         return cnt
      end
   end
   cnt
end

function Base.push!( cnt::JointBiasCounter, bin::Int32 )
   cnt.gc[bin] += one(Int32)
   cnt
end

function Base.push!( cnt::JointBiasCounter, hept::UInt16 )
   if haskey( cnt.map, hept )
      cnt.map[hept] += 1
   else
      cnt.map[hept] = 1
   end
   cnt
end

adjust!( cnt::JointBiasCounter, mod::JointBiasMod ) = !cnt.isadjusted ? primer_adjust!( cnt, mod ) : gc_adjust!( cnt, mod )

@inline function primer_adjust!( cnt::JointBiasCounter, mod::JointBiasMod )
   for (k,v) in cnt.map
      cnt.count += value!( mod, k ) * v
   end
   cnt.count /= mod.forenpos
   cnt.isadjusted = true
end

@inline function gc_adjust!( cnt::JointBiasCounter, mod::JointBiasMod )
   weight = 0.0
   for i in 1:length(cnt.gc)
      weight += value!( mod, i ) * cnt.gc[i]
   end
   if sum(cnt.gc) > 0
      weight /= sum(cnt.gc)
      cnt.count *= weight
   end
   cnt
end


text:
module Taskwarrior

import JSON
using Dates

export colnames, tasktable, taskjson

colspecs = Dict{String, DataType}()
for line in readlines(`task columns`)
	result = match(r"([[:alnum:]]+)( +)([[:alnum:]]+)", line)
	if typeof(result) == Nothing
		continue
	elseif result.captures[3] in ["string", "date", "numeric"]
		if result.captures[3] == "string"
			push!(colspecs, result.captures[1] => String)
		elseif result.captures[3] == "date"
			push!(colspecs, result.captures[1] => DateTime)
		else
			push!(colspecs, result.captures[1] => Float64)
		end
	end
end

function colnames()
	collect(keys(colspecs))
end

function taskDateTime(dateString)
	if ismissing(dateString)
		return missing
	end
	year = parse(Int, dateString[1:4])
	month = parse(Int, dateString[5:6])
	day = parse(Int, dateString[7:8])
	hour = parse(Int, dateString[10:11])
	minute = parse(Int, dateString[12:13])
	second = parse(Int, dateString[14:15])
	DateTime(year, month, day, hour, minute, second)
end

function taskjson()
	JSON.parse(read(`task export`, String))
end

"""
    tasktable(x = colnames())

Return a table (as a Dictionary) of Taskwarrior
tasks.

# Arguments
-  'x::Array{String,1}=colnames()' : select which task attributes are returned.

# Examples
tasktable(["uuid", "entry", "description"])
"""
function tasktable(colnames::Array{String, 1} = colnames())
	tasks = taskjson()
	columns = Dict{String, Array}()
	for colname in colnames
		coltype = get(colspecs, colname, String)
		if colname == "tags"
			coltype = Array{String, 1}
		end
		column = Union{coltype, Missing}[]
		for task in tasks
			value = get(task, colname, missing)
			if coltype == DateTime
				value = taskDateTime(value)
			end
			push!(column, value)
		end
		push!(columns, colname => column)
	end
	columns
end

end # module


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

mutable struct P2SVpnConnectionRequest <: SwaggerModel
    vpnConnectionIds::Any # spec type: Union{ Nothing, Vector{String} } # spec name: vpnConnectionIds

    function P2SVpnConnectionRequest(;vpnConnectionIds=nothing)
        o = new()
        validate_property(P2SVpnConnectionRequest, Symbol("vpnConnectionIds"), vpnConnectionIds)
        setfield!(o, Symbol("vpnConnectionIds"), vpnConnectionIds)
        o
    end
end # type P2SVpnConnectionRequest

const _property_map_P2SVpnConnectionRequest = Dict{Symbol,Symbol}(Symbol("vpnConnectionIds")=>Symbol("vpnConnectionIds"))
const _property_types_P2SVpnConnectionRequest = Dict{Symbol,String}(Symbol("vpnConnectionIds")=>"Vector{String}")
Base.propertynames(::Type{ P2SVpnConnectionRequest }) = collect(keys(_property_map_P2SVpnConnectionRequest))
Swagger.property_type(::Type{ P2SVpnConnectionRequest }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_P2SVpnConnectionRequest[name]))}
Swagger.field_name(::Type{ P2SVpnConnectionRequest }, property_name::Symbol) =  _property_map_P2SVpnConnectionRequest[property_name]

function check_required(o::P2SVpnConnectionRequest)
    true
end

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