text:
abstract type Scattered{T,N} end

# TODO: tests for the data collators and the subsampling
# TODO: Write a version of this that accomodates things that are too big to
# hold in memory
# TODO: write a version of this that uses sparse matrices
# TODO: adapt for complex input and/or complex wavelets
struct ScatteredFull{T,N} <: Scattered{T,N}
    m::Int # number of layers, counting the zeroth layer
    k::Int # the meta-dimension of the signals (should be either 1 or 2)
    data # original last dimension is time, new path/frequency
    output # The final averaged results; this is the output from the entire system
    function ScatteredFull{T,N}(m, k, data, output) where {T<:Number,N}
        @assert N > k # The size of the arrays must be at least 1 larger than the
        # input dimension
        m = Int64(m)
        k = Int64(k)
        # check data
        singlentry = eltype(data[1])
        arrayType = eltype(data)
        @assert arrayType <: AbstractArray{singlentry,N}
        @assert size(data, 1) == m + 1

        # check output
        singlentry = eltype(output[1])
        arrayType = eltype(output)
        @assert arrayType <: AbstractArray{singlentry,N}
        @assert size(output, 1) == m + 1
        new(m, k, data, output)
    end
end
function ScatteredFull(data, output, k=1)
    @assert eltype(data) == eltype(output)
    ScatteredOut{eltype(output),length(output)}(length(output) - 1, k,
        data, output)
end
"""
A simple wrapper for the results of the scattering transform. Its one field
`result` contains a tuple of the results from each layer, in the order zero,
one, two. Say we have an example `s`. You can access layer i by `s[i]`, so
`s[0]` gives the zeroth layer
"""
struct ScatteredOut{T,N} <: Scattered{T,N}
    m::Int# number of layers, counting the zeroth layer
    k::Int# the meta-dimension of the signals (should be either 1 or 2)
    output
end
function ScatteredOut(output, k=1)
    ScatteredOut{eltype(output),length(output)}(length(output) - 1, k,
        output)
end

arrayType(::ScatteredOut{T,N}) where {T,N} = T
noutputs(::ScatteredOut{T,N}) where {T,N} = N

@doc """
      ScatteredOut(m, k, fixDim::Array{<:Real, 1}, n::Array{<:Real, 2}, q::Array{<:Real, 1}, T::DataType)

The resulting output of a scattering transform. `m` is the number of layers, while `k` is the actual signal dimension. `fixDim` is the size of any auxillary dimensions for the input (for example if the input is a list of `nEx` examples of length 100 each, then `fixDim=(nEx,)`, while `n` is a list of the subsampled signal sizes and `q` is a list of the paths, as generated by the scatteringTransform). Finally, the element type is `T`. Initializes all entries to zero.

      ScatteredOut(output, k = 1)
A less involved constructor given just a list or tuple of the output from each layer. `k` gives the signal dimension, as above, with the default that `k=1`.
  """
function ScatteredOut(m, k, fixDim, n, q, T)
    @assert m + 1 == size(n, 1)
    @assert m + 1 == length(q)
    @assert k == size(n, 2)
    fixDim = Int.(fixDim)
    n = Int.(n)
    q = Int.(q)
    N = k + length(fixDim) + 1
    output = [zeros(T, n[i, :]..., prod(q[1:(i-1)] .- 1), fixDim...) for i = 1:m+1]
    return ScatteredOut{T,N}(m, k, output)
end

function ScatteredFull(m, k, fixDim::Array{<:Real,1}, n::Array{<:Real,2},
    q::Array{<:Real,1}, T)
    @assert m + 1 == size(n, 1)
    @assert m + 1 == length(q)
    @assert k == size(n, 2)
    fixDim = Int.(fixDim)
    n = Int.(n)
    q = Int.(q)
    N = k + length(fixDim) + 1
    data = [zeros(T, n[i, :]..., prod(q[1:i] .- 1), fixDim...) for i = 1:m+1]
    output = [zeros(T, n[i, :]..., prod(q[1:(i-1)] .- 1), fixDim...) for i = 1:m+1]
    return ScatteredFull{T,N}(m, k, data, output)
end

import Base.==, Base.+, Base.-, Base.broadcast, Base.length, Base.size, Base.iterate, Base.broadcastable

==(a::ScatteredOut, b::ScatteredOut) = all(a.output .== b.output)
==(a::ScatteredFull, b::ScatteredFull) = all(a.data .== b.data) && all(a.output .== b.output)
function +(a::ScatteredOut, b::ScatteredOut)
    @assert a.m == b.m
    @assert a.k == b.k
    ScatteredOut{arrayType(a),ndims(a)}(a.m, a.k, map((x, y) -> x .+ y, a.output, b.output))
end
function +(a::ScatteredFull, b::ScatteredFull)
    @assert a.m == b.m
    @assert a.k == b.k
    ScatteredFull(a.m, a.k, map((x, y) -> x .+ y, a.data, b.data), map((x, y) -> x .+ y, a.output, b.output))
end
function -(a::ScatteredOut, b::ScatteredOut)
    @assert a.m == b.m
    @assert a.k == b.k
    ScatteredOut{arrayType(a),ndims(a)}(a.m, a.k, map((x, y) -> x - y, a.output, b.output))
end
function -(a::ScatteredFull, b::ScatteredFull)
    @assert a.m == b.m
    @assert a.k == b.k
    ScatteredFull(a.m, a.k, map((x, y) -> x - y, a.data, b.data), map((x, y) -> x - y, a.output, b.output))
end

import Base.iterate, Base.broadcastable, Base.similar
broadcast(f, a::ScatteredOut, varargs...) =
    ScatteredOut(a.m, a.k, map(x -> broadcast(f, x, varargs...), a.output))
broadcast(f, a::ScatteredFull, varargs...) =
    ScatteredFull(a.m, a.k, map(x -> broadcast(f, x, varargs...), a.data),
        map(x -> broadcast(f, x, varargs...), a.output))
broadcastable(a::Scattered) = a.output
iterate(a::Scattered) = iterate(a.output)
iterate(a::Scattered, b) = iterate(a.output, b)
length(sct::Scattered) = length(sct.output) - 1
ndims(sct::Scattered) = sct.k
size(sct::ScatteredOut) = map(x -> size(x), sct.output)
size(sct::ScatteredFull) = map(x -> size(x), sct.data)
similar(sct::ScatteredOut) = ScatteredOut(map(x -> similar(arrayType(sct), axes(x)), sct.output), sct.k)
similar(sct::ScatteredFull) = ScatteredFull(map(x -> similar(x), sct.data), map(x -> similar(x), sct.output), sct.k)

import Statistics.mean
function mean(a::ScatteredOut)
    ScatteredOut(a.m, a.k, map(x -> mean(x, dims=ndims(x)), a.output))
end
import Base.cat

function cat(sc::ScatteredOut, sc1::ScatteredOut, sc2::ScatteredOut...; dims=-1) where {N}
    @assert !any([sc.m != sc1.m, (sc.m .!= map(s -> s.m, sc2))])# check they're
    @assert !any([sc.k != sc1.k, (sc.k .!= map(s -> s.k, sc2))])# all equal
    outputs = map(x -> x.output, (sc, sc1, sc2...))
    output = map(x -> cat(x..., dims=ndims(x[1])), zip(outputs...))
    return ScatteredOut(sc.m, sc.k, tuple(output...))
end
function cat(sc::ScatteredFull, sc1::ScatteredFull, sc2::ScatteredFull...; dims=-1) where {N}
    @assert !any([sc.m != sc1.m, (sc.m .!= map(s -> s.m, sc2))])# check they're
    @assert !any([sc.k != sc1.k, (sc.k .!= map(s -> s.k, sc2))])# all equal
    outputs = map(x -> x.output, (sc, sc1, sc2...))
    output = map(x -> cat(x..., dims=ndims(x[1])), zip(outputs...))
    datas = map(x -> x.data, (sc, sc1, sc2...))
    data = map(x -> cat(x..., dims=ndims(x[1])), zip(datas...))
    return ScatteredFull(sc.m, sc.k, tuple(data...), tuple(output...))
end

function Base.show(io::IO, s::ScatteredOut{T,N}) where {T,N}
    print(io, "ScatteredOut{$T,$N} $(s.k) dim. OutputSizes:")
    outputSizes = map(x -> size(x), s.output)
    for x in outputSizes
        print(io, "\n    $(x)")
    end
end

function Base.show(io::IO, s::ScatteredFull{T,N}) where {T,N}
    print(io, "ScatteredFull{T,N} $(s.k) dim. Data sizes:")
    dataSizes = map(x -> size(x), s.data)
    for x in dataSizes
        print(io, "\n    $(x)")
    end
    print(io, "\nOutputSizes:")
    outputSizes = map(x -> size(x), s.output)
    for x in outputSizes
        print(io, "\n    $(x)")
    end
end

# pathLocs constructors that use a Scattered type
function pathLocs(s::Scattered)
    return pathLocs{length(s.output)}(map(x -> axes(x), s.output))
end

function pathLocs(ii, s::Scattered)
    return pathLocs{length(s.output)}(map(x -> axes(x), s.output[ii]))
end

# access methods
import Base: getindex
Base.getindex(X::Scattered, i::Union{Tuple,<:AbstractArray}) = X.output[i.+1]
function Base.getindex(X::Scattered, ::Colon)
    if ndims(X.output[1]) > ndims(X)
        nEx = size(X.output[1])[end]
        return cat([reshape(lay, (:, nEx)) for lay in X.output]..., dims=1)
    else
        return cat(lay[:] for lay in X.output)
    end

end

Base.getindex(X::Scattered, i::Integer) = X.output[i+1]
function Base.getindex(X::ScatteredOut{T,N}, c::Colon,
    i::Union{<:AbstractArray,<:Integer}) where {T,N}
    ScatteredOut{T,N}(X.m, X.k, map(x -> x[axes(x)[1:end-1]..., i...], X.output))
end
function Base.getindex(X::ScatteredFull{T,N}, c::Colon,
    i::Union{<:AbstractArray,<:Integer}) where {T,N}
    ScatteredFull{T,N}(X.m, X.k, map(x -> x[axes(x)[1:end-1]..., i], X.data),
        map(x -> x[axes(x)[1:end-1]..., i], X.output))
end

function Base.getindex(X::Scattered, p::pathLocs{m}) where {m}
    ijk = p.indices
    λ(ii, access) = X[ii][access...]
    λ(ii, access::Nothing) = nothing
    λ(ii, access::AbstractArray) = X[ii][access]
    res = filter(x -> x != nothing, map(λ, (0:length(ijk)-1), ijk))
    if length(res) == 1
        return res[1]
    else
        return res
    end
end

# setting methods
function Base.setindex!(X::Scattered, v::Number, p::pathLocs{m}) where {m}
    res = X.output
    ijk = p.indices
    for (mm, ind) in enumerate(ijk)
        if typeof(ind) <: BitArray
            res[mm][ind] .= v
        elseif typeof(ind) <: Tuple{Vararg{<:Integer}}
            res[mm][ind...] = v
        elseif ind != nothing
            res[mm][ind...] .= v
        end
    end
end
function Base.setindex!(X::Scattered, v::Tuple, p::pathLocs{m}) where {m}
    res = X.output
    ijk = p.indices
    inputInd = 1
    for (mm, ind) in enumerate(ijk)
        if typeof(ind) <: Tuple{Vararg{<:Integer}}
            res[mm][ind...] = v[inputInd]
            inputInd += 1
        elseif typeof(ind) <: Union{BitArray,Array{Bool}}
            netSize = count(ind)
            res[mm][ind] = v[inputInd:inputInd+netSize-1]
            inputInd += netSize
        elseif ind != nothing
            res[mm][ind...] .= v[inputInd]
            inputInd += 1
        end
    end
end
function Base.setindex!(X::Scattered, v, p::pathLocs{m}) where {m}
    res = X.output
    ijk = p.indices
    inputInd = 1
    for (mm, ind) in enumerate(ijk)
        if typeof(ind) <: Tuple{Vararg{<:Integer}}
            res[mm][ind...] = v[inputInd]
            inputInd += 1
        elseif typeof(ind) <: Union{BitArray,Array{Bool}}
            netSize = count(ind)
            res[mm][ind] = v[inputInd:inputInd+netSize-1]
            inputInd += netSize
        elseif ind != nothing && size(v) != size(res[mm][ind...])
            netSize = prod(size(res[mm][ind...]))
            res[mm][ind...] = reshape(v[inputInd:inputInd+netSize-1], size(res[mm][ind...]))
            inputInd += 1
        elseif ind != nothing
            res[mm][ind...] = v
        end
    end
end

# can't normally union with a colon, which I think is just odd
import Base.union, Base.cat
union(a, b::Colon) = Colon()
union(a::Colon, b) = Colon()
union(a::Colon, b::Colon) = Colon()
# if we want to join several pathLocs
function cat(ps::pathLocs...)
    pairwise(a::Nothing, b::Nothing) = nothing
    pairwise(a::Nothing, b) = b
    pairwise(a, b::Nothing) = a
    pairwise(a, b) = union.(a, b)
    reduce((p1, p2) -> pathLocs(map(pairwise, p1.indices, p2.indices)...), ps)
end

"""
    paths = nonZeroPaths(sc; wholePath=true, allTogetherInOne=false)
given a Scattered, return the pathLocs where the Scattered is nonzero. `wholePath` is true if it returns the whole path, and not just the specific location in the signal. For example, if only `sc(pathLocs(1,(30,2)))` is nonzero, if `wholePath` is true, then `pathLocs(1,(2,))` will be returned while if it is false, `pathLocs(1,(30,2))` will be returned instead.
if `allTogetherInOne` is false, then each location is returned separately, otherwise they are joined into a single `pathLocs`.
"""
function nonZeroPaths(sc; wholePath=true, allTogetherInOne=false)
    if wholePath
        return wholeNonzeroPaths(sc, allTogetherInOne)
    else
        return partNonzeroPaths(sc, allTogetherInOne)
    end
end

function partNonzeroPaths(sc, allTogetherInOne)
    if !allTogetherInOne
        error("wholePath false and each path separate not currently implemented because it makes the backend gross")
    end
    sz = size(sc)
    paths = Tuple{Vararg{pathLocs,0}}()
    Nd = ndims(sc)
    # zeroth layer
    nonZeroLocs = abs.(sc[0]) .> 0
    if any(nonZeroLocs)
        paths = (paths..., pathLocs(0, nonZeroLocs, d=Nd))
    end
    # first layer
    nonZero1 = map(i -> abs.(sc[pathLocs(1, i, d=Nd)]) .> 0, 1:sz[2][end-1])
    nonZero1 = map(x -> reshape(x, (size(x)[1:Nd]..., 1, size(x)[end])), nonZero1)
    nonZero1 = cat(nonZero1..., dims=Nd + 1)
    if any(nonZero1)
        paths = (paths..., pathLocs(1, nonZero1, d=Nd))
    end
    # second layer
    nonZero2 = abs.(sc[pathLocs(2, :)]) .> 0
    if any(nonZero2)
        paths = (paths..., pathLocs(2, nonZero2, d=Nd))
    end
    return cat(paths...)
end

function wholeNonzeroPaths(sc, allTogetherInOne=false)
    sz = size(sc)
    paths = Tuple{Vararg{pathLocs,0}}()
    # zeroth layer
    if maximum(sc[0] .> 0)
        paths = (paths..., pathLocs(0, :))
    end
    # first layer
    nonZero1 = filter(i -> maximum(abs.(sc[pathLocs(1, i)])) > 0, 1:sz[2][end-1])
    paths = (paths..., map(x -> pathLocs(1, x), nonZero1)...)

    # second layer
    nonZero2 = [pathLocs(2, (i, j)) for i = 1:sz[3][end-2], j = 1:sz[3][end-1] if maximum(abs.(sc[pathLocs(2, (i, j))])) > 0]
    paths = (paths..., nonZero2...)
    if allTogetherInOne
        return cat(paths...)
    else
        return paths
    end
end


text:
######################################################################
# pf.jl
# Basic filter described in Chapter 4 of Probo Robo
######################################################################

type PF <: AbstractFilter
	n::Int    # number of particles
	X::Vector{LocTuple}
	W::Vector{Float64}
	Xnew::Vector{LocTuple}
	Wnew::Vector{Float64}
	length::Float64

	function PF(m::SearchDomain, n::Int)
		X = Array(LocTuple, n)
		for i = 1:n
			xi = m.length*rand()
			yi = m.length*rand()
			X[i] = (xi, yi)
		end
		Xnew = deepcopy(X)
		W = zeros(n)
		Wnew = zeros(n)
		return new(n, X, W, Xnew, Wnew, m.length)
	end
end

function update!(pf::PF, x::Vehicle, o::Float64)
	w_sum = 0.0
	p = (x.x, x.y, x.heading)
	for i = 1:pf.n
		# don't need to sample from transition, just calculate weights
		 #pf.W[i] = O(x, pf.X[i], o) # calculate probability of obs
		 pf.W[i] = O(x.sensor, pf.X[i], p, o) # calculate probability of obs
		 w_sum = pf.W[i]
	end

	# Due to deprecation
	# TODO: check that Weights works the same way.
	#  I think it does
	#wvec = WeightVec(pf.W)
	wvec = Weights(pf.W)

	for i = 1:pf.n
		# draw with probability
		#pf.Xnew[i] = sample(pf.X, WeightVec(pf.W))
		j = sample(wvec)
		pf.Xnew[i] = pf.X[j]
		pf.Wnew[i] = pf.W[j]
	end
	copy!(pf.X, pf.Xnew)
	copy!(pf.W, pf.Wnew)
end

function centroid(f::PF)
	wsum = 0.0
	xmean = 0.0
	ymean = 0.0
	for i = 1:f.n
		wval = f.W[i]
		xmean += f.X[i][1] * wval
		ymean += f.X[i][2] * wval
		wsum += wval
	end
	return xmean/wsum, ymean/wsum
end

function reset!(f::PF)
	for i = 1:f.n
		xi = f.length*rand()
		yi = f.length*rand()
		f.X[i] = (xi, yi)
		f.Xnew[i] = (xi, yi)
	end
	fill!(f.W, 0.0)
end

ext: jl


text:
module ShirleyRayTracer

using StaticArrays
using LinearAlgebra
using Images
using CUDA
using KernelAbstractions
using CUDAKernels
using Tullio
import Adapt

const Vec3 = SVector{3, Float64}
const Point3 = SVector{3, Float64}
const Color = RGB{Float64}

export Scene, Camera, Point3, Vec3
export trace_scanline, render_pixel, render, render!
export magnitude, add!, randf

magnitude(x,y) = sqrt(x^2 + y^2)
magnitude(x,y,z) = sqrt(x^2 + y^2 + z^2)
magnitude(v) = magnitude(v...)

magnitudesq(x,y) = x^2 + y^2
magnitudesq(x,y,z) = x^2 + y^2 + z^2
magnitudesq(v) = magnitudesq(v...)

randf(fmin, fmax) = fmin + (fmax-fmin)*rand()
near_zero(v) = v.x < 1e-8 && v.y < 1e-8 && v.z < 1e-8

function random_in_unit_disk()
	x,y = randf(-1, 1), randf(-1, 1)
	while magnitudesq(x,y) >= 1
     	x,y = randf(-1, 1), randf(-1, 1)
	end
	x,y
end

function random_in_unit_sphere()
	x,y,z = randf(-1,1), randf(-1,1), randf(-1,1)
	while magnitudesq(x,y,z) >= 1
		x,y,z = randf(-1,1), randf(-1,1), randf(-1,1)
	end
	Point3(x,y,z)
end

random_unit_vector() = normalize(random_in_unit_sphere())

function random_in_hemisphere(normal)
    in_unit_sphere = random_in_unit_sphere()
    dot(in_unit_sphere, normal) > 0.0 ? in_unit_sphere : -in_unit_sphere
end

function refract(uv, n, etai_over_etat)
    cos_theta = min(dot(-uv, n), 1.0)
    r_out_perp = etai_over_etat * (uv + cos_theta*n)
    r_out_parallel = -sqrt(abs(1.0 - magnitudesq(r_out_perp))) * n
    r_out_perp + r_out_parallel
end

struct Ray
	origin::Point3
	direction::Vec3
	udirection::Vec3
	tm::Float64
	Ray(o, d, m) = new(o, d, normalize(d), m)
	Ray(o, d) = Ray(o, d, 0)
	Ray() = new(Vec3(Inf, Inf, Inf), Vec3(Inf, Inf, Inf)) # a type stable sentinel instead of using nothing
end

at(r::Ray, t) = r.origin + t * r.direction

struct Camera
	origin::Point3
	lower_left_corner::Point3
	horizontal::Vec3
	vertical::Vec3
	u::Vec3
	v::Vec3
	w::Vec3
	lens_radius::Float64
	time0::Float64
	time1::Float64
	function Camera(lookfrom, lookat, vup, vfov, aspect_ratio, aperture, focus_dist, time0=0, time1=0)
		viewport_height = 2.0 * tan(deg2rad(vfov)/2)
		viewport_width = aspect_ratio * viewport_height

		w = normalize(lookfrom - lookat)
		u = normalize(cross(vup, w))
		v = cross(w, u)

		origin = lookfrom
		horizontal = focus_dist * viewport_width * u
		vertical = focus_dist * viewport_height * v
		lower_left_corner = origin - horizontal/2 - vertical/2 - focus_dist*w

		lens_radius = aperture / 2
		new(origin, lower_left_corner, horizontal, vertical, u, v, w, lens_radius, time0, time1)
	end
	Camera() = Camera(Point3(0,0,-1), Point3(0,0,0), Vec3(0,1,0), 40, 1, 0, 10)
end

struct Hit
	p::Point3
	normal::Vec3
	t::Float64
	front_face::Bool
end

abstract type Material end
abstract type Hitable end

include("Materials.jl")
include("Hitables.jl")

struct Scene{H}
	camera::Camera
	hitables::H
end

Adapt.adapt_structure(to, x::Scene) = Scene(Adapt.adapt(to, x.camera), Adapt.adapt(to, x.hitables))

get_ray(scene::Scene, s, t) = get_ray(scene.camera, s, t)

function get_ray(cam::Camera, s, t)
	x, y = cam.lens_radius .* random_in_unit_disk()
	offset = cam.u * x + cam.v * y
	Ray(cam.origin + offset, cam.lower_left_corner + s * cam.horizontal + t * cam.vertical - cam.origin - offset, randf(cam.time0, cam.time1))
end

function trace(scene::Scene, ray::Ray, t_min::Float64, t_max::Float64)
	closest_t = Inf
	struck = 0

	for i in 1:length(scene.hitables)
		t = trace(scene.hitables[i], ray, t_min, closest_t)
		if t >= 0
			closest_t = t
			struck = i
		end
	end
	struck, closest_t
end

function ray_color(scene::Scene, ray::Ray, depth)
	rgb = RGB(1.0,1.0,1.0)

	while depth > 0
		struck, t = trace(scene, ray, 0.001, Inf)
		if struck == 0
			t = 0.5*(ray.udirection.y + 1.0)
			t1m = 1.0 - t
			rgb = rgb ⊙ RGB(t1m + 0.5t, t1m + 0.7t, t1m + t)
			return rgb
		end

		obj = scene.hitables[struck]
		ray, a = scatter(obj.material, ray, hit(obj, t, ray))
		if ray.origin.x == Inf
			break
		end
		rgb = rgb ⊙ a

		depth -= 1
	end
	return RGB(0.0,0.0,0.0)
end

rgb(r, g, b) = RGB(clamp(sqrt(r), 0, 1), clamp(sqrt(g), 0, 1), clamp(sqrt(b), 0, 1))

function render_pixel(x, y, scene, height, width, nsamples, max_depth)
	c = RGB(0.0,0.0,0.0)
	for _ in 1:nsamples
		c += ray_color(scene, get_ray(scene, (x + rand()) / width, (y + rand()) / height), max_depth)
	end
	c /= nsamples
	rgb(c.r, c.g, c.b)
end

function render!(image, scene, width, height, nsamples=10, max_depth=50)
	@tullio image[y, x] = render_pixel(x, height-y+1, scene, height, width, nsamples, max_depth)
	image
end

function render(A::Type{<:AbstractArray}, scene::Scene, width, height, nsamples=10, max_depth=50)
	image = A{RGB{Float64}, 2}(undef, height, width)
	render!(image, scene, width, height, nsamples, max_depth)
	image
end

function render(scene::Scene, width, height, nsamples=10, max_depth=50)
	render(Array, image, scene, width, height, nsamples, max_depth)
end

###
end


text:
# Autogenerated wrapper script for SDL2_ttf_jll for aarch64-linux-musl
export libsdl2_ttf

using Bzip2_jll
using FreeType2_jll
using Glib_jll
using Graphite2_jll
using libpng_jll
using PCRE_jll
using SDL2_jll
JLLWrappers.@generate_wrapper_header("SDL2_ttf")
JLLWrappers.@declare_library_product(libsdl2_ttf, "libSDL2_ttf-2.0.so.0")
function __init__()
    JLLWrappers.@generate_init_header(Bzip2_jll, FreeType2_jll, Glib_jll, Graphite2_jll, libpng_jll, PCRE_jll, SDL2_jll)
    JLLWrappers.@init_library_product(
        libsdl2_ttf,
        "lib/libSDL2_ttf.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
const unary_operators = [:(Base.(:+)),
                         :(Base.(:-)),
                         :(Base.(:!)),
                         :(Base.(:*))]

const numeric_unary_operators = [:(Base.(:+)),
                                 :(Base.(:-))]

const logical_unary_operators = [:(Base.(:!))]

const elementary_functions = [:(Base.abs),
                              :(Base.sign),
                              :(Base.acos),
                              :(Base.acosh),
                              :(Base.asin),
                              :(Base.asinh),
                              :(Base.atan),
                              :(Base.atanh),
                              :(Base.sin),
                              :(Base.sinh),
                              :(Base.conj),
                              :(Base.cos),
                              :(Base.cosh),
                              :(Base.tan),
                              :(Base.tanh),
                              :(Base.ceil),
                              :(Base.floor),
                              :(Base.round),
                              :(Base.trunc),
                              :(Base.exp),
                              :(Base.exp2),
                              :(Base.expm1),
                              :(Base.log),
                              :(Base.log10),
                              :(Base.log1p),
                              :(Base.log2),
                              :(Base.exponent),
                              :(Base.sqrt),
                              :(Base.gamma),
                              :(Base.lgamma),
                              :(Base.digamma),
                              :(Base.erf),
                              :(Base.erfc)]

const two_argument_elementary_functions = [:(Base.round),
                                           :(Base.ceil),
                                           :(Base.floor),
                                           :(Base.trunc)]

const special_comparison_operators = [:(Base.isless)]

const comparison_operators = [:(Base.(:(==))),
                              :(Base.(:(.==))),
                              :(Base.(:(!=))),
                              :(Base.(:(.!=))),
                              :(Base.(:(>))),
                              :(Base.(:(.>))),
                              :(Base.(:(>=))),
                              :(Base.(:(.>=))),
                              :(Base.(:(<))),
                              :(Base.(:(.<))),
                              :(Base.(:(<=))),
                              :(Base.(:(.<=)))]

const scalar_comparison_operators = [:(Base.(:(==))),
                                     :(Base.(:(!=))),
                                     :(Base.(:(>))),
                                     :(Base.(:(>=))),
                                     :(Base.(:(<))),
                                     :(Base.(:(<=)))]

const array_comparison_operators = [:(Base.(:(.==))),
                                    :(Base.(:(.!=))),
                                    :(Base.(:(.>))),
                                    :(Base.(:(.>=))),
                                    :(Base.(:(.<))),
                                    :(Base.(:(.<=)))]

const vectorized_comparison_operators = [:(Base.(:(.==))),
                                         :(Base.(:(==))),
                                         :(Base.(:(.!=))),
                                         :(Base.(:(!=))),
                                         :(Base.(:(.>))),
                                         :(Base.(:(>))),
                                         :(Base.(:(.>=))),
                                         :(Base.(:(>=))),
                                         :(Base.(:(.<))),
                                         :(Base.(:(<))),
                                         :(Base.(:(.<=))),
                                         :(Base.(:(<=)))]

const binary_operators = [:(Base.(:+)),
                          :(Base.(:.+)),
                          :(Base.(:-)),
                          :(Base.(:.-)),
                          :(Base.(:*)),
                          :(Base.(:.*)),
                          :(Base.(:/)),
                          :(Base.(:./)),
                          :(Base.(:.^)),
                          :(Base.div),
                          :(Base.mod),
                          :(Base.fld),
                          :(Base.rem)]

const induced_binary_operators = [:(Base.(:^))]

const arithmetic_operators = [:(Base.(:+)),
                              :(Base.(:.+)),
                              :(Base.(:-)),
                              :(Base.(:.-)),
                              :(Base.(:*)),
                              :(Base.(:.*)),
                              :(Base.(:/)),
                              :(Base.(:./)),
                              :(Base.(:.^)),
                              :(Base.div),
                              :(Base.mod),
                              :(Base.fld),
                              :(Base.rem)]

const induced_arithmetic_operators = [:(Base.(:^))]

const biscalar_operators = [:(Base.maximum),
                            :(Base.minimum)]

const scalar_arithmetic_operators = [:(Base.(:+)),
                                     :(Base.(:-)),
                                     :(Base.(:*)),
                                     :(Base.(:/)),
                                     :(Base.div),
                                     :(Base.mod),
                                     :(Base.fld),
                                     :(Base.rem)]

const induced_scalar_arithmetic_operators = [:(Base.(:^))]

const array_arithmetic_operators = [:(Base.(:+)),
                                    :(Base.(:.+)),
                                    :(Base.(:-)),
                                    :(Base.(:.-)),
                                    :(Base.(:.*)),
                                    :(Base.(:.^))]

const bit_operators = [:(Base.(:&)),
                       :(Base.(:|)),
                       :(Base.(:$))]

const unary_vector_operators = [:(Base.minimum),
                                :(Base.maximum),
                                :(Base.prod),
                                :(Base.sum),
                                :(Base.mean),
                                :(Base.median),
                                :(Base.std),
                                :(Base.var),
                                :(StatsBase.mad),
                                :(Base.norm),
                                :(StatsBase.skewness),
                                :(StatsBase.kurtosis)]

# TODO: dist, iqr

const pairwise_vector_operators = [:(Base.diff)]
                                   #:(Base.reldiff),
                                   #:(Base.percent_change)]

const cumulative_vector_operators = [:(Base.cumprod),
                                     :(Base.cumsum),
                                     :(Base.cumsum_kbn),
                                     :(Base.cummin),
                                     :(Base.cummax)]

const ffts = [:(Base.fft)]

const binary_vector_operators = [:(Base.dot),
                                 :(Base.cor),
                                 :(Base.cov),
                                 :(StatsBase.corspearman)]

const rowwise_operators = [:rowminimums,
                           :rowmaxs,
                           :rowprods,
                           :rowsums,
                           :rowmeans,
                           :rowmedians,
                           :rowstds,
                           :rowvars,
                           :rowffts,
                           :rownorms]

const columnar_operators = [:colminimums,
                            :colmaxs,
                            :colprods,
                            :colsums,
                            :colmeans,
                            :colmedians,
                            :colstds,
                            :colvars,
                            :colffts,
                            :colnorms]

const boolean_operators = [:(Base.any),
                           :(Base.all)]

# Swap arguments to fname() anywhere in AST. Returns the number of
# arguments swapped
function swapargs(ast::Expr, fname::Union(Expr, Symbol))
    if ast.head == :call &&
       (ast.args[1] == fname ||
        (isa(ast.args[1], Expr) && ast.args[1].head == :curly &&
         ast.args[1].args[1] == fname)) &&
       length(ast.args) == 3

        ast.args[2], ast.args[3] = ast.args[3], ast.args[2]
        1
    else
        n = 0
        for arg in ast.args
            n += swapargs(arg, fname)
        end
        n
    end
end
function swapargs(ast, fname::Union(Expr, Symbol))
    ast
    0
end

# Return a block consisting of both the given function and a copy of
# the function in which arguments to the function itself and any
# 2-argument calls to a function of the same name are swapped
macro swappable(func, syms...)
    if (func.head != :function && func.head != :(=)) ||
       func.args[1].head != :call || length(func.args[1].args) != 3
        error("@swappable may only be applied to functions of two arguments")
    end

    func2 = deepcopy(func)
    fname = func2.args[1].args[1]
    if isa(fname, Expr) && fname.head == :curly
        fname = fname.args[1]
    end

    for s in unique([fname, syms...])
        if swapargs(func2, s) < 1
            error("No argument swapped")
        end
    end
    esc(Expr(:block, func, func2))
end

#
# Unary operator macros for DataArrays
#

# Apply unary operator to non-NA members of a DataArray or
# AbstractDataArray
macro dataarray_unary(f, intype, outtype, N...)
    esc(quote
        function $(f){T<:$(intype)}(d::$(isempty(N) ? :(DataArray{T}) : :(DataArray{T,$(N[1])})))
            data = d.data
            res = similar(data, $(outtype))
            @bitenumerate d.na i na begin
                if !na
                    @inbounds res[i] = $(f)(data[i])
                end
            end
            DataArray(res, copy(d.na))
        end
        function $(f){T<:$(intype)}(adv::$(isempty(N) ? :(AbstractDataArray{T}) : :(AbstractDataArray{T,$(N[1])})))
            res = similar(adv, $(outtype))
            for i = 1:length(adv)
                res[i] = ($f)(adv[i])
            end
            res
        end
    end)
end

#
# Binary operator macros for DataArrays
#

# Binary operators with one scalar argument
macro dataarray_binary_scalar(vectorfunc, scalarfunc, outtype, swappable)
    esc(Expr(:block,
        # DataArray and AbstractDataArray with scalar
        # XXX It would be really nice to make this work with arbitrary
        # types, but doing so results in a bunch of method ambiguity
        # warnings
        {
            begin
                fns = {
                    :(function $(vectorfunc)(a::DataArray, b::$t)
                        data = a.data
                        res = similar(data, $outtype)
                        @bitenumerate a.na i na begin
                            if !na
                                @inbounds res[i] = $(scalarfunc)(data[i], b)
                            end
                        end
                        DataArray(res, copy(a.na))
                    end),
                    :(function $(vectorfunc)(a::AbstractDataArray, b::$t)
                        res = similar(a, $outtype)
                        for i = 1:length(a)
                            res[i] = $(scalarfunc)(a[i], b)
                        end
                        res
                    end)
                }
                if swappable
                    # For /, Array/Number is valid but not Number/Array
                    # All other operators should be swappable
                    map!(x->Expr(:macrocall, symbol("@swappable"), x, scalarfunc), fns)
                end
                Expr(:block, fns...)
            end
            for t in (:String, :Number)
        }...
    ))
end

# Binary operators with two array arguments
macro dataarray_binary_array(vectorfunc, scalarfunc, outtype)
    esc(Expr(:block,
        # DataArray with other array
        {
            quote
                function $(vectorfunc)(a::$atype, b::$btype)
                    data1 = $(atype == :DataArray || atype == :(DataArray{Bool}) ? :(a.data) : :a)
                    data2 = $(btype == :DataArray || btype == :(DataArray{Bool}) ? :(b.data) : :b)
                    res = Array($outtype, promote_shape(size(a), size(b)))
                    resna = $narule
                    @bitenumerate resna i na begin
                        if !na
                            @inbounds res[i] = $(scalarfunc)(data1[i], data2[i])
                        end
                    end
                    DataArray(res, resna)
                end
            end
            for (atype, btype, narule) in ((:(DataArray{Bool}), :(DataArray{Bool}), :(a.na | b.na)),
                                           (:(DataArray{Bool}), :(AbstractArray{Bool}), :(copy(a.na))),
                                           (:(AbstractArray{Bool}), :(DataArray{Bool}), :(copy(b.na))),
                                           (:DataArray, :DataArray, :(a.na | b.na)),
                                           (:DataArray, :AbstractArray, :(copy(a.na))),
                                           (:AbstractArray, :DataArray, :(copy(b.na))))
        }...,
        # AbstractDataArray with other array
        # Definitinons with DataArray necessary to avoid ambiguity
        {
            quote
                function $(vectorfunc)(a::$atype, b::$btype)
                    res = similar($(asim ? :a : :b), $outtype, promote_shape(size(a), size(b)))
                    for i = 1:length(a)
                        res[i] = $(scalarfunc)(a[i], b[i])
                    end
                    res
                end
            end
            for (asim, atype, btype) in ((true, :(DataArray{Bool}), :(AbstractDataArray{Bool})),
                                         (false, :(AbstractDataArray{Bool}), :(DataArray{Bool})),
                                         (true, :(AbstractDataArray{Bool}), :(AbstractDataArray{Bool})),
                                         (true, :(AbstractDataArray{Bool}), :(AbstractArray{Bool})),
                                         (false, :(AbstractArray{Bool}), :(AbstractDataArray{Bool})),
                                         (true, :DataArray, :AbstractDataArray),
                                         (false, :AbstractDataArray, :DataArray),
                                         (true, :AbstractDataArray, :AbstractDataArray),
                                         (true, :AbstractDataArray, :AbstractArray),
                                         (false, :AbstractArray, :AbstractDataArray))
        }...,
    ))
end

# Unary operators, NA
for f in unary_operators
    @eval $(f)(d::NAtype) = NA
end

# Unary operators, DataArrays.
@dataarray_unary Base.(:(-)) Bool Int
@dataarray_unary Base.(:(-)) Any T
@dataarray_unary Base.(:(!)) Bool T

# Treat ctranspose and * in a special way
for (f, elf) in ((:(Base.ctranspose), :conj), (:(Base.transpose), :identity))
    @eval begin
        function $(f){T}(d::Union(DataVector{T}, DataMatrix{T}))
            # (c)transpose in Base uses a cache-friendly algorithm for
            # numeric arrays, which is faster than our naive algorithm,
            # but chokes on undefined values in the data array.
            # Fortunately, undefined values can only be present in
            # arrays of non-bits types.
            if isbits(T)
                DataArray($(f)(d.data), d.na.')
            else
                data = d.data
                sz = (size(data, 1), size(data, 2))
                res = similar(data, size(data, 2), size(data, 1))
                @bitenumerate d.na i na begin
                    if !na
                        jnew, inew = ind2sub(sz, i)
                        @inbounds res[inew, jnew] = $(elf)(data[i])
                    end
                end
                DataArray(res, d.na.')
            end
        end
    end
end

# Propagates NA's
# For a dissenting view,
# http://radfordneal.wordpress.com/2011/05/21/slowing-down-matrix-multiplication-in-r/
# But we're getting 10x R while maintaining NA's
for (adata, bdata) in ((true, false), (false, true), (true, true))
    @eval begin
        function Base.(:*)(a::$(adata ? :(Union(DataVector, DataMatrix)) : :(Union(Vector, Matrix))),
                           b::$(bdata ? :(Union(DataVector, DataMatrix)) : :(Union(Vector, Matrix))))
            c = $(adata ? :(a.data) : :a) * $(bdata ? :(b.data) : :b)
            res = DataArray(c, falses(size(c)))
            # Propagation can be made more efficient by storing record of corrupt
            # rows and columns, then doing fast edits.
            $(if adata
                quote
                    n1 = size(a, 1)
                    p1 = size(a, 2)
                    corrupt_rows = falses(n1)
                    for j in 1:p1, i in 1:n1
                        # Propagate NA's
                        # Corrupt all rows based on i
                        corrupt_rows[i] |= a.na[i, j]
                    end
                    res.na[corrupt_rows, :] = true
                end
            end)
            $(if bdata
                quote
                    n2 = size(b, 1)
                    p2 = size(b, 2)
                    corrupt_cols = falses(p2)
                    for j in 1:p2, i in 1:n2
                        # Propagate NA's
                        # Corrupt all columns based on j
                        corrupt_cols[j] |= b.na[i, j]
                    end
                    res.na[:, corrupt_cols] = true
                end
            end)
            res
        end
    end
end

#
# Elementary functions
#
# XXX: The below should be revisited once we have a way to infer what
# the proper return type of an array should be.

# One-argument elementary functions that return the same type as their
# inputs
for f in (:(Base.abs), :(Base.conj), :(Base.sign))
    @eval begin
        $(f)(::NAtype) = NA
        @dataarray_unary $(f) Number T
    end
end

# One-argument elementary functions that always return floating points
for f in (:(Base.acos), :(Base.acosh), :(Base.asin), :(Base.asinh), :(Base.atan), :(Base.atanh),
          :(Base.sin), :(Base.sinh), :(Base.cos), :(Base.cosh), :(Base.tan), :(Base.tanh),
          :(Base.exp), :(Base.exp2), :(Base.expm1), :(Base.log), :(Base.log10), :(Base.log1p),
          :(Base.log2), :(Base.exponent), :(Base.sqrt), :(Base.gamma), :(Base.lgamma),
          :(Base.digamma), :(Base.erf), :(Base.erfc))
    @eval begin
        ($f)(::NAtype) = NA
        @dataarray_unary $(f) FloatingPoint T
        @dataarray_unary $(f) Real Float64
    end
end

# Elementary functions that take varargs
for f in (:(Base.round), :(Base.ceil), :(Base.floor), :(Base.trunc))
    @eval begin
        ($f)(::NAtype, args::Integer...) = NA

        # ambiguity
        @dataarray_unary $(f) Real T 1
        @dataarray_unary $(f) Real T 2
        @dataarray_unary $(f) Real T

        function $(f){T<:Real}(d::DataArray{T}, args::Integer...)
            data = similar(d.data)
            @bitenumerate d.na i na begin
                if !na
                    @inbounds data[i] = $(f)(d[i], args...)
                end
            end
            DataArray(data, copy(d.na))
        end
        function $(f){T<:Real}(adv::AbstractDataArray{T}, args::Integer...)
            res = similar(adv)
            for i = 1:length(adv)
                res[i] = ($f)(adv[i], args...)
            end
            res
        end
    end
end

#
# Bit operators
#

@swappable Base.(:&)(a::NAtype, b::Bool) = b ? NA : false
@swappable Base.(:|)(a::NAtype, b::Bool) = b ? true : NA
@swappable Base.(:$)(a::NAtype, b::Bool) = NA

# To avoid ambiguity warning
@swappable Base.(:|)(a::NAtype, b::Function) = NA

for f in (:(Base.(:&)), :(Base.(:|)), :(Base.(:$)))
    @eval begin
        # Scalar with NA
        ($f)(::NAtype, ::NAtype) = NA
        @swappable ($f)(::NAtype, b::Integer) = NA
    end
end

# DataArray with DataArray
Base.(:&)(a::DataArray{Bool}, b::DataArray{Bool}) =
    DataArray(a.data & b.data, (a.na & b.data) | (b.na & a.data))
Base.(:|)(a::DataArray{Bool}, b::DataArray{Bool}) =
    DataArray(a.data | b.data, (a.na & !b.data) | (b.na & !a.data))
Base.(:$)(a::DataArray{Bool}, b::DataArray{Bool}) =
    DataArray(a.data $ b.data, a.na | b.na)

# DataArray with non-DataArray
# Need explicit definition for BitArray to avoid ambiguity
for t in (:(BitArray), :(Union(AbstractArray{Bool}, Bool)))
    @eval begin
        @swappable Base.(:&)(a::DataArray{Bool}, b::$t) = DataArray(convert(Array{Bool}, a.data & b), a.na & b)
        @swappable Base.(:|)(a::DataArray{Bool}, b::$t) = DataArray(convert(Array{Bool}, a.data | b), a.na & !b)
        @swappable Base.(:$)(a::DataArray{Bool}, b::$t) = DataArray(convert(Array{Bool}, a.data $ b), copy(a.na))
    end
end

#
# Comparison operators
#

Base.isless(::NAtype, ::NAtype) = false
Base.isless(::NAtype, b) = true
Base.isless(a, ::NAtype) = false

# This is for performance only; the definition in Base is sufficient
# for AbstractDataArrays
function Base.isequal(a::DataArray, b::DataArray)
    if size(a) != size(b) || a.na != b.na
        return false
    end
    @bitenumerate a.na i na begin
        @inbounds if !na && !isequal(a.data[i], b.data[i])
            return false
        end
    end
    return true
end

# ambiguity
@swappable ==(a::DataArray{Bool}, b::BitArray) =
    invoke(==, (DataArray, AbstractArray), a, b)
@swappable ==(a::DataArray, b::BitArray) =
    invoke(==, (DataArray, AbstractArray), a, b)
@swappable ==(a::AbstractDataArray{Bool}, b::BitArray) =
    invoke(==, (DataArray, AbstractArray), a, b)
@swappable ==(a::AbstractDataArray, b::BitArray) =
    invoke(==, (DataArray, AbstractArray), a, b)

function Base.(:(==))(a::DataArray, b::DataArray)
    size(a) == size(b) || return false
    adata = a.data
    bdata = b.data
    bchunks = b.na.chunks
    has_na = false
    @bitenumerate a.na i na begin
        if na || Base.unsafe_bitgetindex(bchunks, i)
            has_na = true
        else
            @inbounds adata[i] == bdata[i] || return false
        end
    end
    has_na ? NA : true
end

# ambiguity
@swappable ==(a::DataArray, b::AbstractDataArray) =
    invoke(==, (AbstractDataArray, AbstractDataArray), a, b)

@swappable function Base.(:(==))(a::DataArray, b::AbstractArray)
    size(a) == size(b) || return false
    adata = a.data
    has_na = false
    @bitenumerate a.na i na begin
        if na
            has_na = true
        else
            @inbounds adata[i] == b[i] || return false
        end
    end
    has_na ? NA : true
end

function Base.(:(==))(a::AbstractDataArray, b::AbstractDataArray)
    size(a) == size(b) || return false
    has_na = false
    for i = 1:length(a)
        if isna(a[i]) || isna(b[i])
            has_na = true
        else
            a[i] == b[i] || return false
        end
    end
    has_na ? NA : true
end

@swappable function Base.(:(==))(a::AbstractDataArray, b::AbstractArray)
    size(a) == size(b) || return false
    has_na = false
    for i = 1:length(a)
        if isna(a[i])
            has_na = true
        else
            a[i] == b[i] || return false
        end
    end
    has_na ? NA : true
end

for (sf,vf) in zip(scalar_comparison_operators, array_comparison_operators)
    @eval begin
        # Array with NA
        @swappable ($(vf)){T,N}(::NAtype, b::AbstractArray{T,N}) =
            DataArray(Array(Bool, size(b)), trues(size(b)))

        # Scalar with NA
        ($(vf))(::NAtype, ::NAtype) = NA
        ($(sf))(::NAtype, ::NAtype) = NA
        @swappable ($(vf))(::NAtype, b) = NA
        @swappable ($(sf))(::NAtype, b) = NA

        @dataarray_binary_scalar $(vf) $(sf) Bool true
        @dataarray_binary_array $(vf) $(sf) Bool
    end
end

#
# Binary operators
#

# Necessary to avoid ambiguity warnings
Base.(:.^)(::MathConst{:e}, B::DataArray) = exp(B)
Base.(:.^)(::MathConst{:e}, B::AbstractDataArray) = exp(B)

for f in (:(Base.(:+)), :(Base.(:.+)), :(Base.(:-)), :(Base.(:.-)),
          :(Base.(:*)), :(Base.(:.*)), :(Base.(:.^)), :(Base.div),
          :(Base.mod), :(Base.fld), :(Base.rem))
    @eval begin
        # Scalar with NA
        ($f)(::NAtype, ::NAtype) = NA
        @swappable ($f)(d::NAtype, x::Number) = NA
    end
end

# Define methods for UniformScaling. Otherwise we get ambiguity
# warnings...
function +{TA,TJ}(A::DataArray{TA,2},J::UniformScaling{TJ})
    n = Base.LinAlg.chksquare(A)
    B = similar(A,promote_type(TA,TJ))
    copy!(B,A)
    @inbounds for i = 1:n
        if !B.na[i,i]
            B.data[i,i] += J.λ
        end
    end
    B
end
+{TA}(J::UniformScaling,A::DataArray{TA,2}) = A + J

function -{TA,TJ<:Number}(A::DataArray{TA,2},J::UniformScaling{TJ})
    n = Base.LinAlg.chksquare(A)
    B = similar(A,promote_type(TA,TJ))
    copy!(B,A)
    @inbounds for i = 1:n
        if !B.na[i,i]
            B.data[i,i] -= J.λ
        end
    end
    B
end
function -{TA,TJ<:Number}(J::UniformScaling{TJ},A::DataArray{TA,2})
    n = Base.LinAlg.chksquare(A)
    B = -A
    @inbounds for i = 1:n
        if !B.na[i,i]
            B.data[i,i] += J.λ
        end
    end
    B
end

+(A::DataArray{Bool,2},J::UniformScaling{Bool}) =
    invoke(+, (AbstractArray{Bool,2}, UniformScaling{Bool}), A, J)
+(J::UniformScaling{Bool},A::DataArray{Bool,2}) =
    invoke(+, (UniformScaling{Bool}, AbstractArray{Bool,2}), J, A)
-(A::DataArray{Bool,2},J::UniformScaling{Bool}) =
    invoke(-, (AbstractArray{Bool,2}, UniformScaling{Bool}), A, J)
-(J::UniformScaling{Bool},A::DataArray{Bool,2}) =
    invoke(-, (UniformScaling{Bool}, AbstractArray{Bool,2}), J, A)

+{TA,TJ}(A::AbstractDataArray{TA,2},J::UniformScaling{TJ}) =
    invoke(+, (AbstractArray{TA,2}, UniformScaling{TJ}), A, J)
+{TA}(J::UniformScaling,A::AbstractDataArray{TA,2}) =
    invoke(+, (UniformScaling, AbstractArray{TA,2}), J, A)
-{TA,TJ<:Number}(A::AbstractDataArray{TA,2},J::UniformScaling{TJ}) =
    invoke(-, (AbstractArray{TA,2}, UniformScaling{TJ}), A, J)
-{TA,TJ<:Number}(J::UniformScaling{TJ},A::AbstractDataArray{TA,2}) =
    invoke(-, (UniformScaling{TJ}, AbstractArray{TA,2}), J, A)

+(A::AbstractDataArray{Bool,2},J::UniformScaling{Bool}) =
    invoke(+, (AbstractArray{Bool,2}, UniformScaling{Bool}), A, J)
+(J::UniformScaling{Bool},A::AbstractDataArray{Bool,2}) =
    invoke(+, (UniformScaling{Bool}, AbstractArray{Bool,2}), J, A)
-(A::AbstractDataArray{Bool,2},J::UniformScaling{Bool}) =
    invoke(-, (AbstractArray{Bool,2}, UniformScaling{Bool}), A, J)
-(J::UniformScaling{Bool},A::AbstractDataArray{Bool,2}) =
    invoke(-, (UniformScaling{Bool}, AbstractArray{Bool,2}), J, A)

for f in (:(Base.(:.+)), :(Base.(:.-)), :(Base.(:*)), :(Base.(:.*)),
          :(Base.(:.^)), :(Base.div), :(Base.mod), :(Base.fld), :(Base.rem))
    @eval begin
        # Array with NA
        @swappable $(f){T,N}(::NAtype, b::AbstractArray{T,N}) =
            DataArray(Array(T, size(b)), trues(size(b)))

        # DataArray with scalar
        @dataarray_binary_scalar $f $f promote_type(eltype(a), eltype(b)) true
    end
end

Base.(:^)(::NAtype, ::NAtype) = NA
Base.(:^)(a, ::NAtype) = NA
Base.(:^)(::NAtype, ::Integer) = NA
Base.(:^)(::NAtype, ::Number) = NA

for (vf, sf) in ((:(Base.(:+)), :(Base.(:+))),
                 (:(Base.(:.+)), :(Base.(:+))),
                 (:(Base.(:-)), :(Base.(:-))),
                 (:(Base.(:.-)), :(Base.(:-))),
                 (:(Base.(:.*)), :(Base.(:*))),
                 (:(Base.(:.^)), :(Base.(:^))))
    @eval begin
        # Necessary to avoid ambiguity warnings
        @swappable ($vf)(A::BitArray, B::AbstractDataArray{Bool}) = ($vf)(bitunpack(A), B)
        @swappable ($vf)(A::BitArray, B::DataArray{Bool}) = ($vf)(bitunpack(A), B)

        @dataarray_binary_array $vf $sf promote_type(eltype(a), eltype(b))
    end
end

@swappable Base.(:./)(A::BitArray, B::AbstractDataArray{Bool}) = ./(bitunpack(A), B)
@swappable Base.(:./)(A::BitArray, B::DataArray{Bool}) = ./(bitunpack(A), B)

# / and ./ are defined separately since they promote to floating point
for f in ((:(Base.(:/)), :(Base.(:./))))
    @eval begin
        ($f)(::NAtype, ::NAtype) = NA
        @swappable ($f)(d::NAtype, x::Number) = NA
    end
end

Base.(:/){T,N}(b::AbstractArray{T,N}, ::NAtype) =
    DataArray(Array(T, size(b)), trues(size(b)))
@dataarray_binary_scalar Base.(:/) Base.(:/) eltype(a) <: FloatingPoint || typeof(b) <: FloatingPoint ?
                                      promote_type(eltype(a), typeof(b)) : Float64 false
@swappable Base.(:./){T,N}(::NAtype, b::AbstractArray{T,N}) =
    DataArray(Array(T, size(b)), trues(size(b)))
@dataarray_binary_scalar Base.(:./) Base.(:/) eltype(a) <: FloatingPoint || typeof(b) <: FloatingPoint ?
                                      promote_type(eltype(a), typeof(b)) : Float64 true
@dataarray_binary_array Base.(:./) Base.(:/) eltype(a) <: FloatingPoint || eltype(b) <: FloatingPoint ?
                                             promote_type(eltype(a), eltype(b)) : Float64

for f in biscalar_operators
    @eval begin
        ($f)(::NAtype, ::NAtype) = NA
        @swappable $(f)(::Number, ::NAtype) = NA
    end
end

for f in pairwise_vector_operators
    @eval function ($f)(dv::DataVector)
        n = length(dv)
        new_data = ($f)(dv.data)
        new_na = falses(n - 1)
        new_na[1] = dv.na[1]
        for i = 2:(n - 1)
            if dv.na[i]
                new_na[i - 1] = true
                new_na[i] = true
            end
        end
        new_na[n - 1] = new_na[n - 1] || dv.na[n]
        return DataArray(new_data, new_na)
    end
end

for f in cumulative_vector_operators
    @eval function ($f)(dv::DataVector)
        new_data = ($f)(dv.data)
        new_na = falses(length(dv))
        hitna = false
        @bitenumerate dv.na i na begin
            hitna |= na
            if hitna
                new_na[i] = true
            end
        end
        return DataArray(new_data, new_na)
    end
end

for f in [unary_vector_operators; ffts]
    @eval ($f)(dv::DataVector) = any(dv.na) ? NA : ($f)(dv.data)
end

for f in binary_vector_operators
    @eval ($f)(dv1::DataVector, dv2::DataVector) =
            any(dv1.na) || any(dv2.na) ? NA : ($f)(dv1.data, dv2.data)
end

for f in (:(Base.minimum), :(Base.maximum), :(Base.prod), :(Base.sum),
          :(Base.mean), :(Base.median), :(Base.std), :(Base.var),
          :(Base.norm))
    colf = symbol("col$(f)s")
    rowf = symbol("row$(f)s")
    @eval begin
        function ($colf)(dm::AbstractDataMatrix)
            n, p = nrow(dm), ncol(dm)
            res = datazeros(p)
            for j in 1:p
                res[j] = ($f)(dm[:, j])
            end
            return res
        end
        function ($rowf)(dm::DataMatrix)
            n, p = nrow(dm), ncol(dm)
            res = datazeros(n)
            for i in 1:n
                res[i] = ($f)(DataArray(reshape(dm.data[i, :], p), reshape(dm.na[i, :], p)))
            end
            return res
        end
    end
end

#
# Boolean operators
#

function Base.all(dv::DataArray{Bool})
    data = dv.data
    has_na = false
    @bitenumerate dv.na i na begin
        if !na
            data[i] || return false
        else
            has_na = true
        end
    end
    has_na ? NA : true
end

function Base.all(dv::AbstractDataArray{Bool})
    has_na = false
    for i in 1:length(dv)
        x = dv[i]
        if !isna(x)
            x || return false
        else
            has_na = true
        end
    end
    has_na ? NA : true
end

function Base.any(dv::DataArray{Bool})
    has_na = false
    @bitenumerate dv.na i na begin
        if !na
            if dv.data[i]
                return true
            end
        else
            has_na = true
        end
    end
    has_na ? NA : false
end

function Base.any(dv::AbstractDataArray{Bool})
    has_na = false
    for i in 1:length(dv)
        if !isna(dv[i])
            if dv[i]
                return true
            end
        else
            has_na = true
        end
    end
    has_na ? NA : false
end

StatsBase.range{T}(dv::AbstractDataVector{T}) = maximum(dv) - minimum(dv)

function rle{T}(v::AbstractVector{T})
    n = length(v)
    current_value = v[1]
    current_length = 1
    values = similar(v, n)
    total_values = 1
    lengths = Array(Int16, n)
    total_lengths = 1
    for i in 2:n
        if v[i] == current_value
            current_length += 1
        else
            values[total_values] = current_value
            total_values += 1
            lengths[total_lengths] = current_length
            total_lengths += 1
            current_value = v[i]
            current_length = 1
        end
    end
    values[total_values] = current_value
    lengths[total_lengths] = current_length
    return (values[1:total_values], lengths[1:total_lengths])
end

function rle{T}(v::AbstractDataVector{T})
    n = length(v)
    current_value = v[1]
    current_length = 1
    values = DataArray(T, n)
    total_values = 1
    lengths = Array(Int16, n)
    total_lengths = 1
    for i in 2:n
        if isna(v[i]) || isna(current_value)
            if isna(v[i]) && isna(current_value)
                current_length += 1
            else
                values[total_values] = current_value
                total_values += 1
                lengths[total_lengths] = current_length
                total_lengths += 1
                current_value = v[i]
                current_length = 1
            end
        else
            if v[i] == current_value
                current_length += 1
            else
                values[total_values] = current_value
                total_values += 1
                lengths[total_lengths] = current_length
                total_lengths += 1
                current_value = v[i]
                current_length = 1
            end
        end
    end
    values[total_values] = current_value
    lengths[total_lengths] = current_length
    return (values[1:total_values], lengths[1:total_lengths])
end

## inverse run-length encoding
function inverse_rle{T, I <: Integer}(values::AbstractVector{T},
                                      lengths::Vector{I})
    total_n = sum(lengths)
    pos = 0
    res = similar(values, total_n)
    n = length(values)
    for i in 1:n
        v = values[i]
        l = lengths[i]
        for j in 1:l
            pos += 1
            res[pos] = v
        end
    end
    return res
end


text:
# basics - data types
# working in command line
# julia> is the prompt
# version: 1.0

# numero uno
julia> 1
1

# numero dos
julia> 2
2

# float
julia> 2.3
2.3

# typeof
julia> typeof(1)
Int64

julia> typeof(2.3)
Float64

# Julia Sys
julia> Sys
Base.Sys

# Julia Sys Length
julia> Sys.length(1)
1

# Julia String declaration
julia> a = "praveen"
"praveen"

julia> a
"praveen"

julia> Sys.length(a)
7


text:
module Example2018

export Point
struct Point{T}
    x::T
    y::T
end

function Base. +(p::Point, q::Point)
    Point(p.x + q.x, p.y + q.y)
end

export MPoint
mutable struct MPoint{T}
    x::T
    y::T
end

###############################################################################

export Polynomial
struct Polynomial{T}
    coeffs::Vector{T}
end

export evaluate
function evaluate(p::Polynomial, x)
    sum(c * x^(i-1) for (i,c) in enumerate(p.coeffs))
end

function Base. ==(p::Polynomial, q::Polynomial)
    d = p - q
    # all(c -> c==0, d.coeffs)
    all(d.coeffs .==0)
end

function Base. +(p::Polynomial, q::Polynomial)
    l = max(length(p.coeffs), length(q.coeffs))
    rcoeffs = zeros(eltype(p.coeffs),l)
    for (i,c) in enumerate(p.coeffs)
        rcoeffs[i] += c
    end
    for (i,c) in enumerate(q.coeffs)
        rcoeffs[i] += c
    end
    Polynomial(rcoeffs)
end

function Base. *(a, p::Polynomial)
    Polynomial(map(x->a*x, p.coeffs))
end
Base. *(p::Polynomial, a) = a*p

function Base. -(p::Polynomial)
    (-1) * p
end

function Base. -(p::Polynomial, q::Polynomial)
    p + (-q)
end

export deriv
function deriv(p::Polynomial)
    Polynomial([(i-1) * c for (i,c) in enumerate(p.coeffs)][2:end])
end

end


text:
# ------------------------------------------------------------------
# Copyright (c) 2018, Júlio Hoffimann Mendes <juliohm@stanford.edu>
# Licensed under the ISC License. See LICENCE in the project root.
# ------------------------------------------------------------------

"""
    PointSetData(data, coords)

Spatial `data` georeferenced with coordinate matrix `coords`.
The `data` argument is a dictionary mapping variable names to Julia
arrays with the actual data.

See also: [`PointSet`](@ref)
"""
struct PointSetData{T<:Real,N} <: AbstractSpatialData{T,N}
  data::Dict{Symbol,<:AbstractArray}
  coords # AbstractMatrix

  function PointSetData{T,N}(data, coords) where {N,T<:Real}
    lens = [length(array) for array in values(data)]
    @assert all(lens .== size(coords, 2)) "data and coords must have the same number of points"
    new(data, coords)
  end
end

PointSetData(data::Dict{Symbol,<:AbstractArray}, coords::AbstractMatrix{T}) where {T<:Real} =
  PointSetData{T,size(coords,1)}(data, coords)

coordinates(geodata::PointSetData{T,N}) where {N,T<:Real} = Dict(Symbol("x$i") => T for i=1:N)

variables(geodata::PointSetData) = Dict(var => eltype(array) for (var,array) in geodata.data)

npoints(geodata::PointSetData) = size(geodata.coords, 2)

function coordinates!(buff::AbstractVector{T}, geodata::PointSetData{T,N}, ind::Int) where {N,T<:Real}
  for i in 1:N
    @inbounds buff[i] = geodata.coords[i,ind]
  end
end

value(geodata::PointSetData, ind::Int, var::Symbol) = geodata.data[var][ind]

function Base.view(geodata::PointSetData{T,N}, inds::AbstractVector{Int}) where {N,T<:Real}
  data = Dict(var => view(array, inds) for (var,array) in geodata.data)
  coords = view(geodata.coords, :, inds)

  PointSetData(data, coords)
end

# ------------
# IO methods
# ------------
function Base.show(io::IO, geodata::PointSetData{T,N}) where {N,T<:Real}
  npts = size(geodata.coords, 2)
  print(io, "$npts PointSetData{$T,$N}")
end

function Base.show(io::IO, ::MIME"text/plain", geodata::PointSetData{T,N}) where {N,T<:Real}
  println(io, geodata)
  println(io, "  variables")
  varlines = ["    └─$var ($(eltype(array)))" for (var, array) in geodata.data]
  print(io, join(varlines, "\n"))
end


text: include("io.jl")


text:
@doc raw"""
    SubGradientProblem <: Problem

A structure to store information about a subgradient based optimization problem

# Fields
* `manifold` – a [Manifold](https://juliamanifolds.github.io/Manifolds.jl/stable/interface.html#ManifoldsBase.Manifold)
* `cost` – the function $F$ to be minimized
* `subgradient` – a function returning a subgradient $\partial F$ of $F$

# Constructor

    SubGradientProblem(M, f, ∂f)

Generate the [`Problem`] for a subgradient problem, i.e. a function `f` on the
manifold `M` and a function `∂f` that returns an element from the subdifferential
at a point.
"""
mutable struct SubGradientProblem{mT <: Manifold} <: Problem
    M::mT
    cost::Function
    subgradient::Function
end
"""
    get_subgradient(p,x)

Evaluate the (sub)gradient of a [`SubGradientProblem`](@ref)` p` at the point `x`.
"""
function get_subgradient(p::P, x) where {P <: SubGradientProblem{M} where M <: Manifold}
    return p.subgradient(x)
end
"""
    SubGradientMethodOptions <: Options
stories option values for a [`subgradient_method`](@ref) solver

# Fields
* `retraction` – the retration to use within
* `stepsize` – a [`Stepsize`](@ref)
* `stop` – a [`StoppingCriterion`](@ref)
* `x` – (initial or current) value the algorithm is at
* `xOptimal` – optimal value
* `∂` the current element from the possivle subgradients at `x` that is used
"""
mutable struct SubGradientMethodOptions{P,T} <: Options where {P,T}
    retract!::Function
    stepsize::Stepsize
    stop::StoppingCriterion
    x::P
    xOptimal::P
    ∂::T
    function SubGradientMethodOptions{P}(
        M::TM,
        x::P,
        sC::StoppingCriterion,
        s::Stepsize,
        retr!::Function=exp!
        ) where {TM <: Manifold, P}
        o = new{P,typeof(zero_tangent_vector(M,x))}();
        o.x = x;
        o.xOptimal = x;
        o.stepsize = s; o.retract! = retr!;
        o.stop = sC;
        return o
    end
end
function SubGradientMethodOptions(
    M::Manifold,
    x::P,
    sC::StoppingCriterion,
    s::Stepsize,
    retr!::Function=exp!
) where {P}
    return SubGradientMethodOptions{P}(M, x, sC, s, retr!)
end