text:
# * Derivatives
# ** Three-point stencils
α(::FiniteDifferences{T}, ::Integer) where T = one(T)
β(::FiniteDifferences{T}, ::Integer) where T = one(T)
γ(::FiniteDifferences{T}, ::Integer) where T = zero(T)

# Three-point stencil with variationally derived coefficients for
# staggered uniform grid.
α( ::RadialDifferences,    j::Integer)         = j^2/(j^2 - 1/4)
β(B::RadialDifferences{T}, j::Integer) where T = (j^2 - j + 1/2)/(j^2 - j + 1/4) + (j == 1 ? B.δβ₁ : zero(T))
γ( ::RadialDifferences{T},  ::Integer) where T = zero(T)

α(B::StaggeredFiniteDifferences,    j::Integer)         = B.α[j]
β(B::StaggeredFiniteDifferences{T}, j::Integer) where T = B.β[j]
γ( ::StaggeredFiniteDifferences{T},  ::Integer) where T = zero(T)

α(B::StaggeredFiniteDifferences)            = B.α
β(B::StaggeredFiniteDifferences)            = B.β
γ(B::StaggeredFiniteDifferences{T}) where T = zeros(T, length(B.r))

α( ::NumerovFiniteDifferences{T},  ::Integer) where T = one(T)
β(B::NumerovFiniteDifferences{T}, j::Integer) where T = one(T) + (j == 1 ? B.δβ₁ : zero(T))
γ(B::NumerovFiniteDifferences{T}, j::Integer) where T = (j == 1 ? B.λ : zero(T))

α(B::AbstractFiniteDifferences) = α.(Ref(B), B.j[1:end-1])
β(B::AbstractFiniteDifferences) = β.(Ref(B), B.j)
γ(B::AbstractFiniteDifferences) = γ.(Ref(B), B.j)

function α(B̃::RestrictedFiniteDifferences)
    j = indices(B̃,2)
    α(parent(B̃))[j[1]:(j[end]-1)]
end

for f in [:β, :γ]
    @eval begin
        function $f(B̃::RestrictedFiniteDifferences)
            $f(parent(B̃))[indices(B̃,2)]
        end
    end
end

for f in [:α, :β, :γ]
    @eval begin
        function $f(Ã::BasisOrRestricted{<:AbstractFiniteDifferences},
                    B̃::BasisOrRestricted{<:AbstractFiniteDifferences},
                    band=0)
            i,i′ = extrema(indices(Ã,2))
            j,j′ = extrema(indices(B̃,2))
            k = max(i+min(0,band),j-max(0,band)):min(i′+min(0,band),j′-max(0,band))
            $f(parent(B̃))[k]
        end
    end
end

# ** Dispatch aliases

const FlatFirstDerivative{Basis<:AbstractFiniteDifferences} = Mul{<:Any, <:Tuple{
    <:QuasiAdjoint{<:Any, <:Basis},
    <:Derivative,
    <:Basis}}
const LazyFirstDerivative{Basis<:AbstractFiniteDifferences} = Mul{<:Any, <:Tuple{
    <:Mul{<:Any, <:Tuple{
        <:QuasiAdjoint{<:Any, <:Basis},
        <:Derivative}},
    <:Basis}}

const FirstDerivative{Basis} = Union{FlatFirstDerivative{Basis}, LazyFirstDerivative{Basis}}

const FlatSecondDerivative{Basis<:AbstractFiniteDifferences} = Mul{<:Any, <:Tuple{
    <:QuasiAdjoint{<:Any, <:Basis},
    <:QuasiAdjoint{<:Any, <:Derivative},
    <:Derivative,
    <:Basis}}
const LazySecondDerivative{Basis<:AbstractFiniteDifferences} = Mul{<:Any, <:Tuple{
    <:Mul{<:Any, <:Tuple{
        <:Mul{<:Any, <:Tuple{
            <:QuasiAdjoint{<:Any, <:Basis}, <:QuasiAdjoint{<:Any, <:Derivative}}},
        <:Derivative}},
    <:Basis}}

const SecondDerivative{Basis} = Union{FlatSecondDerivative{Basis},LazySecondDerivative{Basis}}

const FirstOrSecondDerivative{Basis} = Union{FirstDerivative{Basis},SecondDerivative{Basis}}

# ** Materialization

function derivative_matrix(::Type{T}, Ac, B, diff_order) where T
    A = parent(Ac)
    parent(A) == parent(B) ||
        throw(ArgumentError("Cannot multiply functions on different grids"))
    Ai,Bi = last(indices(A)), last(indices(B))
    if Ai == Bi
        n = length(Ai)
        if diff_order == 1
            Tridiagonal(Vector{T}(undef, n-1),
                        Vector{T}(undef, n),
                        Vector{T}(undef, n-1))
        elseif diff_order == 2
            SymTridiagonal(Vector{T}(undef, n),
                           Vector{T}(undef, n-1))
        end
    else
        m,n = length(Ai),length(Bi)
        offset = Ai[1]-Bi[1]
        BandedMatrix{T}(undef, (m,n), (1-offset,1+offset))
    end
end

@materialize function *(Ac::AdjointBasisOrRestricted{<:AbstractFiniteDifferences},
                        D::Derivative,
                        B::BasisOrRestricted{<:AbstractFiniteDifferences})
    FiniteDifferencesStyle
    T -> begin
        derivative_matrix(T, Ac, B, 1)
    end
    dest::Tridiagonal{T} -> begin
        A = parent(Ac)
        parent(A) == parent(B) ||
            throw(ArgumentError("Cannot multiply functions on different grids"))

        # Central difference approximation
        μ = 2step(B)
        a = α(B)/μ
        dest.dl .= -a
        dest.d .= γ(B)/μ
        dest.du .= a
    end
    dest::BandedMatrix{T} -> begin
        A = parent(Ac)
        parent(A) == parent(B) ||
            throw(ArgumentError("Cannot multiply functions on different grids"))

        dl,d,du = bandrange(dest)
        # Central difference approximation
        μ = 2step(B)
        dest[Band(dl)] .= -α(A,B,-1)/μ
        dest[Band(d)] .= γ(A,B)/μ
        dest[Band(du)] .= α(A,B,1)/μ
    end
end

@materialize function *(Ac::AdjointBasisOrRestricted{<:AbstractFiniteDifferences},
                        Dc::QuasiAdjoint{<:Any,<:Derivative},
                        D::Derivative,
                        B::BasisOrRestricted{<:AbstractFiniteDifferences})
    FiniteDifferencesStyle
    T -> begin
        derivative_matrix(T, Ac, B, 2)
    end
    dest::SymTridiagonal{T} -> begin
        A = parent(Ac)
        parent(A) == parent(B) ||
            throw(ArgumentError("Cannot multiply functions on different grids"))

        δ² = step(B)^2
        dest.dv .= -2β(B)/δ²
        dest.ev .= α(B)/δ²
    end
    dest::BandedMatrix{T} -> begin
        A = parent(Ac)
        parent(A) == parent(B) ||
            throw(ArgumentError("Cannot multiply functions on different grids"))

        dl,d,du = bandrange(dest)
        δ² = step(B)^2
        dest[Band(dl)] .= α(A,B,-1)/δ²
        dest[Band(d)] .= -2β(A,B)/δ²
        dest[Band(du)] .= α(A,B,1)/δ²
    end
end

# *** Numerov derivatives

function NumerovDerivative(::Type{T}, Δ::Tri, ∇::FirstDerivative) where {T,Tri}
    B = last(∇.args)
    # M₁ = Diagonal{T}(I, size(B,2))

    M₁ = similar(Δ)
    M₁.dl .= 1
    M₁.d .= 4
    M₁.du .= 1

    # Eq. (20ff), Muller (1999), λ = λ′ = √3 - 2, but only if
    # `singular_origin==true`.
    B.j[1] == 1 && (M₁.d[1] += B.λ)

    M₁ /= 6one(T)

    NumerovDerivative(Δ, M₁, one(T))
end

function NumerovDerivative(::Type{T}, Δ::Tri, ∇²::SecondDerivative) where {T,Tri}
    B = last(∇².args)

    M₂ = similar(Δ)
    M₂.dv .= 10
    M₂.ev .= 1

    # Eq. (17), Muller (1999)
    B.j[1] == 1 && (M₂.dv[1] -= 2B.δβ₁)

    M₂ /= -6one(T)
    NumerovDerivative(Δ, M₂, -2one(T))
end

materialize(M::Mul{FiniteDifferencesStyle, <:Tuple{
    <:AdjointBasisOrRestricted{B},
    <:Derivative,
    <:BasisOrRestricted{B}}}) where {T,B<:NumerovFiniteDifferences{T}} =
        NumerovDerivative(T, copyto!(similar(M, T), M), M)

materialize(M::Mul{FiniteDifferencesStyle, <:Tuple{
    <:AdjointBasisOrRestricted{B},
    <:QuasiAdjoint{T,<:Derivative},
    <:Derivative,
    <:BasisOrRestricted{B}}}) where {T,B<:NumerovFiniteDifferences{T}} =
    NumerovDerivative(T, copyto!(similar(M, T), M), M)


text:
function almost_pi(N)
    # Your code goes here!
    return 0
end

ext: jl

max_stars_repo_name: basimr/lovelace

max_issues_repo_name: basimr/project-lovelace

max_forks_repo_name: basimr/project-lovelace


text:
function Covariance()
    return Covariance(())
end

function Covariance(arg0::RealMatrix)
    return Covariance((RealMatrix,), arg0)
end

function Covariance(arg0::RealMatrix, arg1::jboolean)
    return Covariance((RealMatrix, jboolean), arg0, arg1)
end

function Covariance(arg0::Vector{Vector{jdouble}})
    return Covariance((Vector{Vector{jdouble}},), arg0)
end

function Covariance(arg0::Vector{Vector{jdouble}}, arg1::jboolean)
    return Covariance((Vector{Vector{jdouble}}, jboolean), arg0, arg1)
end

function covariance(obj::Covariance, arg0::Vector{jdouble}, arg1::Vector{jdouble})
    return jcall(obj, "covariance", jdouble, (Vector{jdouble}, Vector{jdouble}), arg0, arg1)
end

function covariance(obj::Covariance, arg0::Vector{jdouble}, arg1::Vector{jdouble}, arg2::jboolean)
    return jcall(obj, "covariance", jdouble, (Vector{jdouble}, Vector{jdouble}, jboolean), arg0, arg1, arg2)
end

function equals(obj::Object, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_class(obj::Object)
    return jcall(obj, "getClass", Class, ())
end

function get_covariance_matrix(obj::Covariance)
    return jcall(obj, "getCovarianceMatrix", RealMatrix, ())
end

function get_n(obj::Covariance)
    return jcall(obj, "getN", jint, ())
end

function hash_code(obj::Object)
    return jcall(obj, "hashCode", jint, ())
end

function notify(obj::Object)
    return jcall(obj, "notify", void, ())
end

function notify_all(obj::Object)
    return jcall(obj, "notifyAll", void, ())
end

function to_string(obj::Object)
    return jcall(obj, "toString", JString, ())
end

function wait(obj::Object)
    return jcall(obj, "wait", void, ())
end

function wait(obj::Object, arg0::jlong)
    return jcall(obj, "wait", void, (jlong,), arg0)
end

function wait(obj::Object, arg0::jlong, arg1::jint)
    return jcall(obj, "wait", void, (jlong, jint), arg0, arg1)
end


text:
using Requires
@init @require CUDAnative = "be33ccc6-a3ff-5ff2-a52e-74243cff1e17" begin
    using .CUDAnative
end
using KernelAbstractions.Extras: @unroll

@kernel function knl_fillsendbuf!(
    ::Val{Np},
    ::Val{nvar},
    sendbuf,
    buf,
    vmapsend,
    nvmapsend,
) where {Np, nvar}

    i = @index(Global, Linear)
    @inbounds begin
        e, n = fldmod1(vmapsend[i], Np)
        @unroll for s in 1:nvar
            sendbuf[s, i] = buf[n, s, e]
        end
    end
end

@kernel function knl_transferrecvbuf!(
    ::Val{Np},
    ::Val{nvar},
    buf,
    recvbuf,
    vmaprecv,
    nvmaprecv,
) where {Np, nvar}

    i = @index(Global, Linear)
    @inbounds begin
        e, n = fldmod1(vmaprecv[i], Np)
        @unroll for s in 1:nvar
            buf[n, s, e] = recvbuf[s, i]
        end
    end
end

max_stars_repo_name: yairchn/CLIMA

max_issues_repo_name: yairchn/CLIMA

max_forks_repo_name: yairchn/CLIMA


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

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 el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ dc35a05c-3907-11eb-27ee-35fa47c0d2f3
begin
	import Pkg
	Pkg.activate("Mage")

	using GameOne
	using HTTP
	using JSON
	using Plots
	using Images
	using PlutoUI

	plotly()

	function ingredients(path::String)
		# this is from the Julia source code (evalfile in base/loading.jl)
		# but with the modification that it returns the module instead of the last object
		name = Symbol(basename(path))
		m = Module(name)
		Core.eval(m,
			Expr(:toplevel,
				 :(eval(x) = $(Expr(:core, :eval))($name, x)),
				 :(include(x) = $(Expr(:top, :include))($name, x)),
				 :(include(mapexpr::Function, x) = $(Expr(:top, :include))(mapexpr, $name, x)),
				 :(include($path))))
		m
	end

	function get_card_img(img_uri::String)
		img_resp = HTTP.get(img_uri)
		card_img = img_resp.body |> IOBuffer |> load
	end

	const pd = projectdir()
	JSON_DIR = "$pd/games/MtG/MtG.jl/json"
	DATA_FILES = readdir(JSON_DIR)
	DATA_FILE = DATA_FILES[end]
	CARD_DATA = JSON.parsefile("$JSON_DIR/$DATA_FILE")

	md"""
	## MtG Database

	Import / define common card properties
	"""
end

# ╔═╡ 3ba3f3fa-3eb6-11eb-0521-0736dc692258
CREATURES = [ c for c in CARD_DATA if haskey(c, "type_line") && occursin("Creature", c["type_line"]) ]

# ╔═╡ 03b2e1e0-3ebb-11eb-320e-45a1bb06baaa
length(CREATURES)

# ╔═╡ 5ee2ef20-4fcc-11eb-0f00-47f0c1f708b8
k = "type_line"

# ╔═╡ ff6ad0e4-3eb9-11eb-13a2-8d558b83f605
res = ([ c for c in CREATURES if haskey(c,k) && occursin("bird", lowercase(c[k])) ])

# ╔═╡ 455f9bc8-4fd0-11eb-05e9-f93f7f1b1d55
@bind i Slider(1:length(res), show_value=true)

# ╔═╡ 99da634e-4fc5-11eb-37e8-a7cefb18bbe2
C = res[i]

# ╔═╡ 020e45b2-4fc3-11eb-1830-fdc88af88183
sort([keys(C)...])

# ╔═╡ 964caa22-4fc3-11eb-2793-d9a6fd4a5467
try
	get_card_img(C["image_uris"]["border_crop"])
catch
	f = get_card_img(C["card_faces"][1]["image_uris"]["border_crop"])
	r = get_card_img(C["card_faces"][2]["image_uris"]["border_crop"])
	b = vcat(f,r)
end

# ╔═╡ e26a05fc-4fd0-11eb-0865-e5ed03832806
plot([ c["cmc"] for c in res ], ylabel="CMC", xlabel="index")

# ╔═╡ 2e6b3562-4fd0-11eb-0d53-c7f4ac8e6082
length(res)

# ╔═╡ 61295c8e-4fce-11eb-0f8a-3f56656381c1
CARD_DATA;

# ╔═╡ Cell order:
# ╟─dc35a05c-3907-11eb-27ee-35fa47c0d2f3
# ╠═3ba3f3fa-3eb6-11eb-0521-0736dc692258
# ╠═03b2e1e0-3ebb-11eb-320e-45a1bb06baaa
# ╠═020e45b2-4fc3-11eb-1830-fdc88af88183
# ╟─964caa22-4fc3-11eb-2793-d9a6fd4a5467
# ╠═99da634e-4fc5-11eb-37e8-a7cefb18bbe2
# ╠═455f9bc8-4fd0-11eb-05e9-f93f7f1b1d55
# ╠═5ee2ef20-4fcc-11eb-0f00-47f0c1f708b8
# ╠═ff6ad0e4-3eb9-11eb-13a2-8d558b83f605
# ╠═e26a05fc-4fd0-11eb-0865-e5ed03832806
# ╠═2e6b3562-4fd0-11eb-0d53-c7f4ac8e6082
# ╟─61295c8e-4fce-11eb-0f8a-3f56656381c1


text:
"""
    LogTransform <: Transform

Logarithmically transform the data through: sign(x) * log(|x| + 1).

This allows transformations of all real numbers, not just positive ones.
"""
struct LogTransform <: Transform end
cardinality(::LogTransform) = OneToOne()

_logtransform(x) = sign(x) * log(abs(x) + one(x))
_invlogtransform(x) = sign(x) * (exp(sign(x) * x) - one(x))

function _apply(A::AbstractArray, transform::LogTransform; inverse=false, kwargs...)
    inverse && return _invlogtransform.(A)
    return _logtransform.(A)
end

"""
    InverseHyperbolicSine <: Transform

Logarithmically transform the data through: log(x + √(x² + 1)).

This is the inverse hyperbolic sine.
"""
struct InverseHyperbolicSine <: Transform end
cardinality(::InverseHyperbolicSine) = OneToOne()

function _apply(A::AbstractArray, transform::InverseHyperbolicSine; inverse=false, kwargs...)
    inverse && return sinh.(A)
    return asinh.(A)
end


text:
module Generator

export execute

function execute(args::Dict{Symbol,Any})
  println("Generating...")
end

end # Generator


text:
module Stat

import Gadfly
import StatsBase
import Contour
using Colors
using Compat
using Compose
using DataArrays
using DataStructures
using Distributions
using Hexagons
using Loess
using CoupledFields # It is registered in METADATA.jl
using IndirectArrays

import Gadfly: Scale, Coord, input_aesthetics, output_aesthetics,
               default_scales, isconcrete, setfield!, discretize_make_ia
import KernelDensity
# import Distributions: Uniform, Distribution, qqbuild
import IterTools: chain, distinct
import Compat.Iterators: cycle, product

include("bincount.jl")

# Apply a series of statistics.
#
# Args:
#   stats: Statistics to apply in order.
#   scales: Scales used by the plot.
#   aes: A Aesthetics instance.
#
# Returns:
#   Nothing, modifies aes.
#
function apply_statistics(stats::Vector{Gadfly.StatisticElement},
                          scales::Dict{Symbol, Gadfly.ScaleElement},
                          coord::Gadfly.CoordinateElement,
                          aes::Gadfly.Aesthetics)
    for stat in stats
        apply_statistic(stat, scales, coord, aes)
    end
    nothing
end

struct Nil <: Gadfly.StatisticElement end

const nil = Nil

struct Identity <: Gadfly.StatisticElement end

apply_statistic(stat::Identity,
                scales::Dict{Symbol, Gadfly.ScaleElement},
                coord::Gadfly.CoordinateElement,
                aes::Gadfly.Aesthetics) = nothing

const identity = Identity

# Determine bounds of bars positioned at the given values.
function barminmax(vals, iscontinuous::Bool)
    minvalue, maxvalue = extrema(vals)
    span_type = typeof((maxvalue - minvalue) / 1.0)
    barspan = one(span_type)

    if iscontinuous && length(vals) > 1
        sorted_vals = sort(vals)
        T = typeof(sorted_vals[2] - sorted_vals[1])
        z = zero(T)
        minspan = z
        for i in 2:length(vals)
            span = sorted_vals[i] - sorted_vals[i-1]
            if span > z && (span < minspan || minspan == z)
                minspan = span
            end
        end
        barspan = minspan
    end
    position_type = promote_type(typeof(barspan/2.0), eltype(vals))
    minvals = Array{position_type}(length(vals))
    maxvals = Array{position_type}(length(vals))

    for (i, x) in enumerate(vals)
        minvals[i] = x - barspan/2.0
        maxvals[i] = x + barspan/2.0
    end

    return minvals, maxvals
end

struct RectbinStatistic <: Gadfly.StatisticElement end

const rectbin = RectbinStatistic

input_aesthetics(stat::RectbinStatistic) = [:x, :y]
output_aesthetics(stat::RectbinStatistic) = [:xmin, :xmax, :ymin, :ymax]

function apply_statistic(stat::RectbinStatistic,
                         scales::Dict{Symbol, Gadfly.ScaleElement},
                         coord::Gadfly.CoordinateElement,
                         aes::Gadfly.Aesthetics)
    Gadfly.assert_aesthetics_defined("RectbinStatistic", aes, :x, :y)

    isxcontinuous = haskey(scales, :x) && isa(scales[:x], Scale.ContinuousScale)
    isycontinuous = haskey(scales, :y) && isa(scales[:y], Scale.ContinuousScale)

    xminvals, xmaxvals = barminmax(aes.x, isxcontinuous)
    yminvals, ymaxvals = barminmax(aes.y, isycontinuous)

    aes.xmin = xminvals
    aes.xmax = xmaxvals
    aes.ymin = yminvals
    aes.ymax = ymaxvals

    if !isxcontinuous
        aes.pad_categorical_x = Nullable(false)
    end
    if !isycontinuous
        aes.pad_categorical_y = Nullable(false)
    end
end

struct BarStatistic <: Gadfly.StatisticElement
    position::Symbol # :dodge or :stack
    orientation::Symbol # :horizontal or :vertical
end
BarStatistic(; position=:stack, orientation=:vertical) = BarStatistic(position, orientation)

input_aesthetics(stat::BarStatistic) = stat.orientation == :vertical ? [:x] : [:y]
output_aesthetics(stat::BarStatistic) =
    stat.orientation == :vertical ? [:xmin, :xmax] : [:ymin, :ymax]

default_scales(stat::BarStatistic) = stat.orientation == :vertical ?
        [Gadfly.Scale.y_continuous()] : [Gadfly.Scale.x_continuous()]

const bar = BarStatistic

function apply_statistic(stat::BarStatistic,
                         scales::Dict{Symbol, Gadfly.ScaleElement},
                         coord::Gadfly.CoordinateElement,
                         aes::Gadfly.Aesthetics)
    if stat.orientation == :horizontal
        in(:y, Gadfly.defined_aesthetics(aes)) || return
        var = :y
        othervar = :x
        minvar = :ymin
        maxvar = :ymax
        viewminvar = :xviewmin
        viewmaxvar = :xviewmax
        other_viewminvar = :yviewmin
        other_viewmaxvar = :yviewmax
        labelvar = :x_label
    else
        in(:x, Gadfly.defined_aesthetics(aes)) || return
        var = :x
        othervar = :y
        minvar = :xmin
        maxvar = :xmax
        viewminvar = :yviewmin
        viewmaxvar = :yviewmax
        other_viewminvar = :xviewmin
        other_viewmaxvar = :xviewmax
        labelvar = :y_label
    end

    vals = getfield(aes, var)
    if isempty(vals)
      setfield!(aes, minvar, Float64[1.0])
      setfield!(aes, maxvar, Float64[1.0])
      setfield!(aes, var, Float64[1.0])
      setfield!(aes, othervar, Float64[0.0])
      return
    end

    iscontinuous = haskey(scales, var) && isa(scales[var], Scale.ContinuousScale)

    if getfield(aes, minvar) == nothing || getfield(aes, maxvar) == nothing
        minvals, maxvals = barminmax(vals, iscontinuous)

        setfield!(aes, minvar, minvals)
        setfield!(aes, maxvar, maxvals)
    end

    z = zero(eltype(getfield(aes, othervar)))
    if getfield(aes, viewminvar) == nothing && z < minimum(getfield(aes, othervar))
        setfield!(aes, viewminvar, z)
    elseif getfield(aes, viewmaxvar) == nothing && z > maximum(getfield(aes, othervar))
        setfield!(aes, viewmaxvar, z)
    end

    if stat.position == :stack && aes.color !== nothing
        groups = Dict{Any,Float64}()
        for (x, y) in zip(getfield(aes, othervar), vals)
            Gadfly.isconcrete(x) || continue

            if !haskey(groups, y)
                groups[y] = Float64(x)
            else
                groups[y] += Float64(x)
            end
        end

        viewmin, viewmax = extrema(values(groups))
        aes_viewminvar = getfield(aes, viewminvar)
        if aes_viewminvar === nothing || aes_viewminvar > viewmin
            setfield!(aes, viewminvar, viewmin)
        end
        aes_viewmaxvar = getfield(aes, viewmaxvar)
        if aes_viewmaxvar === nothing || aes_viewmaxvar < viewmax
            setfield!(aes, viewmaxvar, viewmax)
        end
    end

    if !iscontinuous
        if stat.orientation == :horizontal
            aes.pad_categorical_y = Nullable(false)
        else
            aes.pad_categorical_x = Nullable(false)
        end
    end
end

struct HistogramStatistic <: Gadfly.StatisticElement
    minbincount::Int
    maxbincount::Int
    position::Symbol # :dodge or :stack
    orientation::Symbol
    density::Bool
end

function HistogramStatistic(; bincount=nothing,
                              minbincount=3,
                              maxbincount=150,
                              position=:stack,
                              orientation=:vertical,
                              density=false)
    if bincount != nothing
        HistogramStatistic(bincount, bincount, position, orientation, density)
    else
        HistogramStatistic(minbincount, maxbincount, position, orientation, density)
    end
end

input_aesthetics(stat::HistogramStatistic) = stat.orientation == :vertical ? [:x] : [:y]
output_aesthetics(stat::HistogramStatistic) =
    stat.orientation == :vertical ? [:x, :y, :xmin, :xmax] : [:y, :x, :ymin, :ymax]

default_scales(stat::HistogramStatistic) = stat.orientation == :vertical ?
        [Gadfly.Scale.y_continuous()] : [Gadfly.Scale.x_continuous()]

const histogram = HistogramStatistic

function apply_statistic(stat::HistogramStatistic,
                         scales::Dict{Symbol, Gadfly.ScaleElement},
                         coord::Gadfly.CoordinateElement,
                         aes::Gadfly.Aesthetics)
    if stat.orientation == :horizontal
        var = :y
        othervar = :x
        minvar = :ymin
        maxvar = :ymax
        viewminvar = :xviewmin
        viewmaxvar = :xviewmax
        labelvar = :x_label
    else
        var = :x
        othervar = :y
        minvar = :xmin
        maxvar = :xmax
        viewminvar = :yviewmin
        viewmaxvar = :yviewmax
        labelvar = :y_label
    end

    Gadfly.assert_aesthetics_defined("HistogramStatistic", aes, var)

    vals = getfield(aes, var)
    stat.minbincount > stat.maxbincount && error("Histogram minbincount > maxbincount")

    if isempty(getfield(aes, var))
        setfield!(aes, minvar, Float64[1.0])
        setfield!(aes, maxvar, Float64[1.0])
        setfield!(aes, var, Float64[1.0])
        setfield!(aes, othervar, Float64[0.0])
        return
    end

    if haskey(scales, var) && isa(scales[var], Scale.DiscreteScale)
        isdiscrete = true
        x_min, x_max = extrema(vals)
        d = x_max - x_min + 1
        bincounts = zeros(Int, d)
        for x in vals
            bincounts[x - x_min + 1] += 1
        end
        x_min -= 0.5 # adjust the left side of the bar
        binwidth = 1.0
    else
        x_min = Gadfly.concrete_minimum(vals)

        isdiscrete = false
        if estimate_distinct_proportion(vals) <= 0.9
            value_set = sort!(collect(Set(vals[Bool[Gadfly.isconcrete(v) for v in vals]])))
            d, bincounts, x_max = choose_bin_count_1d_discrete(
                        vals, value_set, stat.minbincount, stat.maxbincount)
        else
            d, bincounts, x_max = choose_bin_count_1d(
                        vals, stat.minbincount, stat.maxbincount)
        end

        if stat.density
            x_min = Gadfly.concrete_minimum(vals)
            span = x_max - x_min
            binwidth = span / d
            bincounts = bincounts ./ (sum(bincounts) * binwidth)
        end

        binwidth = (x_max - x_min) / d
    end

    if aes.color === nothing
        T = typeof(x_min + 1*binwidth)
        setfield!(aes, othervar, Array{Float64}(d))
        setfield!(aes, minvar, Array{T}(d))
        setfield!(aes, maxvar, Array{T}(d))
        setfield!(aes, var, Array{T}(d))
        for j in 1:d
            getfield(aes, minvar)[j] = x_min + (j - 1) * binwidth
            getfield(aes, maxvar)[j] = x_min + j * binwidth
            getfield(aes, var)[j] = x_min + (j - 0.5) * binwidth
            getfield(aes, othervar)[j] = bincounts[j]
        end
    else
        groups = Dict()
        for (x, c) in zip(vals, cycle(aes.color))
            Gadfly.isconcrete(x) || continue

            if !haskey(groups, c)
                groups[c] = Float64[x]
            else
                push!(groups[c], x)
            end
        end
        T = typeof(x_min + 1*binwidth)
        setfield!(aes, minvar, Array{T}(d * length(groups)))
        setfield!(aes, maxvar, Array{T}(d * length(groups)))
        setfield!(aes, var, Array{T}(d * length(groups)))

        setfield!(aes, othervar, Array{Float64}(d * length(groups)))
        colors = Array{RGB{Float32}}(d * length(groups))

        x_span = x_max - x_min
        stack_height = zeros(Int, d)
        for (i, (c, xs)) in enumerate(groups)
            fill!(bincounts, 0)
            for x in xs
                Gadfly.isconcrete(x) || continue
                if isdiscrete
                    bincounts[round(Int,x)] += 1
                else
                    bin = max(1, min(d, (ceil(Int, (x - x_min) / binwidth))))
                    bincounts[bin] += 1
                end
            end

            if stat.density
                binwidth = x_span / d
                bincounts = bincounts ./ (sum(bincounts) * binwidth)
            end

            stack_height += bincounts[1:d]

            if isdiscrete
                for j in 1:d
                    idx = (i-1)*d + j
                    getfield(aes, var)[idx] = j
                    getfield(aes, othervar)[idx] = bincounts[j]
                    colors[idx] = c
                end
            else
                for j in 1:d
                    idx = (i-1)*d + j
                    getfield(aes, minvar)[idx] = x_min + (j - 1) * binwidth
                    getfield(aes, maxvar)[idx] = x_min + j * binwidth
                    getfield(aes, var)[idx] = x_min + (j - 0.5) * binwidth
                    getfield(aes, othervar)[idx] = bincounts[j]
                    colors[idx] = c
                end
            end
        end

        if stat.position == :stack
            viewmax = Float64(maximum(stack_height))
            aes_viewmax = getfield(aes, viewmaxvar)
            if aes_viewmax === nothing || aes_viewmax < viewmax
                setfield!(aes, viewmaxvar, viewmax)
            end
        end

        aes.color = discretize_make_ia(colors)
    end

    getfield(aes, viewminvar) === nothing && setfield!(aes, viewminvar, 0.0)

    if haskey(scales, othervar)
        data = Gadfly.Data()
        setfield!(data, othervar, getfield(aes, othervar))
        setfield!(data, viewmaxvar, getfield(aes, viewmaxvar))
        Scale.apply_scale(scales[othervar], [aes], data)

        # See issue #560. Stacked histograms on a non-linear y scale are a strange
        # thing. After some discussion, the least confusing thing is to make the stack
        # partitioned linearly. Here we make that adjustment.
        if stat.position == :stack && aes.color != nothing
            # A little trickery to figure out the scale stack height.
            data = Gadfly.Data()
            setfield!(data, othervar, stack_height)
            scaled_stackheight_aes = Gadfly.Aesthetics()
            Scale.apply_scale(scales[othervar], [scaled_stackheight_aes], data)
            scaled_stackheight = getfield(scaled_stackheight_aes, othervar)

            othervals = getfield(aes, othervar)
            for j in 1:d
                naive_stackheight = 0
                for i in 1:length(groups)
                    idx = (i-1)*d + j
                    naive_stackheight += othervals[idx]
                end

                naive_stackheight == 0 && continue

                for i in 1:length(groups)
                    idx = (i-1)*d + j
                    othervals[idx] = scaled_stackheight[j] * othervals[idx] / naive_stackheight
                end
            end
        end
    else
        setfield!(aes, labelvar, Scale.identity_formatter)
    end
end

struct Density2DStatistic <: Gadfly.StatisticElement
    n::Tuple{Int,Int} # Number of points sampled
    bw::Tuple{Real,Real} # Bandwidth used for the kernel density estimation
    levels::Union{Int,Vector,Function}
end
Density2DStatistic(; n=(256,256), bandwidth=(-Inf,-Inf), levels=15) =
      Density2DStatistic(n, bandwidth, levels)

const density2d = Density2DStatistic

input_aesthetics(stat::Density2DStatistic) = [:x, :y]
output_aesthetics(stat::Density2DStatistic) = [:x, :y, :z]

default_scales(::Density2DStatistic) = [Gadfly.Scale.y_continuous()]

function apply_statistic(stat::Density2DStatistic,
                         scales::Dict{Symbol, Gadfly.ScaleElement},
                         coord::Gadfly.CoordinateElement,
                         aes::Gadfly.Aesthetics)
    Gadfly.assert_aesthetics_defined("Density2DStatistic", aes, :x, :y)

    window = (stat.bw[1] <= 0.0 ? KernelDensity.default_bandwidth(aes.x) : stat.bw[1],
              stat.bw[2] <= 0.0 ? KernelDensity.default_bandwidth(aes.y) : stat.bw[2])
    k = KernelDensity.kde((aes.x,aes.y), bandwidth=window, npoints=stat.n)
    aes.z = k.density
    aes.x = collect(k.x)
    aes.y = collect(k.y)
    apply_statistic(ContourStatistic(levels=stat.levels), scales, coord, aes)
end

struct DensityStatistic <: Gadfly.StatisticElement
    # Number of points sampled
    n::Int
    # Bandwidth used for the kernel density estimation
    bw::Real
end
DensityStatistic(; n=256, bandwidth=-Inf) = DensityStatistic(n, bandwidth)

const density = DensityStatistic

input_aesthetics(stat::DensityStatistic) = [:x]
output_aesthetics(stat::DensityStatistic) = [:x, :y]

default_scales(::DensityStatistic) = [Gadfly.Scale.y_continuous()]

function apply_statistic(stat::DensityStatistic,
                         scales::Dict{Symbol, Gadfly.ScaleElement},
                         coord::Gadfly.CoordinateElement,
                         aes::Gadfly.Aesthetics)
    Gadfly.assert_aesthetics_defined("DensityStatistic", aes, :x)

    if aes.color === nothing
        isa(aes.x[1], Real) || error("Kernel density estimation only works on Real types.")

        x_f64 = collect(Float64, aes.x)

        window = stat.bw <= 0.0 ? KernelDensity.default_bandwidth(x_f64) : stat.bw
        f = KernelDensity.kde(x_f64, bandwidth=window, npoints=stat.n)
        aes.x = collect(Float64, f.x)
        aes.y = f.density
    else
        groups = Dict()
        for (x, c) in zip(aes.x, cycle(aes.color))
            if !haskey(groups, c)
                groups[c] = Float64[x]
            else
                push!(groups[c], x)
            end
        end

        colors = Array{RGB{Float32}}(0)
        aes.x = Array{Float64}(0)
        aes.y = Array{Float64}(0)
        for (c, xs) in groups
            window = stat.bw <= 0.0 ? KernelDensity.default_bandwidth(xs) : stat.bw
            f = KernelDensity.kde(xs, bandwidth=window, npoints=stat.n)
            append!(aes.x, f.x)
            append!(aes.y, f.density)
            for _ in 1:length(f.x)
                push!(colors, c)
            end
        end
        aes.color = discretize_make_ia(colors)
    end
    aes.y_label = Gadfly.Scale.identity_formatter
end

struct Histogram2DStatistic <: Gadfly.StatisticElement
    xminbincount::Int
    xmaxbincount::Int
    yminbincount::Int
    ymaxbincount::Int
end

function Histogram2DStatistic(; xbincount=nothing,
                                xminbincount=3,
                                xmaxbincount=150,
                                ybincount=nothing,
                                yminbincount=3,
                                ymaxbincount=150)
    if xbincount != nothing
        xminbincount = xbincount
        xmaxbincount = xbincount
    end

    if ybincount != nothing
        yminbincount = ybincount
        ymaxbincount = ybincount
    end

    Histogram2DStatistic(xminbincount, xmaxbincount, yminbincount, ymaxbincount)
end

input_aesthetics(stat::Histogram2DStatistic) = [:x, :y]
output_aesthetics(stat::Histogram2DStatistic) = [:xmin, :ymax, :ymin, :ymax, :color]

default_scales(::Histogram2DStatistic, t::Gadfly.Theme=Gadfly.current_theme()) =
    [t.continuous_color_scale]

const histogram2d = Histogram2DStatistic

function apply_statistic(stat::Histogram2DStatistic,
                         scales::Dict{Symbol, Gadfly.ScaleElement},
                         coord::Gadfly.CoordinateElement,
                         aes::Gadfly.Aesthetics)
    Gadfly.assert_aesthetics_defined("Histogram2DStatistic", aes, :x, :y)

    x_min, x_max = Gadfly.concrete_minimum(aes.x), Gadfly.concrete_maximum(aes.x)
    y_min, y_max = Gadfly.concrete_minimum(aes.y), Gadfly.concrete_maximum(aes.y)

    if haskey(scales, :x) && isa(scales[:x], Scale.DiscreteScale)
        x_categorial = true
        xminbincount = x_max - x_min + 1
        xmaxbincount = xminbincount
    else
        x_categorial = false
        xminbincount = stat.xminbincount
        xmaxbincount = stat.xmaxbincount
    end

    if haskey(scales, :y) && isa(scales[:y], Scale.DiscreteScale)
        y_categorial = true
        yminbincount = y_max - y_min + 1
        ymaxbincount = yminbincount
    else
        y_categorial = false
        yminbincount = stat.yminbincount
        ymaxbincount = stat.ymaxbincount
    end

    dy, dx, bincounts = choose_bin_count_2d(aes.x, aes.y,
                                            xminbincount, xmaxbincount,
                                            yminbincount, ymaxbincount)

    wx = x_categorial ? 1 : (x_max - x_min) / dx
    wy = y_categorial ? 1 : (y_max - y_min) / dy

    n = 0
    for cnt in bincounts
        if cnt > 0
            n += 1
        end
    end

    if x_categorial
        aes.x = Array{Int64}(n)
    else
        aes.xmin = Array{Float64}(n)
        aes.xmax = Array{Float64}(n)
    end

    if y_categorial
        aes.y = Array{Int64}(n)
    else
        aes.ymin = Array{Float64}(n)
        aes.ymax = Array{Float64}(n)
    end

    k = 1
    for i in 1:dy, j in 1:dx
        cnt = bincounts[i, j]
        if cnt > 0
            if x_categorial
                aes.x[k] = x_min + (j - 1)
            else
                aes.xmin[k] = x_min + (j - 1) * wx
                aes.xmax[k] = x_min + j * wx
            end

            if y_categorial
                aes.y[k] = y_min + (i - 1)
            else
                aes.ymin[k] = y_min + (i - 1) * wy
                aes.ymax[k] = y_min + i * wy
            end
            k += 1
        end
    end
    @assert k - 1 == n

    haskey(scales, :color) || error("Histogram2DStatistic requires a color scale.")
    color_scale = scales[:color]
    typeof(color_scale) <: Scale.ContinuousColorScale ||
            error("Histogram2DStatistic requires a continuous color scale.")

    aes.color_key_title = "Count"

    data = Gadfly.Data()
    data.color = Array{Int}(n)
    k = 1
    for cnt in transpose(bincounts)
        if cnt > 0
            data.color[k] = cnt
            k += 1
        end
    end

    if x_categorial
        aes.xmin, aes.xmax = barminmax(aes.x, false)
        aes.x = discretize_make_ia(aes.x)
        aes.pad_categorical_x = Nullable(false)
    end

    if y_categorial
        aes.ymin, aes.ymax = barminmax(aes.y, false)
        aes.y = discretize_make_ia(aes.y)
        aes.pad_categorical_y = Nullable(false)
    end

    Scale.apply_scale(color_scale, [aes], data)
    nothing
end

# Find reasonable places to put tick marks and grid lines.
struct TickStatistic <: Gadfly.StatisticElement
    in_vars::Vector{Symbol}
    out_var::AbstractString

    granularity_weight::Float64
    simplicity_weight::Float64
    coverage_weight::Float64
    niceness_weight::Float64

    # fixed ticks, or nothing
    ticks::Union{Symbol, AbstractArray}
end

@deprecate xticks(ticks) xticks(ticks=ticks)

xticks(; ticks=:auto,
         granularity_weight=1/4,
         simplicity_weight=1/6,
         coverage_weight=1/3,
         niceness_weight=1/4) =
    TickStatistic([:x, :xmin, :xmax, :xintercept], "x",
              granularity_weight, simplicity_weight, coverage_weight, niceness_weight, ticks)

@deprecate yticks(ticks) yticks(ticks=ticks)

yticks(; ticks=:auto,
         granularity_weight=1/4,
         simplicity_weight=1/6,
         coverage_weight=1/3,
         niceness_weight=1/4) =
    TickStatistic([:y, :ymin, :ymax, :yintercept, :middle, :lower_hinge, :upper_hinge,
        :lower_fence, :upper_fence], "y",
        granularity_weight, simplicity_weight, coverage_weight, niceness_weight, ticks)

# Apply a tick statistic.

# Args:
#   stat: statistic.
#   aes: aesthetics.

# Returns:
#   nothing

# Modifies:
#   aes
#
function apply_statistic(stat::TickStatistic,
                         scales::Dict{Symbol, Gadfly.ScaleElement},
                         coord::Gadfly.CoordinateElement,
                         aes::Gadfly.Aesthetics)
    isa(stat.ticks, Symbol) && stat.ticks != :auto &&
            error("Invalid value $(stat.ticks) for ticks parameter.")

    isa(coord, Coord.SubplotGrid) &&
            error("TickStatistic cannot be applied to subplot coordinates.")

    # don't clobber existing ticks
    getfield(aes, Symbol(stat.out_var, "tick")) == nothing || return

    in_group_var = Symbol(stat.out_var, "group")
    minval, maxval = nothing, nothing
    in_vals = Any[]
    categorical = (:x in stat.in_vars && Scale.iscategorical(scales, :x)) ||
                  (:y in stat.in_vars && Scale.iscategorical(scales, :y))

    for var in stat.in_vars
        categorical && !in(var,[:x,:y]) && continue
        vals = getfield(aes, var)
        if vals != nothing && eltype(vals) != Function
            if minval == nothing
                minval = first(vals)
            end
            if maxval == nothing
                maxval = first(vals)
            end
            T = promote_type(typeof(minval), typeof(maxval))
            T = promote_type(T, eltype(vals))
            minval = convert(T, minval)
            maxval = convert(T, maxval)

            if stat.out_var == "x"
                dsize = aes.xsize === nothing ? [nothing] : aes.xsize
            elseif stat.out_var == "y"
                dsize = aes.ysize === nothing ? [nothing] : aes.ysize
            else
                dsize = [nothing]
            end

            size = aes.size === nothing ? [nothing] : aes.size

            minval, maxval = apply_statistic_typed(minval, maxval, vals, size, dsize)
            push!(in_vals, vals)
        end
    end

    isempty(in_vals) && return

    in_vals = chain(in_vals...)

    # consider forced tick marks
    if stat.ticks != :auto
        minval = min(minval, minimum(stat.ticks))
        maxval = max(maxval, maximum(stat.ticks))
    end

    # TODO: handle the outliers aesthetic

    n = Gadfly.concrete_length(in_vals)

    # check the x/yviewmin/max pesudo-aesthetics
    if stat.out_var == "x"
        if aes.xviewmin != nothing
            minval = min(minval, aes.xviewmin)
        end
        if aes.xviewmax != nothing
            maxval = max(maxval, aes.xviewmax)
        end
    elseif stat.out_var == "y"
        if aes.yviewmin != nothing
            minval = min(minval, aes.yviewmin)
        end
        if aes.yviewmax != nothing
            maxval = max(maxval, aes.yviewmax)
        end
    end

    # take into account a forced viewport in cartesian coordinates.
    strict_span = false
    if typeof(coord) == Coord.Cartesian
        if stat.out_var == "x"
            if coord.xmin !== nothing
                minval = coord.xmin
                strict_span = true
            end
            if coord.xmax !== nothing
                maxval = coord.xmax
                strict_span = true
            end
        elseif stat.out_var == "y"
            if coord.ymin !== nothing
                minval = coord.ymin
                strict_span = true
            end
            if coord.ymax !== nothing
                maxval = coord.ymax
                strict_span = true
            end
        end
    end

    # all the input values in order.
    if stat.ticks != :auto
        grids = ticks = stat.ticks
        viewmin = minval
        viewmax = maxval
        tickvisible = fill(true, length(ticks))
        tickscale = fill(1.0, length(ticks))
    elseif categorical
        ticks = Set{Int}()
        for val in in_vals
            val>0 && push!(ticks, round(Int, val))
        end
        ticks = Int[t for t in ticks]
        sort!(ticks)
        grids = (ticks .- 0.5)[2:end]
        viewmin = minimum(ticks)
        viewmax = maximum(ticks)
        tickvisible = fill(true, length(ticks))
        tickscale = fill(1.0, length(ticks))
    else
        minval, maxval = promote(minval, maxval)

        ticks, viewmin, viewmax = Gadfly.optimize_ticks(minval, maxval, extend_ticks=true,
                granularity_weight=stat.granularity_weight,
                simplicity_weight=stat.simplicity_weight,
                coverage_weight=stat.coverage_weight,
                niceness_weight=stat.niceness_weight,
                strict_span=strict_span)
        grids = ticks
        multiticks = Gadfly.multilevel_ticks(viewmin - (viewmax - viewmin),
                                             viewmax + (viewmax - viewmin))
        tickcount = length(ticks) + sum([length(ts) for ts in values(multiticks)])
        tickvisible = Array{Bool}(tickcount)
        tickscale = Array{Float64}(tickcount)
        i = 1
        for t in ticks
            tickscale[i] = 1.0
            tickvisible[i] = viewmin <= t <= viewmax
            i += 1
        end

        for (scale, ts) in multiticks, t in ts
            push!(ticks, t)
            tickvisible[i] = false
            tickscale[i] = scale
            i += 1
        end
    end

    # We use the first label function we find for any of the aesthetics. I'm not
    # positive this is the right thing to do, or would would be.
    labeler = getfield(aes, Symbol(stat.out_var, "_label"))

    setfield!(aes, Symbol(stat.out_var, "tick"), ticks)
    setfield!(aes, Symbol(stat.out_var, "grid"), grids)
    setfield!(aes, Symbol(stat.out_var, "tick_label"), labeler)
    setfield!(aes, Symbol(stat.out_var, "tickvisible"), tickvisible)
    setfield!(aes, Symbol(stat.out_var, "tickscale"), tickscale)

    viewmin_var = Symbol(stat.out_var, "viewmin")
    if getfield(aes, viewmin_var) === nothing || getfield(aes, viewmin_var) > viewmin
        setfield!(aes, viewmin_var, viewmin)
    end

    viewmax_var = Symbol(stat.out_var, "viewmax")
    if getfield(aes, viewmax_var) === nothing || getfield(aes, viewmax_var) < viewmax
        setfield!(aes, viewmax_var, viewmax)
    end

    nothing
end

function apply_statistic_typed(minval::T, maxval::T, vals, size, dsize) where T
#     for (val, s, ds) in zip(vals, cycle(size), cycle(dsize))
    lensize  = length(size)
    lendsize = length(dsize)
    for (i, val) in enumerate(vals)
        (!Gadfly.isconcrete(val) || !isfinite(val)) && continue

        s = size[mod1(i, lensize)]
        ds = dsize[mod1(i, lendsize)]

        minval, maxval = minvalmaxval(minval, maxval, convert(T, val), s, ds)
    end
    minval, maxval
end

function apply_statistic_typed(minval::T, maxval::T, vals::DataArray{T}, size, dsize) where T
    lensize  = length(size)
    lendsize = length(dsize)
    for i = 1:length(vals)
        vals.na[i] && continue

        val::T = vals.data[i]
        s = size[mod1(i, lensize)]
        ds = dsize[mod1(i, lendsize)]

        minval, maxval = minvalmaxval(minval, maxval, val, s, ds)
    end
    minval, maxval
end

function minvalmaxval(minval::T, maxval::T, val, s, ds) where T
    if val < minval || !isfinite(minval)
        minval = val
    end

    if val > maxval || !isfinite(maxval)
        maxval = val
    end

    if s != nothing && typeof(s) <: AbstractFloat
        minval = min(minval, val - s)::T
        maxval = max(maxval, val + s)::T
    end

    if ds != nothing
        minval = min(minval, val - ds)::T
        maxval = max(maxval, val + ds)::T
    end

    minval, maxval
end

struct BoxplotStatistic <: Gadfly.StatisticElement
    method::Union{Symbol, Vector}
end

BoxplotStatistic(; method=:tukey) = BoxplotStatistic(method)

input_aesthetics(stat::BoxplotStatistic) = [:x, :y]
output_aesthetics(stat::BoxplotStatistic) =
    [:x, :middle, :lower_hinge, :upper_hinge, :lower_fence, :upper_fence, :outliers]

const boxplot = BoxplotStatistic

function apply_statistic(stat::BoxplotStatistic,
                         scales::Dict{Symbol, Gadfly.ScaleElement},
                         coord::Gadfly.CoordinateElement,
                         aes::Gadfly.Aesthetics)
    if aes.y === nothing
        Gadfly.assert_aesthetics_defined("BoxplotStatistic", aes,
            :x, :lower_hinge, :upper_hinge, :lower_fence, :upper_fence)

        aes_color = aes.color === nothing ? [nothing] : aes.color
        groups = Any[]
        for (x, c) in zip(aes.x, cycle(aes_color))
            push!(groups, (x, c))
        end

        if aes.color !== nothing
            aes.color = discretize_make_ia([c for (x, c) in groups],
                filter(!ismissing, aes.color.values))
        end

        return
    end

    if aes.x === nothing
        aes_x = [1]
        aes.x_label = x -> fill("", length(x))
    else
        aes_x = aes.x
    end
    aes_color = aes.color === nothing ? [nothing] : aes.color

    T = isempty(aes.y) ? eltype(aes.y) : typeof(aes.y[1] / 1)
    groups = DefaultOrderedDict(() -> T[])

    for (x, y, c) in zip(cycle(aes_x), aes.y, cycle(aes_color))
        push!(groups[(x, c)], y)
    end

    if aes.y != nothing
        m = length(groups)
        aes.x = Array{eltype(aes.x)}(m)
        aes.middle = Array{T}(m)
        aes.lower_hinge = Array{T}(m)
        aes.upper_hinge = Array{T}(m)
        aes.lower_fence = Array{T}(m)
        aes.upper_fence = Array{T}(m)
        aes.outliers = Vector{T}[]

        for (i, ((x, c), ys)) in enumerate(groups)
            sort!(ys)

            aes.x[i] = x

            if stat.method == :tukey
                aes.lower_hinge[i], aes.middle[i], aes.upper_hinge[i] =
                        quantile(ys, [0.25, 0.5, 0.75])
                iqr = aes.upper_hinge[i] - aes.lower_hinge[i]

                idx = searchsortedfirst(ys, aes.lower_hinge[i] - 1.5iqr)
                aes.lower_fence[i] = ys[idx]

                idx = searchsortedlast(ys, aes.upper_hinge[i] + 1.5iqr)
                aes.upper_fence[i] = ys[idx]
            elseif isa(stat.method, Vector)
                qs = stat.method
                if length(qs) != 5
                    error("Stat.boxplot requires exactly five quantiles.")
                end

                aes.lower_fence[i], aes.lower_hinge[i], aes.middle[i],
                aes.upper_hinge[i], aes.upper_fence[i] = quantile!(ys, qs)
            else
                error("Invalid method specified for State.boxplot")
            end

            push!(aes.outliers,
                 filter(y -> y < aes.lower_fence[i] || y > aes.upper_fence[i], ys))
        end
    end

    if length(aes.x) > 1 && (haskey(scales, :x) && isa(scales[:x], Scale.ContinuousScale))
        xmin, xmax = minimum(aes.x), maximum(aes.x)
        minspan = minimum([xj - xi for (xi, xj) in zip(aes.x[1:end-1], aes.x[2:end])])

        xviewmin = xmin - minspan / 2
        xviewmax = xmax + minspan / 2

        if aes.xviewmin === nothing || aes.xviewmin > xviewmin
            aes.xviewmin = xviewmin
        end

        if aes.xviewmax === nothing || aes.xviewmax < xviewmax
            aes.xviewmax = xviewmax
        end
    end

    if isa(aes_x, IndirectArray)
        aes.x = discretize_make_ia(aes.x, aes_x.values)
    end

    if aes.color !== nothing
        aes.color = discretize_make_ia(RGB{Float32}[c for (x, c) in keys(groups)],
            aes.color.values)
    end

    nothing
end

struct SmoothStatistic <: Gadfly.StatisticElement
    method::Symbol
    smoothing::Float64
end
SmoothStatistic(; method=:loess, smoothing=0.75) = SmoothStatistic(method, smoothing)

const smooth = SmoothStatistic

input_aesthetics(::SmoothStatistic) = [:x, :y]
output_aesthetics(::SmoothStatistic) = [:x, :y]

function apply_statistic(stat::SmoothStatistic,
                         scales::Dict{Symbol, Gadfly.ScaleElement},
                         coord::Gadfly.CoordinateElement,
                         aes::Gadfly.Aesthetics)

    Gadfly.assert_aesthetics_defined("Stat.smooth", aes, :x, :y)
    Gadfly.assert_aesthetics_equal_length("Stat.smooth", aes, :x, :y)

    stat.method in [:loess,:lm] ||
            error("The only Stat.smooth methods currently supported are loess and lm.")

    max_num_steps = 750
    aes_color = aes.color === nothing ? [nothing] : aes.color

    groups = Dict(c => (eltype(aes.x)[], eltype(aes.y)[]) for c in unique(aes_color))
    for (x, y, c) in zip(aes.x, aes.y, cycle(aes_color))
        push!(groups[c][1], x)
        push!(groups[c][2], y)
    end

    local xs, ys, xsp
    aes.x = eltype(aes.x)[]
    # For aes.y returning a Float is ok if `y` is an Int or a Float
    # There does not seem to be strong demand for other types of `y`
    aes.y = Float64[]
    colors = eltype(aes_color)[]

    for (c, (xv, yv)) in groups
        x_min, x_max = minimum(xv), maximum(xv)
        x_min == x_max && error("Stat.smooth requires more than one distinct x value")
        try
            xs = Float64.( eltype(xv) <: Dates.TimeType ? Dates.value.(xv) : xv )
            ys = Float64.( eltype(yv) <: Dates.TimeType ? Dates.value.(yv) : yv )
        catch e
            error("Stat.loess and Stat.lm require that x and y be bound to arrays of plain numbers.")
        end

        nudge = 1e-5 * (x_max - x_min)

        dx = (x_max-x_min)*(1/max_num_steps)
        # For a Date, dx might be 0 days, so correct
        # For Ints, correct dx
        if isa(xv[1], Date)
            dx = max(dx, Dates.Day(1))
        elseif isa(xv[1], Int)
            dx = ceil(Int, dx)
            nudge = 0
        end

        steps = collect((x_min + nudge):dx:(x_max - nudge))
        xsp = Float64.( eltype(steps) <: Dates.TimeType ? Dates.value.(steps) : steps )
        if stat.method == :loess
            smoothys = Loess.predict(loess(xs, ys, span=stat.smoothing), xsp)
        elseif stat.method == :lm
            lmcoeff = linreg(xs,ys)
            smoothys = lmcoeff[2].*xsp .+ lmcoeff[1]
        end

    # New aes
        append!(aes.x, steps)
        append!(aes.y, smoothys)
        append!(colors, fill(c, length(steps)))
    end

    if !(aes.color===nothing)
        aes.color = discretize_make_ia(colors)
    end
end

struct HexBinStatistic <: Gadfly.StatisticElement
    xbincount::Int
    ybincount::Int
end
HexBinStatistic(; xbincount=50, ybincount=50) = HexBinStatistic(xbincount, ybincount)

const hexbin = HexBinStatistic

function apply_statistic(stat::HexBinStatistic,
                         scales::Dict{Symbol, Gadfly.ScaleElement},
                         coord::Gadfly.CoordinateElement,
                         aes::Gadfly.Aesthetics)
    xmin, xmax = minimum(aes.x), maximum(aes.x)
    ymin, ymax = minimum(aes.y), maximum(aes.y)
    xspan, yspan = xmax - xmin, ymax - ymin

    xsize = xspan / stat.xbincount
    ysize = yspan / stat.ybincount

    counts = Dict{(Tuple{Int, Int}), Int}()
    for (x, y) in zip(aes.x, aes.y)
        h = convert(HexagonOffsetOddR, cube_round(x - xmin + xspan/2,
                                                  y - ymin + yspan/2,
                                                  xsize, ysize))
        idx = (h.q, h.r)
        if !haskey(counts, idx)
            counts[idx] = 1
        else
            counts[idx] += 1
        end
    end

    N = length(counts)
    aes.x = Array{Float64}(N)
    aes.y = Array{Float64}(N)
    data = Gadfly.Data()
    data.color = Array{Int}(N)
    k = 1
    for (idx, cnt) in counts
        x, y = center(HexagonOffsetOddR(idx[1], idx[2]), xsize, ysize,
                      xmin - xspan/2, ymin - yspan/2)
        aes.x[k] = x
        aes.y[k] = y
        data.color[k] = cnt
        k += 1
    end
    aes.xsize = [xsize]
    aes.ysize = [ysize]

    color_scale = scales[:color]
    typeof(color_scale) <: Scale.ContinuousColorScale ||
            error("HexBinGeometry requires a continuous color scale.")

    aes.color_key_title = "Count"

    Scale.apply_scale(color_scale, [aes], data)
end

default_scales(::HexBinStatistic, t::Gadfly.Theme) = [t.continuous_color_scale]

struct StepStatistic <: Gadfly.StatisticElement
    direction::Symbol
end
StepStatistic(; direction=:hv) = StepStatistic(direction)

const step = StepStatistic

input_aesthetics(::StepStatistic) = [:x, :y]
output_aesthetics(::StepStatistic) = [:x, :y]

function apply_statistic(stat::StepStatistic,
                         scales::Dict{Symbol, Gadfly.ScaleElement},
                         coord::Gadfly.CoordinateElement,
                         aes::Gadfly.Aesthetics)
    Gadfly.assert_aesthetics_defined("StepStatistic", aes, :x)
    Gadfly.assert_aesthetics_defined("StepStatistic", aes, :y)
    Gadfly.assert_aesthetics_equal_length("StepStatistic", aes, :x, :y)

    p = sortperm(aes.x, alg=MergeSort)
    permute!(aes.x, p)
    permute!(aes.y, p)
    aes.group != nothing && permute!(aes.group, p)
    aes.color != nothing && permute!(aes.color, p)

    if aes.group != nothing
        Gadfly.assert_aesthetics_equal_length("StepStatistic", aes, :x, :group)
        permute!(aes.x, p)
        permute!(aes.y, p)
        permute!(aes.group, p)
        aes.color != nothing && permute!(aes.color, p)
    end

    if aes.color != nothing
        Gadfly.assert_aesthetics_equal_length("StepStatistic", aes, :x, :color)
        # TODO: use this when we switch to 0.4
        # sortperm!(p, aes.color, alg=MergeSort, lt=Gadfly.color_isless)
        p = sortperm(aes.color, alg=MergeSort, lt=Gadfly.color_isless)
        permute!(aes.x, p)
        permute!(aes.y, p)
        permute!(aes.color, p)
        aes.group != nothing && permute!(aes.group, p)
    end

    x_step = Array{eltype(aes.x)}(0)
    y_step = Array{eltype(aes.y)}(0)
    color_step = aes.color == nothing ? nothing : Array{eltype(aes.color)}(0)
    group_step = aes.group == nothing ? nothing : Array{eltype(aes.group)}(0)

    i = 1
    i_offset = 1
    while true
        u = i_offset + div(i - 1, 2) + (isodd(i) || stat.direction != :hv ? 0 : 1)
        v = i_offset + div(i - 1, 2) + (isodd(i) || stat.direction != :vh ? 0 : 1)

        (u > length(aes.x) || v > length(aes.y)) && break

        if (aes.color != nothing &&
             (aes.color[u] != aes.color[i_offset] || aes.color[v] != aes.color[i_offset])) ||
           (aes.group != nothing &&
             (aes.group[u] != aes.color[i_offset] || aes.color[v] != aes.group[i_offset]))
            i_offset = max(u, v)
            i = 1
        else
            push!(x_step, aes.x[u])
            push!(y_step, aes.y[v])
            aes.color != nothing && push!(color_step, aes.color[i_offset])
            aes.group != nothing && push!(group_step, aes.group[i_offset])
            i += 1
        end
    end

    aes.x = x_step
    aes.y = y_step
    aes.color = color_step
    aes.group = group_step
end

struct FunctionStatistic <: Gadfly.StatisticElement
    # Number of points to evaluate the function at
    num_samples::Int
end
FunctionStatistic(; num_samples=250) = FunctionStatistic(num_samples)

const func = FunctionStatistic

default_scales(::FunctionStatistic) = [Gadfly.Scale.x_continuous(), Gadfly.Scale.y_continuous()]

input_aesthetics(::FunctionStatistic) = [:y, :xmin, :xmax]
output_aesthetics(::FunctionStatistic) = [:x, :y, :group]

function apply_statistic(stat::FunctionStatistic,
                         scales::Dict{Symbol, Gadfly.ScaleElement},
                         coord::Gadfly.CoordinateElement,
                         aes::Gadfly.Aesthetics)
    Gadfly.assert_aesthetics_defined("FunctionStatistic", aes, :y)
    Gadfly.assert_aesthetics_defined("FunctionStatistic", aes, :xmin)
    Gadfly.assert_aesthetics_defined("FunctionStatistic", aes, :xmax)
    Gadfly.assert_aesthetics_equal_length("FunctionStatistic", aes, :xmin, :xmax)

    aes.x = Array{Float64}(length(aes.y) * stat.num_samples)
    ys = Array{Float64}(length(aes.y) * stat.num_samples)

    i = 1
    for (f, xmin, xmax) in zip(aes.y, cycle(aes.xmin), cycle(aes.xmax))
        for x in linspace(xmin, xmax, stat.num_samples)
            aes.x[i] = x
            ys[i] = f(x)
            i += 1
        end
    end

    # color was bound explicitly
    if aes.color != nothing
        func_color = aes.color
        aes.color = Array{eltype(aes.color)}(length(aes.y) * stat.num_samples)
        groups = DataArray(Int, length(aes.y) * stat.num_samples)
        for i in 1:length(aes.y)
            aes.color[1+(i-1)*stat.num_samples:i*stat.num_samples] = func_color[i]
            groups[1+(i-1)*stat.num_samples:i*stat.num_samples] = i
        end
        aes.group = discretize_make_ia(groups)
    elseif length(aes.y) > 1 && haskey(scales, :color)
        data = Gadfly.Data()
        data.color = Array{AbstractString}(length(aes.y) * stat.num_samples)
        groups = DataArray(Int, length(aes.y) * stat.num_samples)
        for i in 1:length(aes.y)
            fname = "f<sub>$(i)</sub>"
            data.color[1+(i-1)*stat.num_samples:i*stat.num_samples] = fname
            groups[1+(i-1)*stat.num_samples:i*stat.num_samples] = i
        end
        Scale.apply_scale(scales[:color], [aes], data)
        aes.group = discretize_make_ia(groups)
    end

    data = Gadfly.Data()
    data.y = ys
    Scale.apply_scale(scales[:y], [aes], data)
end

struct ContourStatistic <: Gadfly.StatisticElement
    levels::Union{Int,Vector,Function}
    samples::Int
end
ContourStatistic(; levels=15, samples=150) = ContourStatistic(levels, samples)

input_aesthetics(::ContourStatistic) = [:z, :xmin, :xmax, :ymin, :ymax]
output_aesthetics(::ContourStatistic) = [:x, :y, :color, :group]

const contour = ContourStatistic

default_scales(::ContourStatistic, t::Gadfly.Theme=Gadfly.current_theme()) =
        [Gadfly.Scale.z_func(), Gadfly.Scale.x_continuous(), Gadfly.Scale.y_continuous(),
            t.continuous_color_scale]

function apply_statistic(stat::ContourStatistic,
                         scales::Dict{Symbol, Gadfly.ScaleElement},
                         coord::Gadfly.CoordinateElement,
                         aes::Gadfly.Aesthetics)
    xs = aes.x === nothing ? nothing : convert(Vector{Float64}, aes.x)
    ys = aes.y === nothing ? nothing : convert(Vector{Float64}, aes.y)

    if typeof(aes.z) <: Function
        if xs == nothing && aes.xmin != nothing && aes.xmax != nothing
            xs = linspace(aes.xmin[1], aes.xmax[1], stat.samples)
        end

        if ys == nothing && aes.ymin != nothing && aes.ymax != nothing
            ys = linspace(aes.ymin[1], aes.ymax[1], stat.samples)
        end

        zs = Float64[aes.z(x, y) for x in xs, y in ys]

    elseif typeof(aes.z) <: Matrix
        zs = convert(Matrix{Float64}, aes.z)

        if xs == nothing
            xs = collect(Float64, 1:size(zs)[1])
        end
        if ys == nothing
            ys = collect(Float64, 1:size(zs)[2])
        end
        size(zs) != (length(xs), length(ys)) &&
                error("Stat.contour requires dimension of z to be length(x) by length(y)")
    else
        error("Stat.contour requires either a matrix or a function")
    end

    levels = Float64[]
    contour_xs = eltype(xs)[]
    contour_ys = eltype(ys)[]

    stat_levels = typeof(stat.levels) <: Function ? stat.levels(zs) : stat.levels

    groups = discretize_make_ia(Int[])
    group = 0
    for level in Contour.levels(Contour.contours(xs, ys, zs, stat_levels))
        for line in Contour.lines(level)
            xc, yc = Contour.coordinates(line)
            append!(contour_xs, xc)
            append!(contour_ys, yc)
            for _ in 1:length(xc)
                push!(groups, group)
                push!(levels, Contour.level(level))
            end
            group += 1
        end
    end

    aes.group = groups
    color_scale = get(scales, :color, Gadfly.Scale.color_continuous_gradient())
    Scale.apply_scale(color_scale, [aes], Gadfly.Data(color=levels))
    Scale.apply_scale(scales[:x], [aes],  Gadfly.Data(x=contour_xs))
    Scale.apply_scale(scales[:y], [aes], Gadfly.Data(y=contour_ys))
end

struct QQStatistic <: Gadfly.StatisticElement end

input_aesthetics(::QQStatistic) = [:x, :y]
output_aesthetics(::QQStatistic) = [:x, :y]

const qq = QQStatistic

default_scales(::QQStatistic) =
        [Gadfly.Scale.x_continuous(), Gadfly.Scale.y_continuous]

function apply_statistic(stat::QQStatistic,
                         scales::Dict{Symbol, Gadfly.ScaleElement},
                         coord::Gadfly.CoordinateElement,
                         aes::Gadfly.Aesthetics)

    Gadfly.assert_aesthetics_defined("Stat.qq", aes, :x, :y)
    Gadfly.assert_aesthetics_undefined("State.qq", aes, :color)

    # NOTES:
    #
    # apply_scales happens before apply_statistics, so we need to handle in
    # apply_scales the Distributions that might be bound to x and y... By
    # analogy with Stat.func, we can add a check in apply_statistic which defers
    # application.  Stat.func though requires an ARRAY of Functions, and doesn't
    # work on naked functions bound to aes.y.  If we want to bind Distributions,
    # we'd need to extend the types that are allowed for aes.y/.x (e.g. change
    # type of Aesthetics fields x and y).  Right now these are of type
    # NumericalOrCategoricalAesthetic.  The .x and .y fields are the _only_
    # place where this type is used, but I'm not sure if there's a reason that
    # changing this typealias would be a bad idea...for now I've just used a
    # direct `@compat(Union{NumericalOrCategoricalAesthetic, Distribution})`.
    #
    # TODO:
    #
    # Grouping by color etc.?

    # a little helper function to convert either numeric or distribution
    # variables to a format suitable to input to qqbuild.
    toVecOrDist = v -> typeof(v) <: Distribution ? v : convert(Vector{Float64}, v)

    # check and convert :x and :y to proper types for input to qqbuild
    local xs, ys
    try
        (xs, ys) = map(toVecOrDist, (aes.x, aes.y))
    catch e
        error("Stat.qq requires that x and y be bound to either a Distribution or to arrays of plain numbers.")
    end

    qqq = qqbuild(xs, ys)

    aes.x = qqq.qx
    aes.y = qqq.qy

    # apply_scale to Distribution-bound aesthetics is deferred, so re-apply here
    # (but only for Distribution, numeric data is already scaled).  Only one of
    # :x or :y can be a Distribution since qqbuild will throw an error for two
    # Distributions.
    data = Gadfly.Data()
    if typeof(xs) <: Distribution
        data.x = aes.x
        Scale.apply_scale(scales[:x], [aes], data)
    elseif typeof(ys) <: Distribution
        data.y = aes.y
        Scale.apply_scale(scales[:y], [aes], data)
    end
end

struct ViolinStatistic <: Gadfly.StatisticElement
    # Number of points sampled
    n::Int
end
ViolinStatistic() = ViolinStatistic(300)

input_aesthetics(::ViolinStatistic) = [:x, :y, :width]

const violin = ViolinStatistic

function apply_statistic(stat::ViolinStatistic,
                         scales::Dict{Symbol, Gadfly.ScaleElement},
                         coord::Gadfly.CoordinateElement,
                         aes::Gadfly.Aesthetics)

    isa(aes.y[1], Real) || error("Kernel density estimation only works on Real types.")

    grouped_y = Dict(1=>aes.y)
    grouped_color = Dict{Int, Gadfly.ColorOrNothing}(1=>nothing)
    ux = unique(aes.x)
    uxflag = length(ux) < length(aes.x)
    colorflag = aes.color != nothing

    uxflag && (grouped_y = Dict(x=>aes.y[aes.x.==x] for x in ux))

    grouped_color = (colorflag ? Dict(x=>first(aes.color[aes.x.==x]) for x in ux) :
        uxflag && Dict(x=>nothing for x in ux) )

    aes.x     = Array{Float64}(0)
    aes.y     = Array{Float64}(0)
    aes.width = Array{Float64}(0)
    colors = eltype(aes.color)[]

    for (x, ys) in grouped_y
        window = stat.n > 1 ? KernelDensity.default_bandwidth(ys) : 0.1
        f = KernelDensity.kde(ys, bandwidth=window, npoints=stat.n)
        append!(aes.x, fill(x, length(f.x)))
        append!(aes.y, f.x)
        append!(aes.width, f.density)
        append!(colors, fill(grouped_color[x], length(f.x)))
    end

    colorflag  && (aes.color = colors)

    pad = 0.1
    maxwidth = maximum(aes.width)
    broadcast!(*, aes.width, aes.width, 1 - pad)
    broadcast!(/, aes.width, aes.width, maxwidth)
end

struct JitterStatistic <: Gadfly.StatisticElement
    vars::Vector{Symbol}
    range::Float64
    seed::UInt32
end
JitterStatistic(vars; range=0.8, seed=0x0af5a1f7) = JitterStatistic(vars, range, seed)

x_jitter(; range=0.8, seed=0x0af5a1f7) = JitterStatistic([:x], range=range, seed=seed)
y_jitter(; range=0.8, seed=0x0af5a1f7) = JitterStatistic([:y], range=range, seed=seed)

input_aesthetics(stat::JitterStatistic) = stat.vars
output_aesthetics(stat::JitterStatistic) = stat.vars

function minimum_span(vars::Vector{Symbol}, aes::Gadfly.Aesthetics)
    span = nothing
    for var in vars
        data = getfield(aes, var)
        length(data) < 2 && continue
        dataspan = data[2] - data[1]
        T = eltype(data)
        z = convert(T, zero(T))
        sorteddata = sort(data)
        for δ in diff(sorteddata)
            if δ != z && (δ < dataspan || dataspan == z)
                dataspan = δ
            end
        end

        if span == nothing || (dataspan != nothing && dataspan < span)
            span = dataspan
        end
    end

    return span
end

function apply_statistic(stat::JitterStatistic,
                         scales::Dict{Symbol, Gadfly.ScaleElement},
                         coord::Gadfly.CoordinateElement,
                         aes::Gadfly.Aesthetics)
    span = minimum_span(stat.vars, aes)
    span == nothing && return

    rng = MersenneTwister(stat.seed)
    for var in stat.vars
        data = getfield(aes, var)
        outdata = Array{Float64}(size(data))
        broadcast!(+, outdata, data, stat.range * (rand(rng, length(data)) - 0.5) .* span)
        setfield!(aes, var, outdata)
    end
end

# Bin mean returns the mean of x and y in n bins of x
struct BinMeanStatistic <: Gadfly.StatisticElement
    n::Int
end
BinMeanStatistic(; n=20) = BinMeanStatistic(n)

const binmean = BinMeanStatistic

input_aesthetics(::BinMeanStatistic) = [:x, :y]
output_aesthetics(::BinMeanStatistic) = [:x, :y]

function apply_statistic(stat::BinMeanStatistic,
                         scales::Dict{Symbol, Gadfly.ScaleElement},
                         coord::Gadfly.CoordinateElement,
                         aes::Gadfly.Aesthetics)

    Gadfly.assert_aesthetics_defined("Stat.binmean", aes, :x, :y)

    breaks = quantile(aes.x, collect(1:stat.n)/stat.n)

    Tx = eltype(aes.x)
    Ty = eltype(aes.y)

    if aes.color === nothing
        (aes.x, aes.y) = mean_by_group(aes.x, aes.y, breaks)
    else
        groups = Dict()
        for (x, y, c) in zip(aes.x, aes.y, cycle(aes.color))
            if !haskey(groups, c)
                xs = append!(Tx[], collect(Tx, aes.x))
                ys = append!(Ty[], collect(Ty, aes.y))
                groups[c] = Array[xs, ys]
            else
                push!(groups[c][1], x)
                push!(groups[c][2], y)
            end
        end
        colors = Array{RGB{Float32}}(0)
        aes.x = Array{Tx}(0)
        aes.y = Array{Ty}(0)
        for (c, v) in groups
            (fx, fy) = mean_by_group(v[1], v[2], breaks)
            append!(aes.x, fx)
            append!(aes.y, fy)
            for _ in 1:length(fx)
                push!(colors, c)
            end
        end
        aes.color = discretize_make_ia(colors)
    end
end

function mean_by_group(x::Vector{Tx}, y::Vector{Ty}, breaks::Vector{Float64}) where {Tx, Ty}
    count = zeros(Int64, length(breaks))
    totalx = zeros(Tx, length(breaks))
    totaly = zeros(Ty, length(breaks))
    for i in 1:length(x)
        refs = searchsortedfirst(breaks, x[i])
        count[refs] += 1
        totalx[refs] += x[i]
        totaly[refs] += y[i]
    end
    subset = count .> 0
    count = count[subset]
    return (totalx[subset] ./ count, totaly[subset] ./ count)
end

struct EnumerateStatistic <: Gadfly.StatisticElement
    var::Symbol
end

input_aesthetics(stat::EnumerateStatistic) = [stat.var]
output_aesthetics(stat::EnumerateStatistic) = [stat.var]

function default_scales(stat::EnumerateStatistic)
    if stat.var == :y
        return [Gadfly.Scale.y_continuous()]
    elseif stat.var == :x
        return [Gadfly.Scale.x_continuous()]
    else
        return Gadfly.ScaleElement[]
    end
end

const x_enumerate = EnumerateStatistic(:x)
const y_enumerate = EnumerateStatistic(:y)

function apply_statistic(stat::EnumerateStatistic,
                         scales::Dict{Symbol, Gadfly.ScaleElement},
                         coord::Gadfly.CoordinateElement,
                         aes::Gadfly.Aesthetics)
    has_x = aes.x != nothing
    has_y = aes.y != nothing

    if stat.var == :x && !has_x && has_y
        aes.x = collect(1:length(aes.y))
    elseif stat.var == :y && !has_y && has_x
        aes.y = collect(1:length(aes.x))
    end
end

### Vector Field Statistic

struct VecFieldStatistic <: Gadfly.StatisticElement
    smoothness::Float64
    scale::Float64
    samples::Int64
end
VecFieldStatistic(; smoothness=1.0, scale=1.0, samples=20) =
        VecFieldStatistic(smoothness, scale, samples)

input_aesthetics(stat::VecFieldStatistic) = [:z, :x, :y, :color]
output_aesthetics(stat::VecFieldStatistic) = [:x, :y, :xend, :yend, :color]
default_scales(stat::VecFieldStatistic, t::Gadfly.Theme=Gadfly.current_theme()) =
        [Gadfly.Scale.z_func(), Gadfly.Scale.x_continuous(), Gadfly.Scale.y_continuous(),
            t.continuous_color_scale ]

const vectorfield = VecFieldStatistic

function apply_statistic(stat::VecFieldStatistic,
                         scales::Dict{Symbol, Gadfly.ScaleElement},
                         coord::Gadfly.CoordinateElement,
                         aes::Gadfly.Aesthetics)
    xs = aes.x === nothing ? nothing : Float64.(aes.x)
    ys = aes.y === nothing ? nothing : Float64.(aes.y)

    if isa(aes.z, Function)
        if xs == nothing && aes.xmin != nothing && aes.xmax != nothing
            xs = linspace(aes.xmin[1], aes.xmax[1], stat.samples)
        end
        if ys == nothing && aes.ymin != nothing && aes.ymax != nothing
            ys = linspace(aes.ymin[1], aes.ymax[1], stat.samples)
        end

        zs = Float64[aes.z(x, y) for x in xs, y in ys]

    elseif isa(aes.z, Matrix)
        zs = Float64.(aes.z)

        if xs == nothing
            xs = collect(Float64, 1:size(zs)[1])
        end
        if ys == nothing
            ys = collect(Float64, 1:size(zs)[2])
        end
        if size(zs) != (length(xs), length(ys))
            error("Stat.vectorfield requires dimension of z to be length(x) by length(y)")
        end
    else
        error("Stat.vectorfield requires either a matrix or a function")
    end

    X = vcat([[x y] for x in xs, y in ys]...)
    Z = vec(zs)
    # The next two lines make use of the package CoupledFields.jl
    kpars = GaussianKP(X)
    ∇g = hcat(gradvecfield([stat.smoothness -7.0], X, Z[:,1:1], kpars)...)'
    vecf = [X-∇g*stat.scale X+∇g*stat.scale]

    aes.z = nothing
    aes.x = vecf[:,1]
    aes.y = vecf[:,2]
    aes.xend = vecf[:,3]
    aes.yend = vecf[:,4]
    color_scale = get(scales, :color, Gadfly.Scale.color_continuous_gradient())
    Scale.apply_scale(color_scale, [aes], Gadfly.Data(color=Z))

end

### Hair Statistic

struct HairStatistic <: Gadfly.StatisticElement
    intercept
    orientation::Symbol # :horizontal or :vertical like BarStatistic
end
HairStatistic(;intercept=0.0, orientation=:vertical) = HairStatistic(intercept, orientation)

input_aesthetics(stat::HairStatistic) = [:x, :y]

output_aesthetics(stat::HairStatistic) = [:x, :y, :xend, :yend]

default_scales(stat::HairStatistic) = [Gadfly.Scale.x_continuous(), Gadfly.Scale.y_continuous()]

const hair = HairStatistic

function apply_statistic(stat::HairStatistic,
                         scales::Dict{Symbol, Gadfly.ScaleElement},
                         coord::Gadfly.CoordinateElement,
                         aes::Gadfly.Aesthetics)
    if stat.orientation == :vertical
        aes.xend = aes.x
        aes.yend = fill(stat.intercept, length(aes.y))
    else
        aes.yend = aes.y
        aes.xend = fill(stat.intercept, length(aes.x))
    end
end

### Ellipse Statistic

struct EllipseStatistic <: Gadfly.StatisticElement
    distribution::Type{<:ContinuousMultivariateDistribution}
    levels::Vector{<:AbstractFloat}
    nsegments::Int
end

function EllipseStatistic(;
        distribution::(Type{<:ContinuousMultivariateDistribution})=MvNormal,
        levels::Vector{Float64}=[0.95],
        nsegments::Int=51 )
    return EllipseStatistic(distribution, levels, nsegments)
end

Gadfly.input_aesthetics(stat::EllipseStatistic) = [:x, :y]
Gadfly.output_aesthetics(stat::EllipseStatistic) = [:x, :y]
Gadfly.default_scales(stat::EllipseStatistic) = [Gadfly.Scale.x_continuous(), Gadfly.Scale.y_continuous()]

const ellipse = EllipseStatistic

function Gadfly.Stat.apply_statistic(stat::EllipseStatistic,
                         scales::Dict{Symbol, Gadfly.ScaleElement},
                         coord::Gadfly.CoordinateElement,
                         aes::Gadfly.Aesthetics)

    Dat = [aes.x aes.y]
    grouped_xy = Dict(1=>Dat)
    grouped_color = Dict{Int, Gadfly.ColorOrNothing}(1=>nothing)
    colorflag = aes.color != nothing
    aes.group = (colorflag ? aes.color : aes.group)

    if aes.group != nothing
        ug = unique(aes.group)
        grouped_xy = Dict(g=>Dat[aes.group.==g,:] for g in ug)
        grouped_color = Dict(g=>first(aes.group[aes.group.==g]) for g in ug)
    end

    levels = Float64[]
    colors = eltype(aes.color)[]
    ellipse_x = eltype(Dat)[]
    ellipse_y = eltype(Dat)[]

    dfn = 2
    θ = 2π*(0:stat.nsegments)/stat.nsegments
    n = length(θ)
    for (g, data) in grouped_xy
        dfd = size(data,1)-1
        dhat = fit(stat.distribution, data')
        Σ½ = chol(cov(dhat))
        rv = sqrt.(dfn*[quantile(FDist(dfn,dfd), p) for p in stat.levels])
        ellxy =  [cos.(θ) sin.(θ)] * Σ½
        μ = mean(dhat)
        for r in rv
            append!(ellipse_x, r*ellxy[:,1].+μ[1])
            append!(ellipse_y, r*ellxy[:,2].+μ[2])
            append!(colors, fill(grouped_color[g], n))
            append!(levels, fill(r, n))
        end
    end

    aes.group = discretize_make_ia(levels)
    colorflag && (aes.color = colors)
    aes.x = ellipse_x
    aes.y = ellipse_y
end

end # module Stat


text:
using DataStructures
import GlobalSensitivityAnalysis

mutable struct SobolData <: AbstractSimulationData
    calc_second_order::Bool

    function SobolData(;calc_second_order = false)
        return new(calc_second_order)
    end
end

function Base.show(data::SobolData)
    println("Calc 2nd order: $(data.calc_second_order)")
end

const SobolSimulationDef = SimulationDef{SobolData}
const SobolSimulationInstance = SimulationInstance{SobolData}

function _compute_N(sim_inst::SobolSimulationInstance)
    num_rvs = length(sim_inst.sim_def.rvdict)
    factor = (sim_inst.sim_def.data.calc_second_order ? 2 : 1)
    N = sim_inst.trials / (factor * num_rvs + 2)
    return N
end

function _compute_trials(sim_inst::SobolSimulationInstance, N::Int)
    num_rvs = length(sim_inst.sim_def.rvdict)
    factor = (sim_inst.sim_def.data.calc_second_order ? 2 : 1)
    sim_inst.trials = N * (factor * num_rvs + 2)
end

# Use original distribution when resampling from SampleStore
_get_dist(rv::RandomVariable) = (rv.dist isa SampleStore ? rv.dist.dist : rv.dist)

function sample!(sim_inst::SobolSimulationInstance, samplesize::Int)
    rvdict = sim_inst.sim_def.rvdict
    sim_inst.trials = _compute_trials(sim_inst, samplesize)

    # get the samples to plug in to trials
    payload = create_GSA_payload(sim_inst)
    samples = GlobalSensitivityAnalysis.sample(payload)

    for (i, rv) in enumerate(values(rvdict))
        # use underlying distribution, if known
        orig_dist = _get_dist(rv)

        name = rv.name
        values = samples[:, i]
        rvdict[name] = RandomVariable(name, SampleStore(values, orig_dist))
    end
end

function analyze(sim_inst::SobolSimulationInstance, model_output::AbstractArray{<:Number, N}) where N

    if sim_inst.trials == 0
        error("Cannot analyze simulation with 0 trials.")
    end

    payload = create_GSA_payload(sim_inst)
    return GlobalSensitivityAnalysis.analyze(payload, model_output)
end

function create_GSA_payload(sim_inst::SobolSimulationInstance)
    rv_info = OrderedDict{Symbol, Any}([name => _get_dist(rv) for (name, rv) in sim_inst.sim_def.rvdict])

    # back out N
    N = _compute_N(sim_inst)

    # create payload
    return GlobalSensitivityAnalysis.SobolData(params = rv_info, calc_second_order = sim_inst.sim_def.data.calc_second_order, N = N)
end


text:
immutable Alpha <: ContinuousUnivariateDistribution
    location::Float64
    scale::Float64
    function Alpha(l::Real, s::Real)
        if s < 0.0
        	error("scale must be non-negative")
        else
        	new(float64(l), float64(s))
        end
    end
end

Alpha(location::Real) = Alpha(location, 1.0)
Alpha() = Alpha(0.0, 1.0)

const Levy = Alpha


text:
function W=BoundaryWOutput(U,CG,Param)
OP=CG.OrdPoly+1;
NF=Param.Grid.NumFaces;

v1CG=reshape(U(reshape(CG.Glob,OP*OP*NF,1),1,Param.uPos)...
  ,OP,OP,NF);
v2CG=reshape(U(reshape(CG.Glob,OP*OP*NF,1),1,Param.vPos)...
  ,OP,OP,NF,1);

wCG=...
  -(Param.dXdxIC(:,:,:,1,3,1).*v1CG(:,:,:,1)...
  +Param.dXdxIC(:,:,:,1,3,2).*v2CG(:,:,:,1))...
  ./Param.dXdxIC(:,:,:,1,3,3).*Param.JC(:,:,:,1);
W=zeros(size(U,1),1);
for iM=1:size(CG.FaceGlob,1)
  W(reshape(CG.Glob(:,CG.FaceGlob(iM).Ind,:)...
    ,OP*OP*size(CG.FaceGlob(iM).Ind,1),1),:)...
    =W(reshape(CG.Glob(:,CG.FaceGlob(iM).Ind,:)...
    ,OP*OP*size(CG.FaceGlob(iM).Ind,1),1),:)...
    +reshape(wCG(:,:,CG.FaceGlob(iM).Ind,:)...
    ,OP*OP*size(CG.FaceGlob(iM).Ind,1),1);
end
W=W./CG.M(:,1);
end


text:
mutable struct ECOSls <: ECOSWrapper
    ecos_ptr::Ptr{ECOS.Cpwork}
    A1::AbstractMatrix{Float64}
    b1::Vector{Float64}
    A2::AbstractMatrix{Float64}
    b2::Vector{Float64}
    n_ecos::Int
    n::Int
    m1::Int
    m2::Int
end

mutable struct ECOSlssplx <: ECOSWrapper
    ecos_ptr::Ptr{ECOS.Cpwork}
    A::AbstractMatrix{Float64}
    b::Vector{Float64}
    n_ecos::Int
    n::Int
    m::Int
end

"""
Solve the constrained least-squares problem

    minimize  ||A1x-b1||₂ st A2x=b2, x≥0.

Rephrases the problem as linear SOC program:

    minimize  [0]'[x]  st  [A1 0  I][x]    [b1] m1
              [1] [t]      [A2 0  0][t] =  [b2] m2
              [0] [r]               [r]
                            n  1  m1
                               x, (t,r) ∈  Rn₊ x (Qn₊)

with dimensions  len(x) = n   len(b1) = m1
                 len(r) = m1  len(b2) = m2

Solve using ECOS, which formulates linear SOC programs as

    minimize  c'x  st  A*x=b, G*x + s = h

Thus, set the ECOS variables as

  c      = [0 , 1,  0]
  A      = [A1  0   I]
            A2  0   0]
  G      = -Identity(n+1+m1)
  h      = [0   0   0]
"""
function ECOSls(A1::AbstractMatrix, b1::Vector{Float64},
                A2::AbstractMatrix, b2::Vector{Float64};
                kwargs...)
    m1, n1 = size(A1)
    m2, n2 = size(A2)
    n1 != n2 && throw(DimensionMismatch("matrix sizes not compatible"))
    n = n1

    n_ecos = n + 1 + m1  # number of vars == len(x)
    m = n_ecos           # number of inequality constraints == len(h)
    p = m1 + m2          # number of equality constraints == len(b)
    l = n                # number of vars in positive orthant
    ncones = 1           # number of SOC cones
    q = [1+m1]           # number of vars in the SOC
    e = 0                # number of exponential cones
    G = ECOS.ECOSMatrix(sparse(-1.0I,n_ecos,n_ecos))
    A = ECOS.ECOSMatrix([A1  spzeros(m1,1)  sparse(1.0I,m1,m1)
                         A2  spzeros(m2,1)      spzeros(m2,m1)])
    c = zeros(n_ecos); c[n+1] = 1.0
    h = zeros(n_ecos)
    b = vcat(b1, b2)

    ecos_ptr = ECOS.setup(n_ecos, m, p, l, ncones, q, e, G, A, c, h, b; kwargs...)
    # problem = ECOSls(ecos_ptr, A1, b1, A2, b2, n_ecos, n, m1, m2)
    ECOS.solve(ecos_ptr)
    # finalizer(cleanup!, problem)
    # return ecos_ptr
end

"""
    ls_simplex(A, b)

Setup the constrained least-squares problem

    minimize  ||Ax-b||₂  subj to  sum(x)=1, x≥0
"""
function ls_simplex(A::AbstractMatrix, b::Vector{Float64}; kwargs...)
    m, n = size(A)
    lsg = ECOSls(A, b, ones(1,n), [1.0]; kwargs...)
    return ECOSlssplx(lsg.ecos_ptr, A, b, lsg.n_ecos, n, m)
end

function solve!(p::ECOSlssplx)

    exitflag = ECOS.solve(p.ecos_ptr)

    if exitflag != ECOS.ECOS_OPTIMAL
        error("did not solve to optimality --- $exitflag")
    end

    # Extract the primal and dual solutions of the ECOS problem
    ecos_prob = unsafe_wrap(Array, p.ecos_ptr, 1)[1]
    xtr = unsafe_wrap(Array, ecos_prob.x, p.n_ecos)[:]
    du_sol_eq = unsafe_wrap(Array, ecos_prob.y, p.m)[:]
    du_sol_ineq = unsafe_wrap(Array, ecos_prob.z, p.n_ecos)[:]
    ECOS.cleanup(p.ecos_ptr)
    # Extract intended solution
    # x = xtr[1:p.n]
    # r = xtr[p.n+2:end]
    # z = du_sol_ineq[1:p.n]
    # μ = du_sol_ineq[p.n+1]

    return xtr
end


text:
lazymap(f, v) = (f(i) for i in v)

# To be replaced by the equivalent operations in TableOperators
_filter(t) = t
function _filter(t, args...)
    mask = [all(i) for i in zip(args...)]
    map(x -> x[mask], Tables.columntable(t))
end

@enum ColumnType categorical numerical arbitrary

const selectordict = Dict(
    categorical => categoricalselector,
    numerical => rangeselector,
    arbitrary => selector,
)

const selectortypes = [:categoricalselector, :rangeselector, :selector]

function hasdistinct(col, n)
    itr = IterTools.distinct(col)
    for (i, _) in enumerate(itr)
        i >= n && return true
    end
    return false
end

struct DefaultType
    threshold::Int
end

(d::DefaultType)(col::AbstractVector{<:Union{Missing, Real}}) = hasdistinct(col, d.threshold) ? numerical : categorical
(d::DefaultType)(col::AbstractVector) = hasdistinct(col, d.threshold) ? arbitrary : categorical

function selectors(t, obs::AbstractObservable; types = Dict(), threshold = 10, default_type = DefaultType(threshold))
    t isa AbstractObservable || (t = Observable{Any}(t))
    cols = map(Tables.columntable, t)
    output = map(x -> Tables.materializer(t[])(x), cols)

    sel_dict = OrderedDict(sym => Observable{Any}(Widget[]) for sym in selectortypes)

    function update_sels!(x)
        for sym in selectortypes
            empty!(sel_dict[sym][])
        end
        for (name, col) in pairs(x)
            sel_func = selectordict[get(types, name, default_type(col))]
            sel = sel_func(col, lazymap)
            push!(sel_dict[widgettype(sel)][], toggled(sel; label = string(name), readout = false))
        end
        for sym in selectortypes
            sel_dict[sym][] = sel_dict[sym][]
        end
    end

    update_sels!(cols[])
    on(update_sels!, cols)

    wdg = Widget{:selectors}(sel_dict; output = output)

    on(obs) do _
        selwdgs = Iterators.flatten(wdg[seltyp][] for seltyp in selectortypes)
        sels = (i[] for i in selwdgs if i[:toggle][])
        output[] = Tables.materializer(t[])(_filter(cols[], sels...))
    end

    layout!(wdg) do x
        sel_cols = [node(
            :div,
            className = "column",
            string(typ),
            @map(node(:div, &x[typ]...))
        ) for typ in selectortypes]
        filters = node(:div, className = "columns", sel_cols...)
    end
end

function selectors(t; kwargs...)
    btn = button("Filter")
    wdg = selectors(t, btn; kwargs...)
    wdg[:filter] = btn
    Widgets.layout(wdg) do x
        node(:div, wdg[:filter], x)
    end
end