text:
type ConvexShape <: Drawable
    ptr::Ptr{Void}
    _texture::Texture

    function ConvexShape(ptr::Ptr{Void})
        c = new(ptr)
        finalizer(c, destroy)
        c
    end
end

function ConvexShape()
    ConvexShape(ccall((:sfConvexShape_create, libcsfml_graphics), Ptr{Void}, ()))
end

function copy(shape::ConvexShape)
    return ConvexShape(ccall((:sfConvexShape_copy, libcsfml_graphics), Ptr{Void}, (Ptr{Void},), shape.ptr))
end

function destroy(shape::ConvexShape)
    ccall((:sfConvexShape_destroy, libcsfml_graphics), Void, (Ptr{Void},), shape.ptr)
end

function set_position(shape::ConvexShape, pos::Vector2f)
    ccall((:sfConvexShape_setPosition, libcsfml_graphics), Void, (Ptr{Void}, Vector2f,), shape.ptr, pos)
end

function set_rotation(shape::ConvexShape, angle::Real)
    ccall((:sfConvexShape_setRotation, libcsfml_graphics), Void, (Ptr{Void}, Cfloat,), shape.ptr, angle)
end

function set_scale(shape::ConvexShape, scale::Vector2f)
    ccall((:sfConvexShape_setScale, libcsfml_graphics), Void, (Ptr{Void}, Vector2f,), shape.ptr, scale)
end

function set_origin(shape::ConvexShape, origin::Vector2f)
    ccall((:sfConvexShape_setOrigin, libcsfml_graphics), Void, (Ptr{Void}, Vector2f,), shape.ptr, origin)
end

function get_position(shape::ConvexShape)
    return ccall((:sfConvexShape_getPosition, libcsfml_graphics), Vector2f, (Ptr{Void},), shape.ptr)
end

function get_rotation(shape::ConvexShape)
    return Real(ccall((:sfConvexShape_getRotation, libcsfml_graphics), Cfloat, (Ptr{Void},), shape.ptr))
end

function get_scale(shape::ConvexShape)
    return ccall((:sfConvexShape_getScale, libcsfml_graphics), Vector2f, (Ptr{Void},), shape.ptr)
end

function get_origin(shape::ConvexShape)
    return ccall((:sfConvexShape_getOrigin, libcsfml_graphics), Vector2f, (Ptr{Void},), shape.ptr)
end

function move(shape::ConvexShape, offset::Vector2f)
    ccall((:sfConvexShape_move, libcsfml_graphics), Void, (Ptr{Void}, Vector2f), shape.ptr, offset)
end

function rotate(shape::ConvexShape, angle::Real)
    ccall((:sfConvexShape_rotate, libcsfml_graphics), Void, (Ptr{Void}, Cfloat,), shape.ptr, angle)
end

function scale(shape::ConvexShape, factors::Vector2f)
    ccall((:sfConvexShape_scale, libcsfml_graphics), Void, (Ptr{Void}, Vector2f,), shape.ptr, factors)
end

function set_texture(shape::ConvexShape, texture::Texture)
    ccall((:sfConvexShape_setTexture, libcsfml_graphics), Void, (Ptr{Void}, Ptr{Void},), shape.ptr, texture.ptr)
    shape._texture = texture
end

function set_texture_rect(shape::ConvexShape, rect::IntRect)
    ccall((:sfConvexShape_setTextureRect, libcsfml_graphics), Void, (Ptr{Void}, IntRect,), shape.ptr, rect)
end

function set_fillcolor(shape::ConvexShape, color::Color)
    ccall((:sfConvexShape_setFillColor, libcsfml_graphics), Void, (Ptr{Void}, Color,), shape.ptr, color)
end

function set_outlinecolor(shape::ConvexShape, color::Color)
    ccall((:sfConvexShape_setOutlineColor, libcsfml_graphics), Void, (Ptr{Void}, Color,), shape.ptr, color)
end

function set_outline_thickness(shape::ConvexShape, thickness::Real)
    ccall((:sfConvexShape_setOutlineThickness, libcsfml_graphics), Void, (Ptr{Void}, Cfloat,), shape.ptr, thickness)
end

function get_texture(shape::ConvexShape)
    return Texture(ccall((:sfConvexShape_getTexture, libcsfml_graphics), Ptr{Void}, (Ptr{Void},), shape.ptr))
end

function get_texture_rect(shape::ConvexShape)
    return ccall((:sfConvexShape_getTextureRect, libcsfml_graphics), IntRect, (Ptr{Void},), shape.ptr)
end

function get_fillcolor(shape::ConvexShape)
    return ccall((:sfConvexShape_getFillColor, libcsfml_graphics), Color, (Ptr{Void},), shape.ptr)
end

function get_outlinecolor(shape::ConvexShape)
    return ccall((:sfConvexShape_getOutlineColor, libcsfml_graphics), Color, (Ptr{Void},), shape.ptr)
end

function get_outline_thickness(shape::ConvexShape)
    return Real(ccall((:sfConvexShape_getOutlineThickness, libcsfml_graphics), Cfloat, (Ptr{Void},), shape.ptr))
end

function get_pointcount(shape::ConvexShape)
    return Int(ccall((:sfConvexShape_getPointCount, libcsfml_graphics), UInt32, (Ptr{Void},), shape.ptr))
end

function get_point(shape::ConvexShape, index::Integer)
    return ccall((:sfConvexShape_getPoint, libcsfml_graphics), Vector2f, (Ptr{Void}, UInt32,), shape.ptr, index)
end

function set_pointcount(shape::ConvexShape, count::Integer)
    ccall((:sfConvexShape_setPointCount, libcsfml_graphics), Void, (Ptr{Void}, UInt32,), shape.ptr, count)
end

function set_point(shape::ConvexShape, index::Integer, point::Vector2f)
    ccall((:sfConvexShape_setPoint, libcsfml_graphics), Void, (Ptr{Void}, UInt32, Vector2f,), shape.ptr, index, point)
end

function get_localbounds(shape::ConvexShape)
    return ccall((:sfConvexShape_getLocalBounds, libcsfml_graphics), FloatRect, (Ptr{Void},), shape.ptr)
end

function get_globalbounds(shape::ConvexShape)
    return ccall((:sfConvexShape_getGlobalBounds, libcsfml_graphics), FloatRect, (Ptr{Void},), shape.ptr)
end


text:
import LaTeXStrings
import PGFPlotsX

"""
"""
function plotroccurve(
        estimator::AbstractFittable,
        features_df::DataFrames.AbstractDataFrame,
        labels_df::DataFrames.AbstractDataFrame,
        single_label_name::Symbol,
        positive_class::AbstractString;
        kwargs...,
        )
    vectorofestimators = [estimator]
    result = plotroccurve(
        vectorofestimators,
        features_df,
        labels_df,
        single_label_name,
        positive_class;
        kwargs...,
        )
    return result
end

"""
"""
function plotroccurve(
        vectorofestimators::AbstractVector{AbstractFittable},
        features_df::DataFrames.AbstractDataFrame,
        labels_df::DataFrames.AbstractDataFrame,
        single_label_name::Symbol,
        positive_class::AbstractString;
        legend_pos::AbstractString = "outer north east",
        )
    legend_pos::String = convert(String, legend_pos)
    if length(vectorofestimators) == 0
        error("length(vectorofestimators) == 0")
    end
    all_plots_and_legends = []
    for i = 1:length(vectorofestimators)
        estimator_i = vectorofestimators[i]
        metrics_i = singlelabelbinaryclassificationmetrics_resultdict(
            estimator_i,
            features_df,
            labels_df,
            single_label_name,
            positive_class;
            threshold = 0.5,
            )
        ytrue_i = metrics_i[:ytrue]
        yscore_i = metrics_i[:yscore]
        allfpr_i, alltpr_i, allthresholds_i = roccurve(
            ytrue_i,
            yscore_i,
            )
        plot_i = PGFPlotsX.@pgf(
            PGFPlotsX.Plot(
                PGFPlotsX.Coordinates(
                    allfpr_i,
                    alltpr_i,
                    ),
                ),
            )
        legend_i = PGFPlotsX.@pgf(
            PGFPlotsX.LegendEntry(
                LaTeXStrings.LaTeXString(estimator_i.name)
                ),
            )
        push!(all_plots_and_legends, plot_i)
        push!(all_plots_and_legends, legend_i)
    end
    all_plots_and_legends = [all_plots_and_legends...]
    p = PGFPlotsX.@pgf(
        PGFPlotsX.Axis(
            {
                xlabel = "False positive rate",
                ylabel = "True positive rate",
                no_markers,
                legend_pos = legend_pos,
                },
            all_plots_and_legends...,
            ),
        )
    wrapper = PGFPlotsXPlot(p)
    return wrapper
end

ext: jl


text:
#################################################################
#
#    1st order derivation testing for functions
#
#################################################################

tf(x) = sum(x)
@compare tf v0ref
@compare tf v1ref
@compare tf v2ref

rosenbrock(x) = (1 - x[1])^2 + 100(x[2] - x[1]^2)^2   # function to be derived
@compare rosenbrock [1., 1.]
@compare rosenbrock [0., 3.]
@compare rosenbrock [10., -10.]

function tf(x)
    a=zeros(1+3)
    for i in 1:4
        t = 4+3+2
        a[i] += b[i]*x+t
    end
    sum(a)
end
@compare tf v0ref
@compare tf 10.

function tf(x)
    a=0
    for i in 1:10
        for j in 1:10
            a += log(x) * sin(j)
        end
    end
    a
end
@compare tf 10.

function tf(x)
    a = zeros(2)
    c = 0
    for i in 1:2
        a[i] = x
        a[i] += b[i+1]
        c = c * sin(2i)
    end
    sum(a)
end
@compare tf v0ref

function tf(x)
    a = zeros(2)
    c = 0
    for i in 1:2
        d = 3
        for j in 1:2
            a[i] = x
            a[i] += b[i+1]
            d = d * sin(2i)
        end
        c += d
    end
    z = c * sum(a)
    z - sum(a)
end
@compare tf v0ref

function tf(x)
    z = 0
    for i in 1:length(x)
        z = i * x[i]
    end
    z
end
@compare tf v1ref
@compare tf v2ref


text:
import ..Tagging: hastag, traithastag, tag, mergetag
import ..Util: mergef, rmkey
import ..Traits: traits
export AbstractΩ, defΩ, defω, resolve, idtype, replacetags, like, Ω
import Random

# # Sample Space
# A sample space represents a set of possible values.
# Sample spaces are structured in the sense that that are composed of parts.
# Different parts are indicated by an id (see IDS)

"Abstract Sample Space"
abstract type AbstractΩ <: Random.AbstractRNG end

# Maybe this should be Union{AbstractΩ, AbstractRNG}
const Ω = AbstractΩ

# # Tags
"`tag(ω::AbstractΩ, tags)` tag `ω` with `tags`."
function tag(ω::AbstractΩ, tags, mergefunc=mergetag)
  replacetags(ω, mergef(mergefunc, ω.tags, tags))
end

"`rmtag(ω::AbstractΩ, tag)` remove `tag` from `ω`"
rmtag(ω::AbstractΩ, tag) =
  replacetags(ω, rmkey(ω.tags, tag))

updatetag(ω::AbstractΩ, tag, val) =
  replacetags(ω, update(ω.tags, tag, val))

traithastag(t::AbstractΩ, tag) = traithastag(t.tags, tag)
hastag(ω::AbstractΩ, tag) = hastag(ω.tags, tag)

traits(ω::AbstractΩ) = traits(ω.tags) # FIXME: Do this at the type level

# FIXME MOve this somewhere (shouldnt really be in AbstractΩ)

function replacetags end

function resolve end

"`ids(ω)` Collection of ids in `ω`, i.e. domain of ω"
function ids end

function idtype end

function like end

function subspace end

"`proj(ω, ss)` Project `ω` onto subspace `ss`"
function proj end
# Explain what project means

"`update(ω, ss)` Change subspace of  `ω` onto to `ss`"
function updatesubspace end

## Updating Interface
export update
# "`update(ω, k, v)` returns ω' which is equvalent to ω except that `ω'[k] = v`"
# function update(ω, k, v) end

# "Equivalent to setindex!, returns mutated ω"
# function update!(ω, k, v) end
# # Defaults
"Default sample space"
function defΩ end

"Default sample space object"
function defω end

"Default sample space object"
function defsubspace end

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
### Unary Map

# zero-preserving functions (z->z, nz->nz)
for op in [:-, :abs, :abs2, :conj]
    @eval begin
        $(op)(x::AbstractSparseVector) =
            SparseVector(length(x), copy(nonzeroinds(x)), $(op)(nonzeros(x)))
    end
end

# functions f, such that
#   f(x) can be zero or non-zero when x != 0
#   f(x) = 0 when x == 0
#
macro unarymap_nz2z_z2z(op, TF)
    esc(quote
        function $(op){Tv<:$(TF),Ti<:Integer}(x::AbstractSparseVector{Tv,Ti})
            R = typeof($(op)(zero(Tv)))
            xnzind = nonzeroinds(x)
            xnzval = nonzeros(x)
            m = length(xnzind)

            ynzind = Array(Ti, m)
            ynzval = Array(R, m)
            ir = 0
            @inbounds for j = 1:m
                i = xnzind[j]
                v = $(op)(xnzval[j])
                if v != zero(v)
                    ir += 1
                    ynzind[ir] = i
                    ynzval[ir] = v
                end
            end
            resize!(ynzind, ir)
            resize!(ynzval, ir)
            SparseVector(length(x), ynzind, ynzval)
        end
    end)
end

real{T<:Real}(x::AbstractSparseVector{T}) = x
@unarymap_nz2z_z2z real Complex

imag{Tv<:Real,Ti<:Integer}(x::AbstractSparseVector{Tv,Ti}) = SparseVector(length(x), Ti[], Tv[])
@unarymap_nz2z_z2z imag Complex

for op in [:floor, :ceil, :trunc, :round]
    @eval @unarymap_nz2z_z2z $(op) Real
end

for op in [:log1p, :expm1,
           :sin, :tan, :sinpi, :sind, :tand,
           :asin, :atan, :asind, :atand,
           :sinh, :tanh, :asinh, :atanh]
    @eval @unarymap_nz2z_z2z $(op) Number
end

# function that does not preserve zeros

macro unarymap_z2nz(op, TF)
    esc(quote
        function $(op){Tv<:$(TF),Ti<:Integer}(x::AbstractSparseVector{Tv,Ti})
            v0 = $(op)(zero(Tv))
            R = typeof(v0)
            xnzind = nonzeroinds(x)
            xnzval = nonzeros(x)
            n = length(x)
            m = length(xnzind)
            y = fill(v0, n)
            @inbounds for j = 1:m
                y[xnzind[j]] = $(op)(xnzval[j])
            end
            y
        end
    end)
end

for op in [:exp, :exp2, :exp10, :log, :log2, :log10,
           :cos, :csc, :cot, :sec, :cospi,
           :cosd, :cscd, :cotd, :secd,
           :acos, :acot, :acosd, :acotd,
           :cosh, :csch, :coth, :sech,
           :acsch, :asech]
    @eval @unarymap_z2nz $(op) Number
end

### Binary Map

# mode:
# 0: f(nz, nz) -> nz, f(z, nz) -> z, f(nz, z) ->  z
# 1: f(nz, nz) -> z/nz, f(z, nz) -> nz, f(nz, z) -> nz
# 2: f(nz, nz) -> z/nz, f(z, nz) -> z/nz, f(nz, z) -> z/nz

function _binarymap{Tx,Ty}(f::BinaryOp,
                           x::AbstractSparseVector{Tx},
                           y::AbstractSparseVector{Ty},
                           mode::Int)

    0 <= mode <= 2 || throw(ArgumentError("Incorrect mode $mode."))
    R = typeof(call(f, zero(Tx), zero(Ty)))
    n = length(x)
    length(y) == n || throw(DimensionMismatch())

    xnzind = nonzeroinds(x)
    xnzval = nonzeros(x)
    ynzind = nonzeroinds(y)
    ynzval = nonzeros(y)
    mx = length(xnzind)
    my = length(ynzind)
    cap = (mode == 0 ? min(mx, my) : mx + my)::Int

    rind = Array(Int, cap)
    rval = Array(R, cap)
    ir = 0
    ix = 1
    iy = 1

    ir = (
        mode == 0 ? _binarymap_mode_0!(f, mx, my,
            xnzind, xnzval, ynzind, ynzval, rind, rval) :
        mode == 1 ? _binarymap_mode_1!(f, mx, my,
            xnzind, xnzval, ynzind, ynzval, rind, rval) :
        _binarymap_mode_2!(f, mx, my,
            xnzind, xnzval, ynzind, ynzval, rind, rval)
    )::Int

    resize!(rind, ir)
    resize!(rval, ir)
    return SparseVector(n, rind, rval)
end

function _binarymap_mode_0!(f::BinaryOp, mx::Int, my::Int,
                            xnzind, xnzval, ynzind, ynzval, rind, rval)
    # f(nz, nz) -> nz, f(z, nz) -> z, f(nz, z) ->  z
    ir = 0; ix = 1; iy = 1
    @inbounds while ix <= mx && iy <= my
        jx = xnzind[ix]
        jy = ynzind[iy]
        if jx == jy
            v = call(f, xnzval[ix], ynzval[iy])
            ir += 1; rind[ir] = jx; rval[ir] = v
            ix += 1; iy += 1
        elseif jx < jy
            ix += 1
        else
            iy += 1
        end
    end
    return ir
end

function _binarymap_mode_1!{Tx,Ty}(f::BinaryOp, mx::Int, my::Int,
                                   xnzind, xnzval::AbstractVector{Tx},
                                   ynzind, ynzval::AbstractVector{Ty},
                                   rind, rval)
    # f(nz, nz) -> z/nz, f(z, nz) -> nz, f(nz, z) -> nz
    ir = 0; ix = 1; iy = 1
    @inbounds while ix <= mx && iy <= my
        jx = xnzind[ix]
        jy = ynzind[iy]
        if jx == jy
            v = call(f, xnzval[ix], ynzval[iy])
            if v != zero(v)
                ir += 1; rind[ir] = jx; rval[ir] = v
            end
            ix += 1; iy += 1
        elseif jx < jy
            v = call(f, xnzval[ix], zero(Ty))
            ir += 1; rind[ir] = jx; rval[ir] = v
            ix += 1
        else
            v = call(f, zero(Tx), ynzval[iy])
            ir += 1; rind[ir] = jy; rval[ir] = v
            iy += 1
        end
    end
    @inbounds while ix <= mx
        v = call(f, xnzval[ix], zero(Ty))
        ir += 1; rind[ir] = xnzind[ix]; rval[ir] = v
        ix += 1
    end
    @inbounds while iy <= my
        v = call(f, zero(Tx), ynzval[iy])
        ir += 1; rind[ir] = ynzind[iy]; rval[ir] = v
        iy += 1
    end
    return ir
end

function _binarymap_mode_2!{Tx,Ty}(f::BinaryOp, mx::Int, my::Int,
                                   xnzind, xnzval::AbstractVector{Tx},
                                   ynzind, ynzval::AbstractVector{Ty},
                                   rind, rval)
    # f(nz, nz) -> z/nz, f(z, nz) -> z/nz, f(nz, z) -> z/nz
    ir = 0; ix = 1; iy = 1
    @inbounds while ix <= mx && iy <= my
        jx = xnzind[ix]
        jy = ynzind[iy]
        if jx == jy
            v = call(f, xnzval[ix], ynzval[iy])
            if v != zero(v)
                ir += 1; rind[ir] = jx; rval[ir] = v
            end
            ix += 1; iy += 1
        elseif jx < jy
            v = call(f, xnzval[ix], zero(Ty))
            if v != zero(v)
                ir += 1; rind[ir] = jx; rval[ir] = v
            end
            ix += 1
        else
            v = call(f, zero(Tx), ynzval[iy])
            if v != zero(v)
                ir += 1; rind[ir] = jy; rval[ir] = v
            end
            iy += 1
        end
    end
    @inbounds while ix <= mx
        v = call(f, xnzval[ix], zero(Ty))
        if v != zero(v)
            ir += 1; rind[ir] = xnzind[ix]; rval[ir] = v
        end
        ix += 1
    end
    @inbounds while iy <= my
        v = call(f, zero(Tx), ynzval[iy])
        if v != zero(v)
            ir += 1; rind[ir] = ynzind[iy]; rval[ir] = v
        end
        iy += 1
    end
    return ir
end

function _binarymap{Tx,Ty}(f::BinaryOp,
                           x::AbstractVector{Tx},
                           y::AbstractSparseVector{Ty},
                           mode::Int)

    0 <= mode <= 2 || throw(ArgumentError("Incorrect mode $mode."))
    R = typeof(call(f, zero(Tx), zero(Ty)))
    n = length(x)
    length(y) == n || throw(DimensionMismatch())

    ynzind = nonzeroinds(y)
    ynzval = nonzeros(y)
    m = length(ynzind)

    dst = Array(R, n)
    if mode == 0
        ii = 1
        @inbounds for i = 1:m
            j = ynzind[i]
            while ii < j
                dst[ii] = zero(R); ii += 1
            end
            dst[j] = call(f, x[j], ynzval[i]); ii += 1
        end
        @inbounds while ii <= n
            dst[ii] = zero(R); ii += 1
        end
    else # mode >= 1
        ii = 1
        @inbounds for i = 1:m
            j = ynzind[i]
            while ii < j
                dst[ii] = call(f, x[ii], zero(Ty)); ii += 1
            end
            dst[j] = call(f, x[j], ynzval[i]); ii += 1
        end
        @inbounds while ii <= n
            dst[ii] = call(f, x[ii], zero(Ty)); ii += 1
        end
    end
    return dst
end

function _binarymap{Tx,Ty}(f::BinaryOp,
                           x::AbstractSparseVector{Tx},
                           y::AbstractVector{Ty},
                           mode::Int)

    0 <= mode <= 2 || throw(ArgumentError("Incorrect mode $mode."))
    R = typeof(call(f, zero(Tx), zero(Ty)))
    n = length(x)
    length(y) == n || throw(DimensionMismatch())

    xnzind = nonzeroinds(x)
    xnzval = nonzeros(x)
    m = length(xnzind)

    dst = Array(R, n)
    if mode == 0
        ii = 1
        @inbounds for i = 1:m
            j = xnzind[i]
            while ii < j
                dst[ii] = zero(R); ii += 1
            end
            dst[j] = call(f, xnzval[i], y[j]); ii += 1
        end
        @inbounds while ii <= n
            dst[ii] = zero(R); ii += 1
        end
    else # mode >= 1
        ii = 1
        @inbounds for i = 1:m
            j = xnzind[i]
            while ii < j
                dst[ii] = call(f, zero(Tx), y[ii]); ii += 1
            end
            dst[j] = call(f, xnzval[i], y[j]); ii += 1
        end
        @inbounds while ii <= n
            dst[ii] = call(f, zero(Tx), y[ii]); ii += 1
        end
    end
    return dst
end

### Binary arithmetics: +, -, *

for (vop, fun, mode) in [(:_vadd, :AddFun, 1),
                         (:_vsub, :SubFun, 1),
                         (:_vmul, :MulFun, 0)]
    @eval begin
        $(vop)(x::AbstractSparseVector, y::AbstractSparseVector) = _binarymap($(fun)(), x, y, $mode)
        $(vop)(x::StridedVector, y::AbstractSparseVector) = _binarymap($(fun)(), x, y, $mode)
        $(vop)(x::AbstractSparseVector, y::StridedVector) = _binarymap($(fun)(), x, y, $mode)
    end
end

# to workaround the ambiguities with vectorized dates/arithmetic.jl functions
if VERSION > v"0.4-dev"
    +{T<:Dates.TimeType,P<:Dates.GeneralPeriod}(x::StridedVector{P}, y::AbstractSparseVector{T}) = _vadd(x, y)
    -{T<:Dates.TimeType,P<:Dates.GeneralPeriod}(x::StridedVector{P}, y::AbstractSparseVector{T}) = _vsub(x, y)
    +{T<:Dates.TimeType,P<:Dates.GeneralPeriod}(x::AbstractSparseVector{T}, y::StridedVector{P}) = _vadd(x, y)
    -{T<:Dates.TimeType,P<:Dates.GeneralPeriod}(x::AbstractSparseVector{T}, y::StridedVector{P}) = _vsub(x, y)
end

# to workaround the ambiguities with BitVector
.*(x::BitVector, y::AbstractSparseVector{Bool}) = _vmul(x, y)
.*(x::AbstractSparseVector{Bool}, y::BitVector) = _vmul(x, y)

# definition of operators

for (op, vop) in [(:+, :_vadd), (:(.+), :_vadd),
                  (:-, :_vsub), (:(.-), :_vsub),
                  (:.*, :_vmul)]
    @eval begin
        $(op)(x::AbstractSparseVector, y::AbstractSparseVector) = $(vop)(x, y)
        $(op)(x::StridedVector, y::AbstractSparseVector) = $(vop)(x, y)
        $(op)(x::AbstractSparseVector, y::StridedVector) = $(vop)(x, y)
    end
end

# definition of other binary functions

for (op, fun, TF, mode) in [(:max, :MaxFun, :Real, 2),
                            (:min, :MinFun, :Real, 2),
                            (:complex, :ComplexFun, :Real, 1)]
    @eval begin
        $(op){Tx<:$(TF),Ty<:$(TF)}(x::AbstractSparseVector{Tx}, y::AbstractSparseVector{Ty}) =
            _binarymap($(fun)(), x, y, $mode)
        $(op){Tx<:$(TF),Ty<:$(TF)}(x::StridedVector{Tx}, y::AbstractSparseVector{Ty}) =
            _binarymap($(fun)(), x, y, $mode)
        $(op){Tx<:$(TF),Ty<:$(TF)}(x::AbstractSparseVector{Tx}, y::StridedVector{Ty}) =
            _binarymap($(fun)(), x, y, $mode)
    end
end

### Reduction

sum(x::AbstractSparseVector) = sum(nonzeros(x))
sumabs(x::AbstractSparseVector) = sumabs(nonzeros(x))
sumabs2(x::AbstractSparseVector) = sumabs2(nonzeros(x))

function maximum{T<:Real}(x::AbstractSparseVector{T})
    n = length(x)
    n > 0 || throw(ArgumentError("maximum over empty array is not allowed."))
    m = nnz(x)
    (m == 0 ? zero(T) :
     m == n ? maximum(nonzeros(x)) :
     max(zero(T), maximum(nonzeros(x))))::T
end

function minimum{T<:Real}(x::AbstractSparseVector{T})
    n = length(x)
    n > 0 || throw(ArgumentError("minimum over empty array is not allowed."))
    m = nnz(x)
    (m == 0 ? zero(T) :
     m == n ? minimum(nonzeros(x)) :
     min(zero(T), minimum(nonzeros(x))))::T
end

maxabs{T<:Number}(x::AbstractSparseVector{T}) = maxabs(nonzeros(x))
minabs{T<:Number}(x::AbstractSparseVector{T}) = nnz(x) < length(x) ? abs(zero(T)) : minabs(nonzeros(x))

vecnorm(x::AbstractSparseVector, p::Real=2) = vecnorm(nonzeros(x), p)


text:
# test cube
@testset "Primitives" begin
    c1 = Cuboid([1,1,1])
    c2 = Cuboid([1,1,1])
    c3 = Cuboid([0,1,1])
    #@test isequal(c1, c2)
    #@test !isequal(c1,c3)
end


text:
function _svg(h::Hist1D)
    paddingx, paddingy = 0.05, 0.10
    framewidth, frameheight = 250, 200
    fill = "#ffffff00" # 0 alpha
    strokecolor, strokewidth = "black", 1
    strokewidth = 1
    _c, _e = copy(bincounts(h)), copy(binedges(h))

    hasneg = minimum(_c) < 0
    zerolineyfrac = 0.0
    if hasneg
        zerolineyfrac = -minimum(_c)/(maximum(_c)-minimum(_c))
        _c .-= minimum(_c)
    end
    ys = frameheight * ((2*paddingy-1)/maximum(_c) .* _c .+ (1-paddingy))
    xs = framewidth * (
        (1 - 2 * paddingx) / (maximum(_e) - minimum(_e))
        .* (_e .- minimum(_e))
        .+ paddingx
    )
    points = [(paddingx*framewidth, (1-paddingy)*frameheight)] # bottom left
    for i in 1:(length(xs)-1)
        push!(points, (xs[i],ys[i])) # left bin edge
        push!(points, (xs[i+1],ys[i])) # right bin edge
    end
    push!(points, ((1-paddingx)*framewidth,(1-paddingy)*frameheight))
    if !hasneg
        push!(points, points[1]) # close path
    end
    pathstr = join(["$(x),$(y)" for (x,y) in points],",")
    zeroliney = hasneg ? frameheight*(1-zerolineyfrac) : frameheight*(1-paddingy)
    xlow, xhigh = round.(extrema(_e), sigdigits=10)
    return """
    <svg width="$(framewidth)" height="$(frameheight)" version="1.1" xmlns="http://www.w3.org/2000/svg">
        <polyline points="$(pathstr)" stroke="$(strokecolor)" fill="$(fill)" stroke-width="$(strokewidth)"/>
        <polyline points="$(framewidth*paddingx),$(zeroliney),$(framewidth*(1-paddingx)),$(zeroliney)" stroke="black" stroke-width="1"/>
        <text x="$(framewidth*paddingx)" y="$(frameheight*(1-0.5*paddingy))" dominant-baseline="middle" text-anchor="start" fill="black">$(xlow)</text>
        <text x="$(framewidth*(1-paddingx))" y="$(frameheight*(1-0.5*paddingy))" dominant-baseline="middle" text-anchor="end" fill="black">$(xhigh)</text>
    </svg>
    """
end

function _svg(h::Hist2D)
    paddingx1, paddingy1 = 0.15, 0.15 # left, bottom
    paddingx2, paddingy2 = 0.05, 0.05 # right, top
    framewidth, frameheight = 250, 200

    (xlow, xhigh), (ylow, yhigh) = extrema.(binedges(h))
    function transform(x::Real, y::Real)
        xfrac = (x - xlow)/(xhigh - xlow)
        xdraw = xfrac * framewidth*(1-paddingx1-paddingx2) + framewidth*paddingx1
        yfrac = (y - ylow)/(yhigh - ylow)
        ydraw = frameheight - (yfrac * frameheight*(1-paddingy1-paddingy2) + frameheight*paddingy1)
        return xdraw, ydraw
    end

    function colorscale(x::Real)
        x = isnan(x) ? 0. : x
        # x is fraction between [0,1]
        # 6th deg polynomial fit to viridis color palette
        # See https://gist.github.com/aminnj/dba84777718613d3a37291a43659feff
        # to generate other color palettes. Black and white is trivial:
        #   return 255 .* (1-x,1-x,1-x)
        r = @evalpoly(x, 0.2731, 0.1270, -0.3617, -4.7456, 6.7092, 4.2491, -5.2678)
        g = @evalpoly(x, 0.0039, 1.3810, 0.3969, -6.4246, 15.3241, -14.7345, 4.9587)
        b = @evalpoly(x, 0.3305, 1.3726, 0.3948, -20.7431, 59.7287, -68.3565, 27.4051)
        return round.(Int, 255 .* (r,g,b))
    end

    counts = bincounts(h)
    mincount, maxcount = extrema(counts)
    counts = clamp.(counts, 0, maxcount)
    ex, ey = binedges(h)
    sx, sy = size(counts)
    rectlines = String[]
    textstyle = """ dominant-baseline="middle" text-anchor="middle" font-size="85%" font-family="sans-serif" """
    for i in 1:sx, j in 1:sy
        tx1, ty1 = ceil.(Int, transform(ex[i], ey[j+1]))
        tx2, ty2 = ceil.(Int, transform(ex[i+1], ey[j]))
        tw, th = tx2-tx1, ty2-ty1
        c = counts[i,j]
        (c == 0) && continue
        r,g,b = colorscale((c-mincount)/(maxcount-mincount))
        rcolor = "rgb($(r),$(g),$(b))"
        line = """<rect x="$(tx1)" y="$(ty1)" width="$(tw)" height="$(th)" fill="$(rcolor)" stroke="none" />"""
        if (sx <= 15) && (sy <= 15)
            tcolor = ((0.299*r + 0.587*g + 0.114*b) < 0.60*255) ? "#fff" : "#000"
            line = """<g>$(line)<text class="svgplotlabels" x="$(tx1+tw/2)" y="$(ty1+th/2)" fill="$(tcolor)" $(textstyle)>$(c)</text></g>"""
        end
        push!(rectlines, line * "\n")
    end

    return """
    <svg width="$(framewidth)" height="$(frameheight)" version="1.1" xmlns="http://www.w3.org/2000/svg" class="svgplot">

        $(join(rectlines))
        <rect x="$(round(Int,framewidth*paddingx1))" y="$(round(Int,frameheight*paddingy2))" width="$(round(Int,framewidth*(1-paddingx1-paddingx2)))" height="$(frameheight*(1-paddingy1-paddingy2))" fill="none" stroke="#000" stroke-width="1" />
        <text x="$(framewidth*paddingx1)" y="$(frameheight*(1-0.5*paddingy1))" fill="black" $(textstyle)>$(minimum(ex))</text>
        <text x="$(framewidth*(1-paddingx2))" y="$(frameheight*(1-0.5*paddingy1))" fill="black" $(textstyle)>$(maximum(ex))</text>
        <text x="$(framewidth*0.5*paddingx1)" y="$(frameheight*(1-paddingy1))" fill="black" $(textstyle)>$(minimum(ey))</text>
        <text x="$(framewidth*0.5*paddingx1)" y="$(frameheight*paddingy2)" fill="black" $(textstyle)>$(maximum(ey))</text>
    </svg>
    """
end

function _svg(h::Hist3D)
    framewidth = 250
    frameheight = 200

    x1, x2 = 60, 180
    y1, y2 = 30, 130
    z1, z2 = 10, 120

    # cabinet projection
    α = pi/4
    P = [1 0 1/2*cos(α) ; 0 1 1/2*sin(α); 0 0 0]
    projection(points) = map(p->begin
                                 x,y,z = (P*p)
                                 [round(x),round(frameheight-y)]
                             end, points)
    flatstring(points) = join([points...;], " ")

    frontface = projection([[x1,y1,z1],[x2,y1,z1],[x2,y2,z1],[x1,y2,z1],[x1,y1,z1]])
    rightface = projection([[x2,y1,z2],[x2,y2,z2],[x2,y2,z1],[x2,y1,z1],[x2,y1,z2]])
    topface = projection([[x1,y2,z2],[x1,y2,z1],[x2,y2,z1],[x2,y2,z2],[x1,y2,z2]])
    backface = projection([[x1,y1,z2],[x2,y1,z2],[x2,y2,z2],[x1,y2,z2],[x1,y1,z2]])
    bottomface = projection([[x1,y1,z2],[x1,y1,z1],[x2,y1,z1],[x2,y1,z2],[x1,y1,z2]])
    xaxisline = projection([[x1,y1,z1],[x2,y1,z1]])
    yaxisline = projection([[x1,y1,z1],[x1,y2,z1]])
    zaxisline = projection([[x1,y2,z1],[x1,y2,z2]])

    (xlow, xhigh), (ylow, yhigh), (zlow, zhigh) = extrema.(binedges(h))
    function data_to_pixels(x, y, z)
        xt = (x - xlow)/(xhigh - xlow)*(x2-x1)+x1
        yt = (y - ylow)/(yhigh - ylow)*(y2-y1)+y1
        zt = (z - zlow)/(zhigh - zlow)*(z2-z1)+z1
        return xt, yt, zt
    end

    boxes = ""
    if prod(nbins(h)) < 5000
        counts = bincounts(h)
        maxcount = maximum(counts)
        counts = clamp.(counts, 0, maxcount)
        ex, ey, ez = binedges(h)
        sx, sy, sz = size(counts)
        lines = String[]
        for i in 1:sx, j in 1:sy, k in 1:sz
            c = counts[i,j,k]
            (c == 0) && continue
            x1t,y1t,z1t = data_to_pixels(ex[i], ey[j], ez[k])
            x2t,y2t,z2t = data_to_pixels(ex[i+1], ey[j+1], ez[k+1])
            sf = c/maxcount # normalize box size by count relative to maximum
            zf = 1 - (0.5*(ez[k]+ez[k+1])-zlow) / (zhigh-zlow)
            opacity = 0.3 + 0.3*zf # more transparent as you go toward the rear
            dx, dy, dz = x2t - x1t, y2t - y1t, z2t - z1t
            x1t += dx*(1-sf)/2
            x2t -= dx*(1-sf)/2
            y1t += dy*(1-sf)/2
            y2t -= dy*(1-sf)/2
            z1t += dz*(1-sf)/2
            z2t -= dz*(1-sf)/2
            face1 = projection([[x1t,y1t,z1t],[x2t,y1t,z1t],[x2t,y2t,z1t],[x1t,y2t,z1t],[x1t,y1t,z1t]])
            face2 = projection([[x2t,y1t,z2t],[x2t,y2t,z2t],[x2t,y2t,z1t],[x2t,y1t,z1t],[x2t,y1t,z2t]])
            face3 = projection([[x1t,y2t,z2t],[x1t,y2t,z1t],[x2t,y2t,z1t],[x2t,y2t,z2t],[x1t,y2t,z2t]])
            push!(lines, """<polyline points="$(flatstring(face1))" stroke="none" fill="#888" opacity="$(opacity)"/>\n""")
            push!(lines, """<polyline points="$(flatstring(face2))" stroke="none" fill="#ccc" opacity="$(opacity)"/>\n""")
            push!(lines, """<polyline points="$(flatstring(face3))" stroke="none" fill="#aaa" opacity="$(opacity)"/>\n""")
        end
        boxes = join(lines)
    end

    c1, c2, c3 = "#1f77b4", "#d62728", "#2ca02c"
    (x1text, x2text), (y1text, y2text), (z1text, z2text) = map(x->(first(x),last(x)), binedges(h))
    textstyle = """ font-size="85%" font-family="sans-serif" """
    return """
    <svg width="$(framewidth)" height="$(frameheight)" version="1.1" xmlns="http://www.w3.org/2000/svg">
        <polyline points="$(flatstring(backface))" stroke="gray" stroke-width="1" stroke-dasharray="4" fill="none"/>
        <polyline points="$(flatstring(bottomface))" stroke="gray" stroke-width="1" stroke-dasharray="4" fill="none"/>
        $(boxes)
        <polyline points="$(flatstring(frontface))" stroke="black" stroke-width="1" fill="none"/>
        <polyline points="$(flatstring(rightface))" stroke="black" stroke-width="1" fill="none"/>
        <polyline points="$(flatstring(topface))" stroke="black" stroke-width="1" fill="none"/>

        <polyline points="$(flatstring(xaxisline))" stroke="$(c1)" stroke-width="2" fill="none"/>
        <polyline points="$(flatstring(yaxisline))" stroke="$(c2)" stroke-width="2" fill="none"/>
        <polyline points="$(flatstring(zaxisline))" stroke="$(c3)" stroke-width="2" fill="none"/>

        <text x="$(mean(xaxisline[1])+10)" y="$(xaxisline[1][2]+5)" fill="$(c1)" dominant-baseline="hanging" text-anchor="middle" $(textstyle)>x</text>
        <text x="$(xaxisline[1][1])" y="$(xaxisline[1][2]+5)" fill="$(c1)" dominant-baseline="hanging" text-anchor="left" $(textstyle)>$(x1text)</text>
        <text x="$(xaxisline[2][1])" y="$(xaxisline[2][2]+5)" fill="$(c1)" dominant-baseline="hanging" text-anchor="end" $(textstyle)>$(x2text)</text>

        <text x="$(yaxisline[1][1]-5)" y="$(mean(yaxisline[1]))" fill="$(c2)" dominant-baseline="text-bottom" text-anchor="end" $(textstyle)>y</text>
        <text x="$(yaxisline[1][1]-5)" y="$(yaxisline[1][2])" fill="$(c2)" dominant-baseline="text-bottom" text-anchor="end" $(textstyle)>$(y1text)</text>
        <text x="$(yaxisline[2][1]-5)" y="$(yaxisline[2][2])" fill="$(c2)" dominant-baseline="hanging" text-anchor="end" $(textstyle)>$(y2text)</text>

        <text x="$((zaxisline[1][1]+zaxisline[2][1])/2-10)" y="$((zaxisline[1][2]+zaxisline[2][2])/2)" fill="$(c3)" dominant-baseline="text-bottom" text-anchor="end" $(textstyle)>z</text>
        <text x="$(zaxisline[1][1]+15)" y="$(zaxisline[1][2]-4)" fill="$(c3)" dominant-baseline="text-bottom" text-anchor="start" $(textstyle)>$(z1text)</text>
        <text x="$(zaxisline[2][1]+5)" y="$(zaxisline[2][2]+3)" fill="$(c3)" dominant-baseline="hanging" text-anchor="start" $(textstyle)>$(z2text)</text>
    </svg>
         """
end

function Base.show(io::IO, h::Hist1D)
    compact = get(io, :compact, false)
    if compact
        print(io, "Hist1D{$(eltype(bincounts(h)))}, ")
        print(io, "edges=$(repr(binedges(h), context=:limit => true)), ")
        print(io, "integral=$(integral(h))")
    else
        if (nbins(h) < 50) && all(bincounts(h) .>= 0)
            _e = binedges(h)
            _h = Histogram(float.(_e), bincounts(h))
            show(io, UnicodePlots.histogram(_h; width=30, xlabel=""))
        end
        println(io)
        println(io, "edges: ", binedges(h))
        println(io, "bin counts: ", bincounts(h))
        print(io, "total count: ", integral(h))
    end
end

function Base.show(io::IO, m::MIME"text/html", h::Hist1D)
    println(io, """

$(_svg(h))

edges: $(repr(binedges(h), context=:limit => true))

bin counts: $(repr(bincounts(h), context=:limit => true))

maximum count: $(maximum(bincounts(h)))

total count: $(integral(h))

    """)
end

function Base.show(io::IO, h::Hist2D)
    compact = get(io, :compact, false)
    if compact
        print(io, "Hist2D{$(eltype(bincounts(h)))}, ")
        print(io, "edges=$(repr(binedges(h), context=:limit => true)), ")
        print(io, "integral=$(integral(h))")
    else
        ex, ey = binedges(h)
        nx, ny = nbins(h)
        xfact = nx > 1 ? (maximum(ex)-minimum(ex))/(nx-1) : 0.0
        yfact = ny > 1 ? (maximum(ey)-minimum(ey))/(ny-1) : 0.0
        show(io, UnicodePlots.heatmap(bincounts(h)'; xfact=xfact, xoffset=minimum(ex), yfact=yfact, yoffset=minimum(ey)))
        println(io)
        println(io, "edges: ", binedges(h))
        println(io, "bin counts: ", bincounts(h))
        print(io, "total count: ", integral(h))
    end
end

function Base.show(io::IO, m::MIME"text/html", h::Hist2D)
    println(io, """

$(_svg(h))

edges: $(repr(binedges(h), context=:limit => true))

bin counts: $(repr(bincounts(h), context=:limit => true))

maximum count: $(maximum(bincounts(h)))

total count: $(integral(h))

    """)
end

function Base.show(io::IO, h::Hist3D)
    print(io, "Hist3D{$(eltype(bincounts(h)))}, ")
    print(io, "edges=$(repr(binedges(h), context=:limit => true)), ")
    print(io, "integral=$(integral(h))")
end

function Base.show(io::IO, m::MIME"text/html", h::Hist3D)
    println(io, """

$(_svg(h))

edges: $(repr(binedges(h), context=:limit => true))

bin counts: $(repr(bincounts(h), context=:limit => true))

maximum count: $(maximum(bincounts(h)))

total count: $(integral(h))

    """)
end


text:
include("base_data.jl")

include("stream_merger.jl")
include("poisson_stream.jl")
include("regular_pattern_stream.jl")
include("pattern_channel.jl")

include("streams.jl")


text:
"""
`Manopt.jl` – Optimization on Manifolds in Julia.
"""
module Manopt
  using SimpleTraits
  using Markdown
# Manifolds
  include("manifolds/Manifold.jl") #base type
  # errors for false combinations of types or nonimplemented cases
  include("manifolds/defaults/manifoldFallbacks.jl")
  # Extended Vector decorations
  include("manifolds/defaults/extendedData.jl")

  # Traits (properties/decorators)
  include("manifolds/traits/EmbeddedManifold.jl")
  include("manifolds/traits/LieGroup.jl")
  include("manifolds/traits/MatrixManifold.jl")
  # specific manifolds
  include("manifolds/Circle.jl")
  include("manifolds/Euclidean.jl")
  include("manifolds/Graph.jl")
  include("manifolds/Grassmannian.jl")
  include("manifolds/Hyperbolic.jl")
  include("manifolds/Rotations.jl")
  include("manifolds/Sphere.jl")
  include("manifolds/Stiefel.jl")
  include("manifolds/SymmetricPositiveDefinite.jl")
  include("manifolds/Symmetric.jl")
  include("manifolds/TangentBundle.jl")
  # meta
  include("manifolds/Power.jl")
  include("manifolds/Product.jl")
  # ...corresponding plans consisting of problems and options
  include("plans/plan.jl")
  # Functions
  include("functions/adjointDifferentials.jl")
  include("functions/costFunctions.jl")
  include("functions/differentials.jl")
  include("functions/gradients.jl")
  include("functions/jacobiFields.jl")
  include("functions/proximalMaps.jl")
  # ...solvers (1) general framework
  include("solvers/solver.jl")
  # ...solvers (2) specific solvers
  include("solvers/cyclicProximalPoint.jl")
  include("solvers/DouglasRachford.jl")
  include("solvers/steepestDescent.jl")
  include("solvers/subGradientMethod.jl")
  # extended metasolvers
  include("solvers/debugSolver.jl")
  include("solvers/recordSolver.jl")
  # algorithms
  include("algorithms/basicAlgorithms.jl")
  # Plots
  include("plots/SpherePlots.jl")
  # helpers
  include("helpers/errorMeasures.jl")
  include("helpers/imageHelpers.jl")
  # Exports
  include("helpers/exports/Asymptote.jl")
  # data
  include("data/artificialDataFunctions.jl")
end


text:
module Flist                    # feature matrices
using KUparser: SFeature, DFeature # sparse and dense features

gn13cpv = DFeature["s0v","s1v","n0v","n1v",
	           "s0c","s1c","n0c","n1c","s1rc","s1lc","s0lc","s0rc","n0lc",
	           "s0p","s1p","n0p","n1p","s1rp","s1lp","s0lp","s0rp","n0lp"]

# ZN11 dense variants.  The original has:
# (w,p) for (s,n,n1,n2,sh,sh2,sl,sr,nl,sl2,sr2,nl2): 24 features, 1200 dims for v,c; 540 dims for p
# L for (s,sh,sl,sr,nl,sl2,sr2,nl2): 8 features, 96 dims
# sd, sa, sb, sA, sB, na, nA: 7 features, 69 dims
# We replace (w,p) with 8 subsets of (v,c,p)

zn11n = DFeature["sA","nA","sa","na","sB","sb","sd",
                 "sL","shL","slL","srL","nlL","sl2L","sr2L","nl2L"]

zn11p = DFeature["sA","nA","sa","na","sB","sb","sd",
                 "sL","shL","slL","srL","nlL","sl2L","sr2L","nl2L",
                 "sp","np","n1p","n2p","shp","sh2p","slp","srp","nlp","sl2p","sr2p","nl2p"]

zn11v = DFeature["sA","nA","sa","na","sB","sb","sd",
                 "sL","shL","slL","srL","nlL","sl2L","sr2L","nl2L",
                 "sv","nv","n1v","n2v","shv","sh2v","slv","srv","nlv","sl2v","sr2v","nl2v"]

zn11c = DFeature["sA","nA","sa","na","sB","sb","sd",
                 "sL","shL","slL","srL","nlL","sl2L","sr2L","nl2L",
                 "sc","nc","n1c","n2c","shc","sh2c","slc","src","nlc","sl2c","sr2c","nl2c"]

zn11pv = DFeature["sA","nA","sa","na","sB","sb","sd",
                  "sL","shL","slL","srL","nlL","sl2L","sr2L","nl2L",
                  "sp","np","n1p","n2p","shp","sh2p","slp","srp","nlp","sl2p","sr2p","nl2p",
                  "sv","nv","n1v","n2v","shv","sh2v","slv","srv","nlv","sl2v","sr2v","nl2v"]

zn11cv = DFeature["sA","nA","sa","na","sB","sb","sd",
                  "sL","shL","slL","srL","nlL","sl2L","sr2L","nl2L",
                  "sc","nc","n1c","n2c","shc","sh2c","slc","src","nlc","sl2c","sr2c","nl2c",
                  "sv","nv","n1v","n2v","shv","sh2v","slv","srv","nlv","sl2v","sr2v","nl2v"]

zn11cp = DFeature["sA","nA","sa","na","sB","sb","sd",
                  "sL","shL","slL","srL","nlL","sl2L","sr2L","nl2L",
                  "sc","nc","n1c","n2c","shc","sh2c","slc","src","nlc","sl2c","sr2c","nl2c",
                  "sp","np","n1p","n2p","shp","sh2p","slp","srp","nlp","sl2p","sr2p","nl2p"]

zn11cpv = DFeature["sA","nA","sa","na","sB","sb","sd",
                   "sL","shL","slL","srL","nlL","sl2L","sr2L","nl2L",
                   "sc","nc","n1c","n2c","shc","sh2c","slc","src","nlc","sl2c","sr2c","nl2c",
                   "sp","np","n1p","n2p","shp","sh2p","slp","srp","nlp","sl2p","sr2p","nl2p",
                   "sv","nv","n1v","n2v","shv","sh2v","slv","srv","nlv","sl2v","sr2v","nl2v"]

# Sparse feature sets
#
# Each string specifies a particular feature and has the form:
#   [sn]\d?([hlr]\d?)*[wpdLabAB]
#
# The meaning of each letter is below.  In the following i is a single
# digit integer which is optional if the default value is used:
#
# si: i'th stack word, default i=0 means top
# ni: i'th buffer word, default i=0 means first
# hi: i'th degree head, default i=1 means direct head
# li: i'th leftmost child, default i=1 means the leftmost child
# ri: i'th rightmost child, default i=1 means the rightmost child
# w: word
# p: postag
# d: distance to the right.  e.g. s1d is s0s1 distance, s0d is s0n0
#    distance.  encoding: 1,2,3,4,5-10,10+ (from ZN11)
# L: dependency label (0 is ROOT or NONE)
# a: number of left children.  following ZN11, any number is allowed.
# b: number of right children.  following ZN11, any number is allowed.
# A: set of left dependency labels
# B: set of right dependency labels

# The original zn11 set in sparse format with feature conjunctions
zn11orig = SFeature[
["sw"],                  # STw		StackWord
["sp"],                  # STt		StackTag
["sw","sp"],             # STwt		StackWordTag
["nw"],                  # N0w		NextWord
["np"],                  # N0t		NextTag
["nw","np"],             # N0wt		NextWordTag
["n1w"],                 # N1w		Next+1Word
["n1p"],                 # N1t		Next+1Tag
["n1w","n1p"],           # N1wt		Next+1WordTag
["n2w"],                 # N2w		Next+2Word
["n2p"],                 # N2t		Next+2Tag
["n2w","n2p"],           # N2wt		Next+2WordTag
["shw"],                 # STHw		StackHeadWord
["shp"],                 # STHt		StackHeadTag
["sL"],                  # STi		StackLabel
["sh2w"],                # STHHw		StackHeadHeadWord
["sh2p"],                # STHHt		StackHeadHeadTag
["shL"],                 # STHi		StackLabel ?? StackHeadLabel
["slw"],                 # STLDw		StackLDWord
["slp"],                 # STLDt		StackLDTag
["slL"],                 # STLDi		StackLDLabel
["srw"],                 # STRDw		StackRDWord
["srp"],                 # STRDt		StackRDTag
["srL"],                 # STRDi		StackRDLabel
["nlw"],                 # N0LDw		NextLDWord
["nlp"],                 # N0LDt		NextLDTag
["nlL"],                 # N0LDi		NextLDLabel
["sl2w"],                # STL2Dw	StackL2DWord
["sl2p"],                # STL2Dt	StackL2DTag
["sl2L"],                # STL2Di	StackL2DLabel
["sr2w"],                # STR2Dw	StackR2DWord
["sr2p"],                # STR2Dt	StackR2DTag
["sr2L"],                # STR2Di	StackR2DLabel
["nl2w"],                # N0L2Dw	NextL2DWord
["nl2p"],                # N0L2Dt	NextL2DTag
["nl2L"],                # N0L2Di	NextL2DLabel
			## HTw		HeadStackWord
			## HTt		HeadStackTag
			## HTwt		HeadStackWordTag
["sw","sp","nw","np"],	# STwtN0wt	StackWordTagNextWordTag
["sw","sp","nw"],	# STwtN0w	StackWordTagNextWord
["sw","nw","np"],	# STwN0wt	StackWordNextWordTag
["sp","nw","np"],	# STtN0wt	StackTagNextWordTag
["sw","sp","np"],	# STwtN0t	StackWordTagNextTag
["sw","nw"],             # STwN0w	StackWordNextWord
["sp","np"],             # STtN0t	StackTagNextTag
["np","n1p"],            # N0tN1t	NextTagNext+1Tag
["np","n1p","n2p"],	# N0tN1tN2t	NextTagTrigram
["sp","np","n1p"],	# STtN0tN1t	StackTagNextTagNext+1Tag
["sp","np","nlp"],	# STtN0tN0LDt	StackTagNextTagNextLDTag
["np","nlp","nl2p"],     # N0tN0LDtN0L2Dt	StackTagNextTagNextLDTagNextTagNextL2DTag ?? NextTagNextLDTagNextL2DTag
["shp","sp","np"],       # STHtSTtN0t	StackHeadTagStackTagNextTag
			## HTtHT2tN0t	HeadStackTagHeadStack2TagNextTag
["sh2p","shp","sp"],     # STHHtSTHtSTt	StackHeadHeadTagStackHeadTagStackTag
["sp","slp","np"],       # STtSTLDtN0t	StackTagStackLDTagNextTag
["sp","slp","sl2p"],     # STtSTLDtSTL2Dt	StackTagStackLDTagStackL2DTag
["sp","srp","np"],       # STtSTRDtN0t	StackTagStackRDTagNextTag
["sp","srp","sr2p"],     # STtSTRDtSTR2Dt	StackTagStackRDTagStackR2DTag
["sw","sd"],             # STwd		StackWordDist
["sp","sd"],             # STtd		StackTagDist
["nw","sd"],             # N0wd		NextWordDist
["np","sd"],             # N0td		NextTagDist
["sw","nw","sd"],        # STwN0wd	StackWordNextWordDist
["sp","np","sd"],        # STtN0td	StackTagNextTagDist
["sw","sb"],             # STwra		StackWordRightArity
["sp","sb"],             # STtra		StackTagRightArity
["sw","sa"],             # STwla		StackWordLeftArity
["sp","sa"],             # STtla		StackTagLeftArity
["nw","na"],             # N0wla		NextWordRightArity ?? NextWordLeftArity
["np","na"],             # N0tla		NextTagRightArity  ?? NextTagLeftArity
["sw","sB"],             # STwrp		StackWordRightSetoftags
["sp","sB"],             # STtrp		StackTagRightSetoftags
["sw","sA"],             # STwlp		StackWordLeftSetoftags
["sp","sA"],             # STtlp		StackTagLeftSetoftags
["nw","nA"],             # N0wlp		Next0WordLeftSetoftags ?? NextWordLeftSetoftags
["np","nA"],             # N0tlp		Next0TagLeftSetoftags  ?? NextTagLeftSetoftags
			## STl		StackLemma
			## STc		StackCPOS
			## STf		StackFeats
			## N0l		NextLemma
			## N0c		NextCPOS
			## N0f		NextFeats
			## N1l		Next+1Lemma
			## N1c		Next+1CPOS
			## N1f		Next+1Feats
]

# The reduced zn11 set in sparse format with singleton features
zn11single = SFeature[
["sw"],                  # STw		StackWord
["sp"],                  # STt		StackTag
#["sw","sp"],             # STwt		StackWordTag
["nw"],                  # N0w		NextWord
["np"],                  # N0t		NextTag
#["nw","np"],             # N0wt		NextWordTag
["n1w"],                 # N1w		Next+1Word
["n1p"],                 # N1t		Next+1Tag
#["n1w","n1p"],           # N1wt		Next+1WordTag
["n2w"],                 # N2w		Next+2Word
["n2p"],                 # N2t		Next+2Tag
#["n2w","n2p"],           # N2wt		Next+2WordTag
["shw"],                 # STHw		StackHeadWord
["shp"],                 # STHt		StackHeadTag
["sL"],                  # STi		StackLabel
["sh2w"],                # STHHw		StackHeadHeadWord
["sh2p"],                # STHHt		StackHeadHeadTag
["shL"],                 # STHi		StackLabel ?? StackHeadLabel
["slw"],                 # STLDw		StackLDWord
["slp"],                 # STLDt		StackLDTag
["slL"],                 # STLDi		StackLDLabel
["srw"],                 # STRDw		StackRDWord
["srp"],                 # STRDt		StackRDTag
["srL"],                 # STRDi		StackRDLabel
["nlw"],                 # N0LDw		NextLDWord
["nlp"],                 # N0LDt		NextLDTag
["nlL"],                 # N0LDi		NextLDLabel
["sl2w"],                # STL2Dw	StackL2DWord
["sl2p"],                # STL2Dt	StackL2DTag
["sl2L"],                # STL2Di	StackL2DLabel
["sr2w"],                # STR2Dw	StackR2DWord
["sr2p"],                # STR2Dt	StackR2DTag
["sr2L"],                # STR2Di	StackR2DLabel
["nl2w"],                # N0L2Dw	NextL2DWord
["nl2p"],                # N0L2Dt	NextL2DTag
["nl2L"],                # N0L2Di	NextL2DLabel
			## HTw		HeadStackWord
			## HTt		HeadStackTag
			## HTwt		HeadStackWordTag
#["sw","sp","nw","np"],	# STwtN0wt	StackWordTagNextWordTag
#["sw","sp","nw"],	# STwtN0w	StackWordTagNextWord
#["sw","nw","np"],	# STwN0wt	StackWordNextWordTag
#["sp","nw","np"],	# STtN0wt	StackTagNextWordTag
#["sw","sp","np"],	# STwtN0t	StackWordTagNextTag
#["sw","nw"],             # STwN0w	StackWordNextWord
#["sp","np"],             # STtN0t	StackTagNextTag
#["np","n1p"],            # N0tN1t	NextTagNext+1Tag
#["np","n1p","n2p"],	# N0tN1tN2t	NextTagTrigram
#["sp","np","n1p"],	# STtN0tN1t	StackTagNextTagNext+1Tag
#["sp","np","nlp"],	# STtN0tN0LDt	StackTagNextTagNextLDTag
#["np","nlp","nl2p"],     # N0tN0LDtN0L2Dt	StackTagNextTagNextLDTagNextTagNextL2DTag ?? NextTagNextLDTagNextL2DTag
#["shp","sp","np"],       # STHtSTtN0t	StackHeadTagStackTagNextTag
			## HTtHT2tN0t	HeadStackTagHeadStack2TagNextTag
#["sh2p","shp","sp"],     # STHHtSTHtSTt	StackHeadHeadTagStackHeadTagStackTag
#["sp","slp","np"],       # STtSTLDtN0t	StackTagStackLDTagNextTag
#["sp","slp","sl2p"],     # STtSTLDtSTL2Dt	StackTagStackLDTagStackL2DTag
#["sp","srp","np"],       # STtSTRDtN0t	StackTagStackRDTagNextTag
#["sp","srp","sr2p"],     # STtSTRDtSTR2Dt	StackTagStackRDTagStackR2DTag
["sd"],
#["sw","sd"],             # STwd		StackWordDist
#["sp","sd"],             # STtd		StackTagDist
#["nw","sd"],             # N0wd		NextWordDist
#["np","sd"],             # N0td		NextTagDist
#["sw","nw","sd"],        # STwN0wd	StackWordNextWordDist
#["sp","np","sd"],        # STtN0td	StackTagNextTagDist
["sb"],
#["sw","sb"],             # STwra		StackWordRightArity
#["sp","sb"],             # STtra		StackTagRightArity
["sa"],
#["sw","sa"],             # STwla		StackWordLeftArity
#["sp","sa"],             # STtla		StackTagLeftArity
["na"],
#["nw","na"],             # N0wla		NextWordRightArity ?? NextWordLeftArity
#["np","na"],             # N0tla		NextTagRightArity  ?? NextTagLeftArity
["sB"],
#["sw","sB"],             # STwrp		StackWordRightSetoftags
#["sp","sB"],             # STtrp		StackTagRightSetoftags
["sA"],
#["sw","sA"],             # STwlp		StackWordLeftSetoftags
#["sp","sA"],             # STtlp		StackTagLeftSetoftags
["nA"],
#["nw","nA"],             # N0wlp		Next0WordLeftSetoftags ?? NextWordLeftSetoftags
#["np","nA"],             # N0tlp		Next0TagLeftSetoftags  ?? NextTagLeftSetoftags
			## STl		StackLemma
			## STc		StackCPOS
			## STf		StackFeats
			## N0l		NextLemma
			## N0c		NextCPOS
			## N0f		NextFeats
			## N1l		Next+1Lemma
			## N1c		Next+1CPOS
			## N1f		Next+1Feats
]

# Dense feature sets:

# Each string specifies a particular feature and has the form:
#   [sn]\d?([hlr]\d?)*[vcpdLabAB]
#
# The meaning of each letter is below.  In the following i is a single
# digit integer which is optional if the default value is used:
#
# si: i'th stack word, default i=0 means top
# ni: i'th buffer word, default i=0 means first
# hi: i'th degree head, default i=1 means direct head
# li: i'th leftmost child, default i=1 means the leftmost child
# ri: i'th rightmost child, default i=1 means the rightmost child
# v: word vector
# c: context vector
# p: postag
# d: distance to the right.  e.g. s1d is s0s1 distance, s0d is s0n0
#    distance.  encoding: 1,2,3,4,5-10,10+ (from ZN11)
# L: dependency label (0 is ROOT or NONE)
# a: number of left children.  encoding: 0,1,...,8,9+
# b: number of right children.  encoding: 0,1,...,8,9+
# A: set of left dependency labels
# B: set of right dependency labels

# fs11h54h27 (9144): local best for acl11, ArcHybrid13, start=hybrid27, all=hybrid54, bparser10, shuffled
hybrid25 = split("n0A n0a n0lL n0lc n0lv n0c n0v n1c n1v s0A s0B s0a s0b s0d s0c s0rL s0v s1A s1B s1a s1rL s1rc s1rv s1v s2v")

# fs11e64e27 (9098): local best for acl11, ArcEager13, start=eager27=hybrid27, all=eager64, bparser10, shuffled
eager25 = split("n0A n0a n0l2L n0lL n0lc n0lv n0c n0v n1c n1v n2c n2v s0B s0c s0v s1A s1B s1a s1d s1c s1rL s1rc s1v s2c s2v")

# fs11h54h28 (8997): local best in acl11, ArcEager13 (by mistake), start=hybrid28, all=hybrid54
hybrid27 = split("n0A n0a n0lL n0lc n0lv n0c n0v n1c n1v s0A s0B s0a s0d s0c s0v s1A s1B s1a s1b s1d s1c s1rL s1rc s1rv s1v s2c s2v")
eager27 = hybrid27

# fs11h54h14 (8992): local best for acl11, ArcHybrid13, start=hybrid14, all=hybrid54
hybrid25b = split("n0a n0l2c n0l2v n0lc n0lv n0v n1c n1v n2v s0B s0a s0d s0l2L s0lL s0lc s0c s0rc s0v s1B s1b s1lL s1lc s1c s1r2c s1v")

# fs11e64e20 (8926): local best for acl11, ArcEager13, start=eager20, all=eager64
eager20 = split("n0A n0a n0lL n0c n0v n1c n1v n2c n2v s0B s0L s0b s0d s0hc s0hv s0c s0rL s0v s1c s1v")

# fs11e39e21 (8906): local best for acl11, ArcEager13, start=eager21, all=eager39
eager23 = split("n0A n0lL n0lc n0c n0v n1c n1v n2c n2v s0B s0L s0a s0b s0h2c s0hv s0l2L s0l2c s0l2v s0lc s0c s0r2L s0r2c s0v")

# fs07h54h28 (8797): local best for conll07, ArcHybrid13, start=hybrid28, all=hybrid54
hybrid28b = split("n0A n0a n0l2v n0lL n0lc n0c n0v n1c n1v n2v s0A s0B s0a s0b s0d s0c s0v s1A s1a s1b s1d s1c s1rL s1rc s1rv s1v s2c s2v")

# fs07h54h13 (8791): Local best for conll07, ArcHybrid13, start=tacl13hybrid (aka hybrid13), all=hybrid54
hybrid14 = split("n0lc n0lv n0v n1c n1v s0d s0lc s0c s0rc s0v s1B s1lc s1c s1v")

# fs07e39e39 (8725): local best in conll07, ArcEager13, start=acl11eager, all=acl11eager
eager36 = split("n0A n0a n0l2L n0l2c n0l2v n0lL n0lc n0lv n0c n0v n1c n1v n2c n2v s0A s0B s0L s0b s0d s0h2c s0hL s0hc s0hv s0l2L s0l2c s0l2v s0lL s0lc s0lv s0c s0r2c s0r2v s0rL s0rc s0rv s0v")

# fs07e39e21 (8705): Local best for conll07, ArcEager13, start=eager21, all=eager39
eager20b = split("n0l2L n0lv n0c n0v n1v n2c n2v s0B s0L s0a s0h2c s0hv s0l2L s0l2c s0l2v s0lc s0c s0r2L s0r2c s0v")

# fs07e64e20 (8677): Local best for conll07, ArcEager13, start=eager20, all=eager64
eager21b = split("n0A n0a n0l2L n0lL n0c n0v n1c n1v n2c n2v s0B s0L s0b s0d s0hc s0hv s0c s0rL s0v s1c s1v")

# more complete set than acl11eager: added s1 features
eager64 = split("n0A n0a n0l2L n0l2c n0l2v n0lL n0lc n0lv n0c n0v n1c n1v n2c n2v s0A s0B s0L s0a s0b s0d s0h2c s0h2v s0hL s0hc s0hv s0l2L s0l2c s0l2v s0lL s0lc s0lv s0c s0r2L s0r2c s0r2v s0rL s0rc s0rv s0v s1A s1B s1L s1a s1b s1d s1h2c s1h2v s1hL s1hc s1hv s1l2L s1l2c s1l2v s1lL s1lc s1lv s1c s1r2L s1r2c s1r2v s1rL s1rc s1rv s1v")

# Manually constructed starting point
eager21 = split("n0A n0lL n0c n0v n1c n1v n2c n2v s0B s0L s0a s0h2c s0hv s0l2L s0l2c s0l2v s0lc s0c s0r2L s0r2c s0v")

# Manually constructed starting point
hybrid28 = split("n0A n0a n0lL n0lc n0lv n0c n0v n1c n1v s0A s0B s0a s0b s0d s0c s0v s1A s1B s1a s1b s1d s1c s1rL s1rc s1rv s1v s2c s2v")

# Features that may be relevant for ArcHybrid13
hybrid54 = DFeature[
"n0v", "n0c", "n0A", "n0a",
"s0v", "s0c", "s0A", "s0a", "s0B", "s0b", "s0d",
"s1v", "s1c", "s1A", "s1a", "s1B", "s1b", "s1d",
"n1v", "n1c", "n2v", "n2c", "s2v", "s2c",
"n0lv", "n0l2v", "s0lv", "s0l2v", "s0rv", "s0r2v", "s1lv", "s1l2v", "s1rv", "s1r2v",
"n0lc", "n0l2c", "s0lc", "s0l2c", "s0rc", "s0r2c", "s1lc", "s1l2c", "s1rc", "s1r2c",
"n0lL", "n0l2L", "s0lL", "s0l2L", "s0rL", "s0r2L", "s1lL", "s1l2L", "s1rL", "s1r2L"
]

# Goldberg&Nivre TACL13 ArcHybrid features from: tacl2013dynamicoracles/lefttoright/features:HybridFeatures
# s0,s1,s2: stack word forms: a=-1,-2,-3 b=0 c=4
# w0,w1: buffer word forms: a=0,1 b=0 c=4
# Ts0: tag of s0 (-1 0 -4)
# Trcs1: tag of rightmost child of s1 (-2 1 -4)
# use word half of vector (c=+4) for form, context half of vector (c=-4) for tag
# No need for feature combinations, we have a nonlinear model.

tacl13hybrid = DFeature[
	# (1)
"s0v",	# append("s0_%s" % s0)
"s0c",	# append("Ts0_%s" % Ts0)
        # append("Ts0s0_%s_%s" % (Ts0, s0))

"s1v",	# append("s1_%s" % s1)
"s1c",	# append("Ts1_%s" % Ts1)
     	# append("Ts1s1_%s_%s" % (Ts1, s1))

"n0v",	# append("w0_%s" % w0)
"n0c",	# append("Tw0_%s" % Tw0)
	# append("Tw0w0_%s_%s" % (Tw0, w0))
	# +hybrid
"n1v",	# append("w1_%s" % w1)
"n1c",	# append("Tw1_%s" % Tw1)
	# append("Tw1w1_%s_%s" % (Tw1, w1))

	# (2)
	# append("s0s1_%s_%s" % (s0,s1))
	# append("Ts0Ts1_%s_%s" % (Ts0,Ts1))
	# append("Ts0Tw0_%s_%s" % (Ts0,Tw0))
	# append("s0Ts0Ts1_%s_%s_%s" % (s0,Ts0,Ts1))
	# append("Ts0s1Ts1_%s_%s_%s" % (Ts0,s1,Ts1))
	# append("s0s1Ts1_%s_%s_%s" % (s0,s1,Ts1))
	# append("s0Ts0s1_%s_%s_%s" % (s0,Ts0,s1))
	# append("s0Ts0Ts1_%s_%s_%s" % (s0,Ts0,Ts1))
	# +hybrid
	# append("s0w0_%s_%s" % (s0,w0))
	# append("Ts0Tw0_%s_%s" % (Ts0,Tw0))
	# append("Ts0Tw1_%s_%s" % (Ts0,Tw1))
	# append("s0Ts0Tw0_%s_%s_%s" % (s0,Ts0,Tw0))
	# append("Ts0w0Tw0_%s_%s_%s" % (Ts0,w0,Tw0))
	# append("s0w0Tw0_%s_%s_%s" % (s0,w0,Tw0))
	# append("s0Ts0w0_%s_%s_%s" % (s0,Ts0,w0))
	# append("w0Tw0Tw1_%s_%s_%s" % (w0,Tw0,Tw1)) #

	# (3)
	# append("Ts0Tw0Tw1_%s_%s_%s" % (Ts0,Tw0,Tw1))
	# append("Ts1Ts0Tw0_%s_%s_%s" % (Ts1,Ts0,Tw0))
	# append("s0Tw0Tw1_%s_%s_%s" % (s0,Tw0,Tw1))
	# append("Ts1s0Tw0_%s_%s_%s" % (Ts1,s0,Tw0))

	# (4) rc -1  lc 1
"s1rc",	# append("Ts1Trcs1Tw0_%s_%s_%s" % (Ts1, Trcs1, Tw0)) => Trcs1
	# append("Ts1Trcs1Tw0_%s_%s_%s" % (Ts1, Trcs1, Ts0))
"s1lc",	# append("Ts1Tlcs1Ts0_%s_%s_%s" % (Ts1, Tlcs1, Ts0)) => Tlcs1
"s0lc",	# append("Ts1Ts0Tlcs0_%s_%s_%s" % (Ts1, Ts0, Tlcs0)) => Tlcs0
"s0rc",	# append("Ts1Trcs0Ts0_%s_%s_%s" % (Ts1, Trcs0, Ts0)) => Trcs0
	# append("Ts0Tlcs1s0_%s_%s_%s" % (Ts0, Tlcs1, s0))
	# append("Ts1s0Trcs0_%s_%s_%s" % (Ts1, s0, Trcs0))
	# append("Ts0Tlcs1s0_%s_%s_%s" % (Ts0, Trcs1, s0))
	# append("Ts1s0Trcs0_%s_%s_%s" % (Ts1, s0, Tlcs0))

	# +hybrid
	# append("Ts0Trcs0Tw0_%s_%s_%s" % (Ts0,Trcs0,Tw0))
	# append("Ts0Trcs0Tw0_%s_%s_%s" % (Ts0,Trcs0,Tw1))
	# append("Ts0Tlcs0Tw0_%s_%s_%s" % (Ts0,Tlcs0,Tw0))
"n0lc",	# append("Ts0Tw0Tlcw0_%s_%s_%s" % (Ts0,Tw0,Tlcw0)) => Tlcw0
	# append("Ts0Tlcs0w0_%s_%s_%s" % (Ts0,Tlcs0,w0))
	# append("Ts0Tlcs0w0_%s_%s_%s" % (Ts0,Trcs0,w0))
	# append("Ts0w0Tlcw0_%s_%s_%s" % (Ts0,w0,Tlcw0))

	# append("Ts0Tw0Trcs0Tlcw0_%s_%s_%s_%s" % (Ts0,Tw0,Tlcw0,Trcs0))
	# append("Ts0Ts1Trcs1Tlcs0_%s_%s_%s_%s" % (Ts0,Ts1,Tlcs0,Trcs1))

	# (5)
	# append("Ts2Ts1Ts0_%s_%s_%s" % (Ts2,Ts1,Ts0))
	# append("Ts1Ts0Tw0_%s_%s_%s" % (Ts1,Ts0,Tw0))
	# append("Ts0Tw0Tw1_%s_%s_%s" % (Ts0,Tw0,Tw1))
] # tacl13hybrid

# Variants of tacl13hybrid:

tacl13words = DFeature["s0lv", "s0v", "s0rv", "s1lv", "s1v", "s1rv", "n0lv", "n0v", "n1v"]
tacl13tags  = DFeature["s0lc", "s0c", "s0rc", "s1lc", "s1c", "s1rc", "n0lc", "n0c", "n1c"]
tacl13wordtag = [tacl13words; tacl13tags]

# Renamed:
hybrid13 = tacl13hybrid

# Goldberg&Nivre TACL13 ArcEager features from:
# tacl2013dynamicoracles/lefttoright/features:EagerZhangNivre2011Extractor
# (seems to be the ones used in acl11)
# s0: top of stack (a=-1)
# n0,n1,n2: buffer (a=0,1,2)
# s0h,s0h2: parent and grandparent of s0
# s0l,s0r: left/rightmost child of s0
# s0l2,s0r2: second left/rightmost child of s0
# s0w,s0p,s0wp: form, postag, form-postag for s0
# d: refers to s0-n0 distance, encoded 1,2,3,...,9,10+
# s0vl,s0vr: number of left/right children of s0, flat encoding
# s0L: dependency label for s0 (NA for n0!)
# s0sl,s0sr: set of dependency labels to the left/right of s0

acl11eager = DFeature[
	# # Single Words
	# f("s0wp_%s" % (s0wp))
"s0v",	# f("s0w_%s"  % (s0w))
"s0c",	# f("s0p_%s"  % (s0p))
	# f("n0wp_%s" % (n0wp))
"n0v",	# f("n0w_%s"  % (n0w))
"n0c",	# f("n0p_%s"  % (n0p))
	# f("n1wp_%s" % (n1wp))
"n1v",	# f("n1w_%s"  % (n1w))
"n1c",	# f("n1p_%s"  % (n1p))
	# f("n2wp_%s" % (n2wp))
"n2v",	# f("n2w_%s"  % (n2w))
"n2c",	# f("n2p_%s"  % (n2p))
	#
	# # Pairs
	# f("s0wp,n0wp_%s_%s" % (s0wp, n0wp))
	# f("s0wp,n0w_%s_%s" % (s0wp, n0w))
	# f("s0w,n0wp_%s_%s" % (s0w, n0wp))
	# f("s0wp,n0p_%s_%s" % (s0wp, n0p))
	# f("s0p,n0wp_%s_%s" % (s0p, n0wp))
	# f("s0w,n0w_%s_%s" % (s0w, n0w)) #?
	# f("s0p,n0p_%s_%s" % (s0p, n0p))
	# f("n0p,n1p_%s_%s" % (n0p, n1p))
	#
	# # Tuples
	# f("n0p,n1p,n2p_%s_%s_%s" % (n0p, n1p, n2p))
	# f("s0p,n0p,n1p_%s_%s_%s" % (s0p, n0p, n1p))
	# f("s0hp,s0p,n0p_%s_%s_%s" % (s0hp, s0p, n0p)) => s0hp (handled in unigram section)
	# f("s0p,s0lp,n0p_%s_%s_%s" % (s0p, s0lp, n0p)) => s0lp (handled in unigram section)
	# f("s0p,s0rp,n0p_%s_%s_%s" % (s0p, s0rp, n0p)) => s0rp (handled in unigram section)
	# f("s0p,n0p,n0lp_%s_%s_%s" % (s0p, n0p, n0lp)) => n0lp (handled in unigram section)
	#
	# # Distance
"s0d",	# f("s0wd_%s:%s" % (s0w, d)) => d
	# f("s0pd_%s:%s" % (s0p, d))
	# f("n0wd_%s:%s" % (n0w, d))
	# f("n0pd_%s:%s" % (n0p, d))
	# f("s0w,n0w,d_%s:%s:%s" % (s0w, n0w, d))
	# f("s0p,n0p,d_%s:%s:%s" % (s0p, n0p, d))
	#
	# # Valence
"s0b",	# f("s0wvr_%s:%s" % (s0w, s0vr)) => s0vr
	# f("s0pvr_%s:%s" % (s0p, s0vr))
"s0a",	# f("s0wvl_%s:%s" % (s0w, s0vl)) => s0vl
	# f("s0pvl_%s:%s" % (s0p, s0vl))
"n0a",	# f("n0wvl_%s:%s" % (n0w, n0vl)) => n0vl
	# f("n0pvl_%s:%s" % (n0p, n0vl))
	#
	# # Unigrams
"s0hv",	# f("s0hw_%s" % (s0hw))
"s0hc",	# f("s0hp_%s" % (s0hp))
"s0L",	# f("s0L_%s" % (s0L))
	#
"s0lv",	# f("s0lw_%s" % (s0lw))
"s0lc",	# f("s0lp_%s" % (s0lp))
"s0lL",	# f("s0lL_%s" % (s0lL))
	#
"s0rv",	# f("s0rw_%s" % (s0rw))
"s0rc",	# f("s0rp_%s" % (s0rp))
"s0rL",	# f("s0rL_%s" % (s0rL))
	#
"n0lv",	# f("n0lw_%s" % (n0lw))
"n0lc",	# f("n0lp_%s" % (n0lp))
"n0lL",	# f("n0lL_%s" % (n0lL))
	#
	# # Third-order
	# #do we really need the non-grandparent ones?
"s0h2v",	# f("s0h2w_%s" % (s0h2w))
"s0h2c",	# f("s0h2p_%s" % (s0h2p))
"s0hL",	# f("s0hL_%s"  % (s0hL))
"s0l2v",	# f("s0l2w_%s" % (s0l2w))
"s0l2c",	# f("s0l2p_%s" % (s0l2p))
"s0l2L",	# f("s0l2L_%s" % (s0l2L))
"s0r2v",	# f("s0r2w_%s" % (s0r2w))
"s0r2c",	# f("s0r2p_%s" % (s0r2p))
"s0r2L",	# f("s0r2L_%s" % (s0r2L))
"n0l2v",	# f("n0l2w_%s" % (n0l2w))
"n0l2c",	# f("n0l2p_%s" % (n0l2p))
"n0l2L",	# f("n0l2L_%s" % (n0l2L))
	# f("s0p,s0lp,s0l2p_%s_%s_%s" % (s0p, s0lp, s0l2p))
	# f("s0p,s0rp,s0r2p_%s_%s_%s" % (s0p, s0rp, s0r2p))
	# f("s0p,s0hp,s0h2p_%s_%s_%s" % (s0p, s0hp, s0h2p))
	# f("n0p,n0lp,n0l2p_%s_%s_%s" % (n0p, n0lp, n0l2p))
	#
	# # Labels
"s0B",	# f("s0wsr_%s_%s" % (s0w, s0sr)) => s0sr
	# f("s0psr_%s_%s" % (s0p, s0sr))
"s0A",	# f("s0wsl_%s_%s" % (s0w, s0sl)) => s0sl
	# f("s0psl_%s_%s" % (s0p, s0sl))
"n0A",	# f("n0wsl_%s_%s" % (n0w, n0sl)) => n0sl
	# f("n0psl_%s_%s" % (n0p, n0sl))
] # tacl13eager

# Renamed:
eager39 = acl11eager
# fs11e39e39 (8906): eager39 also local best for acl11, ArcEager13, start=eager39, all=eager39

# Old feature format before we switched to ints:
# f::Features is a nx3 matrix whose rows determine which features to extract
# Each row of f consists of the following three values:
# 1. anchor word: 0:n0, 1:n1, 2:n2, ..., -1:s0, -2:s1, -3:s2, ...
# 2. target word: 0:self, 1:rightmost child, 2:second rightmost, -1:leftmost, -2:second leftmost ...
# 3. feature: One of the features listed below.
#
# 0:wvec (word+context if token, dim depends on encoding)
# +-1: exists/doesn't (one bit)
# +-2: right/left child count (4 bits, == encoding: 0,1,2,3+)
# +3: distance to right (4 bits, == encoding: 1,2,3,4+, root is 4+)
# -3: average of in-between tokens to the right (dim:wvec)
# +-4: word/context half of vector (dim:wvec/2, only valid for token encoding, assumes first half=word, second half=context)
# +-5: right/left child count (4 bits, >= encoding: >=1, >=2, >=3, >=4)
# +-6: right/left child count (4 bits, <= encoding: <=0, <=1, <=2, <=3)
# +7: distance to right, >= encoding (8 bits, >= encoding: >=2, >=3, >=4, >=6, >=8, >=12, >=16, >=20)
# -7: distance to right, <= encoding (8 bits, <= encoding: <=1, <=2, <=3, <=4, <=6, <=8, <=12, <=16)
# +-8: average of in-between word/context vectors to the right (dim:wvec/2)
# +-9: head exists/doesn't (one bit)

# 0.0418668 archybrid_conll07EnglishToken_wikipedia2MUNK-100_rbf376678_1014_cache.mat (trn/tst) (rbf-gamma:0.376747095368119)
fv022b = Int8[-3 0 8;-2 -1 -4;-2 -1 5;-2 0 -9;-2 0 -4;-2 0 4;-2 0 7;-2 1 1;-1 -1 -1;-1 -1 4;-1 0 -9;-1 0 -4;-1 0 4;-1 0 6;-1 1 4;0 -1 -4;0 -1 4;0 0 -4;0 0 -2;0 0 4;1 0 -4;1 0 4]

# 0.0325506937033084 archybrid_conllWSJToken_wikipedia2MUNK-100_rbf333140_cache.mat (rbf-gamma:0.333140404105682)
fv021a = Int8[-3 0 4;-3 1 -2;-2 0 4;-2 0 6;-2 1 -2;-1 -1 -4;-1 -1 2;-1 -1 4;-1 0 -4;-1 0 -1;-1 0 2;-1 0 4;-1 1 2;-1 1 4;0 -1 -4;0 -1 -1;0 -1 4;0 0 -4;0 0 4;1 0 -4;1 0 4]

# 0.0443168273919168 archybrid_conll07EnglishToken_wikipedia2MUNK-100_rbf372064_cache.mat (rbf-gamma:0.372063662109375)
fv019 = Int8[-3 0 8;-2 -1 -4;-2 0 -9;-2 0 -4;-2 0 4;-2 1 1;-1 -1 -1;-1 -1 4;-1 0 -4;-1 0 4;-1 0 6;-1 1 4;0 -1 -4;0 -1 4;0 0 -4;0 0 -2;0 0 4;1 0 -4;1 0 4]

# 0.03492252681764 archybrid_conllWSJToken_wikipedia2MUNK-100_rbf339506_cache.mat (rbf-gamma:0.339506144311523)
fv022a = Int8[-3 0 4;-3 1 -2;-2 0 -4;-2 0 4;-2 0 6;-2 1 -2;-1 -1 -4;-1 -1 2;-1 -1 4;-1 0 -4;-1 0 -1;-1 0 2;-1 0 4;-1 1 2;-1 1 4;0 -1 -4;0 -1 -1;0 -1 4;0 0 -4;0 0 4;1 0 -4;1 0 4]

# 0.0444685231336006 archybrid13_conll07EnglishToken_wikipedia2MUNK-100_rbf3.720637e+05_cache.mat
fv023a = Int8[-3 0 8;-2 -1 -4;-2 -1 -1;-2 0 -9;-2 0 -4;-2 0 4;-2 1 1;-2 1 4;-1 -1 -1;-1 -1 4;-1 0 -4;-1 0 4;-1 0 5;-1 0 6;-1 1 -4;-1 1 4;0 -1 -4;0 -1 4;0 0 -4;0 0 -2;0 0 4;1 0 -4;1 0 4]

# BUGGY: 0.0325637 Best with rbf=0.339506144311523, archybrid_conllWSJToken_wikipedia2MUNK-100_fv130_dump.mat.
fv023 = Int8[-3 0 4;-3 1 -2;-2 -1 -4;-2 0 -4;-2 0 4;-2 0 6;-2 1 -2;-1 -1 -4;-1 -1 2;-1 -1 4;-1 0 -4;-1 0 -1;-1 0 2;-1 0 4;-1 1 2;-1 1 4;0 -1 -4;0 -1 -1;0 -1 4;0 0 -4;0 0 4;1 0 -4;1 0 4]

# BUGGY: 0.0443493 Best with run_featselect_rbf('archybrid13', 'conll07EnglishToken_wikipedia2MUNK-100', 'fv136', 'fv021', 'rbf3721');
fv022 = Int8[-3 0 8;-2 -1 -1;-2 0 -9;-2 0 -4;-2 0 4;-2 1 1;-2 1 4;-1 -1 -1;-1 -1 4;-1 0 -4;-1 0 4;-1 0 6;-1 1 -4;-1 1 4;0 -1 -4;0 -1 4;0 -1 5;0 0 -4;0 0 -2;0 0 4;1 0 -4;1 0 4]

# 0.0445990 Best with rbf372111 archybrid_conll07EnglishToken_wikipedia2MUNK-100_fv136_dump.mat
fv021 = Int8[-3 0 8;-2 -1 -1;-2 0 -9;-2 0 -4;-2 0 4;-2 1 1;-1 -1 -1;-1 -1 4;-1 0 -4;-1 0 4;-1 0 5;-1 0 6;-1 1 -4;-1 1 4;0 -1 -4;0 -1 4;0 0 -4;0 0 -2;0 0 4;1 0 -4;1 0 4]

# Best on run_featselect22('archybrid', 'conllWSJToken_wikipedia2MUNK-100', 'fv130', 'fv018');
fv017a = Int8[-3 0 -4;-3 0 4;-2 0 -4;-2 0 4;-2 1 -2;-1 -1 -4;-1 -1 4;-1 0 -4;-1 0 4;-1 1 -2;-1 1 4;0 -1 -4;0 -1 4;0 0 -4;0 0 4;1 0 -4;1 0 4]

# This one does best on archybrid_conll07EnglishToken_wikipedia2MUNK-100_d5 (degree=5 kernel)
fv018a = Int8[-3 0 8;-2 0 -9;-2 0 -4;-2 0 4;-2 1 1;-1 -1 1;-1 -1 4;-1 0 -4;-1 0 4;-1 1 -4;-1 1 4;0 -1 -4;0 -1 4;0 0 -4;0 0 -2;0 0 4;1 0 -4;1 0 4]

# This one does best on archybrid_conll07EnglishToken_wikipedia2MUNK-100 (degree=3 kernel)
fv015a = Int8[-3 0 8;-2 0 -4;-2 0 4;-2 1 1;-1 -1 4;-1 0 -4;-1 0 4;-1 1 -4;-1 1 4;0 -1 -4;0 -1 4;0 0 -4;0 0 4;1 0 -4;1 0 4]

# This one is a local minimum on archybrid_conll07EnglishToken_wikipedia2MUNK-100 (degree=3 kernel)
fv031a = Int8[-3 0 -8;-3 0 -4;-3 0 -1;-3 0 8;-3 1 -4;-3 1 4;-2 -1 -4;-2 -1 4;-2 0 -8;-2 0 -4;-2 0 4;-2 0 8;-2 1 1;-2 1 2;-1 -1 -4;-1 -1 1;-1 -1 4;-1 0 -8;-1 0 -4;-1 0 2;-1 0 4;-1 0 8;-1 1 -4;-1 1 4;0 -1 -4;0 -1 1;0 -1 4;0 0 -4;0 0 4;1 0 -4;1 0 4]

fv034 = zeros(Int8, 0, 3)
for a=-3:1
  for b=-1:1
    if ((a >= 0) && (b > 0)) continue; end   # no rdeps for buffer words
    if ((a > 0)  && (b < 0)) continue; end   # no ldeps for buffer words other than n0
    if ((a <= -2) && (b < 0)) continue; end  # not interested in ldeps beyond s0
    if ((a <= -3) && (b > 0)) continue; end  # not interested in rdeps beyond s1
    for c=-9:9
      if ((c == 0) || (c == -3)) continue; end  # do not use the word+context combination, keep them separate
      if (abs(c) == 5 || abs(c) == 6) continue; end  # just use = encoding for child count
      if ((abs(c) == 2) && (b != 0)) continue; end # not interested in grand children
      if (abs(c) == 1) continue; end # existence not useful
      if (abs(c) == 9) continue; end # head not useful in archybrid
      if ((c == 3) || (c == -7)) continue; end  # just use >= encoding for distance
      if ((a > 0) && !in(c, [1,-1,4,-4])) continue; end      # no deps/dist/in-between/head for a>0
      if ((a == 0) && (b == 0) && !in(c, [1,-1,4,-4,-2,-5,-6])) continue; end  # no rdeps/dist/in-between/head for a=0
      if ((b != 0) && in(c, [3,-3,7,-7,8,-8,9,-9])) continue; end  # no dist/in-between/head for deps
      fv034 = [fv034; Int8[a b c]]
    end
  end
end

fv039 = Int8[
0 0 4; 		# [-196]    'n0w'
-1 0 4; 	# [-119]    's0w'
1 0 4; 		# [ -73]    'n1w'
-2 0 4; 	# [ -43]    's1w'
0 -1 4; 	# [ -20]    'n0l1w'
0 0 -4; 	# [ -20]    'n0c'
1 0 -4; 	# [ -19]    'n1c'
-1 0 -4; 	# [ -16]    's0c'
-1 -1 4; 	# [ -11]    's0l1w'
-2 0 -4; 	# [ -10]    's1c'
-3 -1 -4;     	# [  -8]    's2l1c'
-1 -1 -4;     	# [  -8]    's0l1c'
0 -1 -4; 	# [  -8]    'n0l1c'
-3 -1 2; 	# [  -7]    's2l1r='
-3 0 8; 	# [  -7]    's2aw'
-2 1 1; 	# [  -7]    's1r1+'
-1 0 -8; 	# [  -7]    's0ac'
-1 0 7; 	# [  -7]    's0d>'
-1 -1 1; 	# [  -6]    's0l1+'
-1 0 8; 	# [  -6]    's0aw'
-3 0 -4; 	# [  -5]    's2c'
-3 0 4; 	# [  -5]    's2w'
-2 -1 -4;     	# [  -5]    's1l1c'
-1 1 4; 	# [  -5]    's0r1w'
-3 0 -8; 	# [  -4]    's2ac'
0 -1 1; 	# [  -4]    'n0l1+'
-3 0 -1; 	# [  -3]    's2-'
-3 0 7; 	# [  -3]    's2d>'
-3 1 4; 	# [  -3]    's2r1w'
-2 -1 4; 	# [  -3]    's1l1w'
-2 0 8; 	# [  -3]    's1aw'
-2 1 2; 	# [  -3]    's1r1r='
-2 -1 -2;     	# [  -2]    's1l1l='
-1 -1 -2;     	# [  -2]    's0l1l='
-2 0 -8; 	# [  -1]    's1ac'
-1 -1 2; 	# [  -1]    's0l1r='
-1 1 -4; 	# [  -1]    's0r1c'
-3 1 -4; 	# [   0]    's2r1c'
-1 0 2; 	# [   0]    's0r='
]

fv012 = Int8[
0 0 4; 		# [-203]    n0w
-1 0 4; 	# [-125]    s0w
-2 0 4; 	# [ -59]    s1w
1 0 4; 		# [ -59]    n1w
0 -1 4; 	# [ -41]    n0l1w
-1 -1 4;	# [ -21]    s0l1w

0 0 -4; 	# [ -17]    n0c
-1 0 -4;	# [ -16]    s0c
-2 0 -4;	# [ -13]    s1c
1 0 -4; 	# [ -15]    n1c
0 -1 -4; 	# [ -11]    n0l1c
-1 -1 -4;	# [ -12]    s0l1c
]

fv017 = Int8[       # 0.0476108, 'archybrid', 'conll07EnglishToken_wikipedia2MUNK-100'
0 0 4; 		# [-203]    n0w
-1 0 4; 	# [-125]    s0w
-2 0 4; 	# [ -59]    s1w
1 0 4; 		# [ -59]    n1w
0 -1 4; 	# [ -41]    n0l1w
-2 1 -2;	# [ -30]    s1r1l=
-1 -1 4;	# [ -21]    s0l1w
-1 1 4; 	# [ -21]    s0r1w
-1 1 -2;	# [ -19]    s0r1l=
0 0 -4; 	# [ -17]    n0c
-1 0 -4;	# [ -16]    s0c
1 0 -4; 	# [ -15]    n1c
-2 0 -4;	# [ -13]    s1c
-3 0 4; 	# [ -12]    s2w
-1 -1 -4;	# [ -12]    s0l1c
-1 1 -4; 	# [ -12]    s0r1c
0 -1 -4; 	# [ -11]    n0l1c
]

fv018 = Int8[       # another version of fv808 that splits token features to word + context
    0 0 4;
    0 0 -4;
    -1 0 4;
    -1 0 -4;
    -2 0 4;
    -2 0 -4;
    1 0 4;
    1 0 -4;
    0 -1 4;
    0 -1 -4;
    -1 1 4;
    -1 1 -4;
    -2 1 -2;
    -1 -1 4;
    -1 -1 -4;
    -1 1 -2;
    -3 0 4;
    -3 0 -4;
]

fv084 = zeros(Int8, 0, 3)
for a=-3:1
  for b=-1:1
    if ((a >= 0) && (b > 0)) continue; end   # no rdeps for buffer words
    if ((a > 0)  && (b < 0)) continue; end   # no ldeps for buffer words other than n0
    for c=-9:9
      if ((c == 0) || (c == -3)) continue; end  # do not use the word+context combination, keep them separate
      if (abs(c) == 5 || abs(c) == 6) continue; end  # just use = encoding for child count
      if ((c == 3) || (c == -7)) continue; end  # just use >= encoding for distance
      if ((a > 0) && !in(c, [1,-1,4,-4])) continue; end      # no deps/dist/in-between/head for a>0
      if ((a == 0) && (b == 0) && !in(c, [1,-1,4,-4,-2,-5,-6])) continue; end  # no rdeps/dist/in-between/head for a=0
      if ((b != 0) && in(c, [3,-3,7,-7,8,-8,9,-9])) continue; end  # no dist/in-between/head for deps
      fv084 = [fv084; Int8[a b c]]
    end
  end
end

# Better initial starting point
fv008w = Int8[
    0  0 4; # n0w
    -1 0 4; # s0w
    1  0 4; # n1w
    -2 0 4; # s1w
    0 -1 4; # n0l1w
    -1 1 4; # s0r1w
    -2 1 4; # s1r1w
    -1 -1 4; # s0l1w
         ]

# Good initial starting point
fv008 = Int8[
    0 0 4;      # n0
    -1 0 4;     # s0
    1 0 4;      # n1
    -2 0 4;     # s1
    0 -1 4;     # n0l1
    -1 1 4;     # s0r1
    0 0 -4;     # n0c
    -1 0 -4;	# s0c
]

# Include the head feature for stack words:
fv136 = zeros(Int8, 0, 3)
for a=-3:2
  for b=-1:1
    if ((a >= 0) && (b > 0)) continue; end   # no rdeps for buffer words
    if ((a > 0)  && (b < 0)) continue; end   # no ldeps for buffer words other than n0
    for c=-9:9
      if ((c == 0) || (c == -3)) continue; end  # do not use the word+context combination, keep them separate
      if ((a > 0) && !in(c, [1,-1,4,-4])) continue; end      # no deps/dist/in-between/head for a>0
      if ((a == 0) && (b == 0) && !in(c, [1,-1,4,-4,-2,-5,-6])) continue; end  # no rdeps/dist/in-between/head for a=0
      if ((b != 0) && in(c, [3,-3,7,-7,8,-8,9,-9])) continue; end  # no dist/in-between/head for deps
      fv136 = [fv136; Int8[a b c]]
    end
  end
end

# All legal features for s2..n2:
fv130 = zeros(Int8, 0, 3)
for a=-3:2
  for b=-1:1
    if ((a >= 0) && (b > 0)) continue; end   # no rdeps for buffer words
    if ((a > 0)  && (b < 0)) continue; end   # no ldeps for a>0
    for c=-8:8
      if ((a > 0) && in(c, [2,-2,3,-3,5,-5,6,-6,7,-7,8,-8])) continue; end      # no deps/dist/in-between for a>0
      if ((a == 0) && (b == 0) && in(c, [2,3,-3,5,6,7,-7,8,-8])) continue; end  # no rdeps/dist/in-between for a=0
      if ((b != 0) && in(c, [3,-3,7,-7,8,-8])) continue; end  # no dist/in-between for deps
      if ((c == 0) || (c == -3)) continue; end  # do not use the word+context combination, keep them separate
      fv130 = [fv130; Int8[a b c]]
    end
  end
end

fv102 = zeros(Int8, 0, 3)
for a=-2:1
  for b=-1:1
    if ((a >= 0) && (b > 0)) continue; end   # no rdeps for buffer words
    if ((a > 0)  && (b < 0)) continue; end   # no ldeps for a>0
    for c=-8:8
      if ((a > 0) && in(c, [2,-2,3,-3,5,-5,6,-6,7,-7,8,-8])) continue; end      # no deps/dist/in-between for a>0
      if ((a == 0) && (b == 0) && in(c, [2,3,-3,5,6,7,-7,8,-8])) continue; end  # no rdeps/dist/in-between for a=0
      if ((b != 0) && in(c, [3,-3,7,-7,8,-8])) continue; end  # no dist/in-between for deps
      fv102 = [fv102; Int8[a b c]]
    end
  end
end

# fv1768 = Int8[
# #n0            n1             n2             s0             s1             s2             s0l1           s0l2           s0r1           s0r2           s1l1           s1l2           s1r1           s1r2           n0l1           n0l2          n0s0 s0s1
# 0  0  0  0  0  1  1  1  1  1  2  2  2  2  2 -1 -1 -1 -1 -1 -2 -2 -2 -2 -2 -3 -3 -3 -3 -3 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -2 -2 -2 -2 -2 -2 -2 -2 -2 -2 -2 -2 -2 -2 -2 -2 -2 -2 -2 -2  0  0  0  0  0  0  0  0  0  0 -1   -2;
# 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0 -1 -1 -1 -1 -1 -2 -2 -2 -2 -2  1  1  1  1  1  2  2  2  2  2 -1 -1 -1 -1 -1 -2 -2 -2 -2 -2  1  1  1  1  1  2  2  2  2  2 -1 -1 -1 -1 -1 -2 -2 -2 -2 -2  0    0;
# 0 -2  2  1 -1  0 -2  2  1 -1  0 -2  2  1 -1  0 -2  2  1 -1  0 -2  2  1 -1  0 -2  2  1 -1  0 -2  2  1 -1  0 -2  2  1 -1  0 -2  2  1 -1  0 -2  2  1 -1  0 -2  2  1 -1  0 -2  2  1 -1  0 -2  2  1 -1  0 -2  2  1 -1  0 -2  2  1 -1  0 -2  2  1 -1  3    3;
# ]';

fv1768 =
Int8[0 0 0
     0 0 -2
#    0 0 2   # buffer words dont have right children
     0 0 1
     0 0 -1
     1 0 0
#    1 0 -2  # buffer words other than n0 do not have ldeps
#    1 0 2   # buffer words dont have right children
     1 0 1
     1 0 -1
     2 0 0
#    2 0 -2  # buffer words other than n0 do not have ldeps
#    2 0 2   # buffer words dont have right children
     2 0 1
     2 0 -1
     -1 0 0
     -1 0 -2
     -1 0 2
     -1 0 1
     -1 0 -1
     -2 0 0
     -2 0 -2
     -2 0 2
     -2 0 1
     -2 0 -1
     -3 0 0
     -3 0 -2
     -3 0 2
     -3 0 1
     -3 0 -1
     -1 -1 0
     -1 -1 -2
     -1 -1 2
     -1 -1 1
     -1 -1 -1
     -1 -2 0
     -1 -2 -2
     -1 -2 2
     -1 -2 1
     -1 -2 -1
     -1 1 0
     -1 1 -2
     -1 1 2
     -1 1 1
     -1 1 -1
     -1 2 0
     -1 2 -2
     -1 2 2
     -1 2 1
     -1 2 -1
     -2 -1 0
     -2 -1 -2
     -2 -1 2
     -2 -1 1
     -2 -1 -1
     -2 -2 0
     -2 -2 -2
     -2 -2 2
     -2 -2 1
     -2 -2 -1
     -2 1 0
     -2 1 -2
     -2 1 2
     -2 1 1
     -2 1 -1
     -2 2 0
     -2 2 -2
     -2 2 2
     -2 2 1
     -2 2 -1
     0 -1 0
     0 -1 -2
     0 -1 2
     0 -1 1
     0 -1 -1
     0 -2 0
     0 -2 -2
     0 -2 2
     0 -2 1
     0 -2 -1
     -1 0 3
     -2 0 3]

fv804 = Int8[
#n0 s0 s1 n1 n0l1 s0r1 s0l1 s1r1 s0r
  0 -1 -2  1  0   -1   -1   -2   -1;
  0  0  0  0 -1    1   -1    1    0;
  0  0  0  0  0    0    0    0    2;
]';

fv708 = Int8[
#n0 n0l1 n0l2 n1 s0 s0r s0r1 s0s1 s1
  0  0    0    1 -1 -1  -1   -2   -2;
  0 -1   -2    0  0  0   1    0    0;
  0  0    0    0  0  2   0    3    0;
]';

fv808 = Int8[
#n0 s0 s1 n1 n0l1 s0r1 s1r1l s0l1 s0r1l s2
  0 -1 -2  1  0   -1   -2    -1   -1    -3;
  0  0  0  0 -1    1    1    -1    1     0;
  0  0  0  0  0    0   -2     0   -2     0;
]';

end # module Flist


text:
# Calculate velocities for tensors

"""
    phase_vels(C, az, inc) -> vp, vs1, vs2, pol, avs

Calculate the phase velocities for the 6x6 elasticity matrix `C` along the direction
(`az`, `inc`), in degrees, and return P-wave velocity `vp`, the fast and slow shear
wave velocities, `vs1` and `vs2`, the polarisation of the fast shear wave `pol`,
and the shear wave velocity anisotropy, `avs`.  Velocities are in m/s if the tensor
`C` is in m^2/s^2 (i.e., is a density-normalised tensor, sometimes called A).

`az` is the azimuth in degrees measured from the x1 towards the -x2 axis.

`inc` is the inclination in degrees from the x1-x2 plane towards the x3 axis.

`pol` is measured when looking towards the (0,0,0) point along the ray (in
the negative propagation direction).  `pol` increases clockwise away from the vector
normal to the propagation direction which points towards to the x3-axis.
"""
function phase_vels(C, az, inc)
    x = incaz2cart(inc, az)
    # Create the Christoffel matrix
    T = make_T(C, x)
    # Find eigenvectors and values, which correspond to velocities^2 (val[i])
    # and polarisation vectors (vec[:,i]).
    # They seem to be sorted into descending order of eigenvalue, but we double
    # check anyway.
    vals, vecs = LinearAlgebra.eigen(T)
    ip = argmax(vals)
    is2 = argmin(vals)
    is1 = 6 - ip - is2
    vp  = sqrt(vals[ip])
    vs1 = sqrt(vals[is1])
    vs2 = sqrt(vals[is2])
    # xp  = vec(vecs[:,ip])
    xs1 = vecs[:,is1]
    # xs2 = vec(vecs[:,is2])
    # Calculate S1 polarisation and amount of shear wave anisotropy
    pol = get_pol(inc, az, x, xs1)
    avs = 200*(vs1 - vs2)/(vs1 + vs2)
    return (vp=vp, vs1=vs1, vs2=vs2, pol=pol, avs=avs)
end

function group_vels(C, az, inc)
    # Return the group velocities in the az, inc direction for a 6x6
    # Voigt elasticity matrix C, in an elastic medium.
    # INPUT:
    #    C(6,6)  : Elasticity matrix
    #    az      : Azimuth (degrees) away from x1 towards -x2 axis
    #    inc     : Inclination (degrees) away from x1-x2 plane towards x3

    x = incaz2cart(inc, az)
    # Create the Christoffel matrix
    T = make_T(C, x)
    # Find eigenvectors, giving phase velocities
    vals, vecs = LinearAlgebra.eigen(T)
    ip = indmax(vals)
    is2 = indmin(vals)
    is1 = 6 - ip - is2
    vp = sqrt(vals[ip])
    vs1 = sqrt(vals[is1])
    vs2 = sqrt(vals[is2])
    xp = vec(vecs[:,ip])
    xs1 = vec(vecs[:,is1])
    xs2 = vec(vecs[:,is2])
    # Calculate group velocities, which are phase velocities projected onto
    # group velocity directions
    pp = xp./vp
    ps1 = xs1./vs1
    ps2 = xs2./vs2
    p = [pp'; ps1'; ps2']
    vg = zero(MVector{3,float(eltype(C))})
    ijkl = VOIGT_CONTRACTION_MATRIX
    for i=1:3, j=1:3, k=1:3, l=1:3
        m = ijkl[i,j]
        n = ijkl[k,l]
        vg[i] = vg[i] + C[m,n]*p[i,k]*x[j]*x[l]
    end
    return vg[1], vg[2], vg[3]
end

"Projection of fast shear wave polarisation onto wavefront plane"
function get_pol(inc, az, x, xs1)
    xs1p = cross(x, cross(x, xs1))
    xs1p /= norm(xs1p)
    # Local up vector
    u = incaz2up(inc, az)
    # Angle between the local up vector and the projected fast orientation
    v = cross(xs1p, u)
    pol = rad2deg(atan(norm(v), dot(xs1p, u)))
    # If v is codirectional with x, then the angle is correct; otherwise, we're
    # the wrong way round
    if dot(x, v) < 0
        pol = -pol
    end
    # Put in range -90° to 90°
    return mod(pol + 90.0, 180.0) - 90.0
end

"""
    incaz2up(inc, az) -> ::SVector{3}

For a vector pointing along a given azimuth `az` and inclination `inc`,
return the local 'up' direction, which is normal to the original direction and
within the plane defined by it and the z-axis.
"""
function incaz2up(inc, az)
    a = deg2rad(az) + pi
    i = deg2rad(inc)
    sini = sin(i)
    return @SVector [cos(a)*sini, -sin(a)*sini, cos(i)]
end

"""
    make_T(C, x) -> T

Create the Christoffel matrix, `T`, given the Voigt matrix `C` and unit vector `x`.
"""
function make_T(C, x)
    T = zero(MMatrix{3,3,float(eltype(C))})
    ijkl = VOIGT_CONTRACTION_MATRIX
    @inbounds for i = 1:3, j = 1:3, k = 1:3, l = 1:3
        m = ijkl[i,j]
        n = ijkl[k,l]
        T[i,k] = T[i,k] + C[m,n]*x[j]*x[l]
    end
    return LinearAlgebra.Hermitian(T) # Real-symmetric
end

"Return a 3-vector which is the cartesian direction corresponding to
 the azimuth `az` and inclination `inc` (degrees)"
function incaz2cart(inc, az)
    i = deg2rad(inc)
    a = deg2rad(az)
    cosi = cos(i)
    return @SVector [cos(a)*cosi, -sin(a)*cosi, sin(i)]
end


text:
# Residual block from Putzky and Welling (2019): https://arxiv.org/abs/1911.10914
# Author: Philipp Witte, pwitte3@gatech.edu
# Date: January 2020

export FluxBlock

"""
    FB = FluxBlock(model::Chain)

 Create a (non-invertible) neural network block from a Flux network.

 *Input*:

 - `model`: Flux neural network of type `Chain`

 *Output*:

 - `FB`: residual block layer

 *Usage:*

 - Forward mode: `Y = FB.forward(X)`

 - Backward mode: `ΔX = FB.backward(ΔY, X)`

 *Trainable parameters:*

 - Network parameters given by `Flux.parameters(model)`

 See also:  [`Chain`](@ref), [`get_params`](@ref), [`clear_grad!`](@ref)
"""
mutable struct FluxBlock <: NeuralNetLayer
    model::Chain
    params::Array{Parameter, 1}
end

@Flux.functor FluxBlock

#######################################################################################################################
# Constructor

function FluxBlock(model::Chain)

    # Collect Flux parameters
    model_params = Flux.params(model)
    nparam = length(model_params.order)
    params = Array{Parameter}(undef, nparam)

    # Create InvertibleNetworks parameter
    for j=1:nparam
        params[j] = Parameter(model_params.order[j])
    end
    return FluxBlock(model, params)
end

#######################################################################################################################
# Functions

# Forward
forward(X::AbstractArray{Float32, 4}, FB::FluxBlock) = FB.model(X)

# Backward 2D
function backward(ΔY::AbstractArray{Float32, 4}, X::AbstractArray{Float32, 4}, FB::FluxBlock; set_grad::Bool=true)

    # Backprop using Zygote
    θ = Flux.params(X, FB.model)
    back = Zygote.pullback(() -> FB.model(X), θ)[2]
    grad = back(ΔY)

    # Set gradients
    ΔX = grad[θ[1]]
    if set_grad
        for j=1:length(FB.params)
            FB.params[j].grad = grad[θ[j+1]]
        end
    else
        Δθ = Array{Parameter, 1}(undef, length(FB.params))
        for j=1:length(FB.params)
            Δθ[j] = grad[θ[j+1]]
        end
    end

    set_grad ? (return ΔX) : (ΔX, Δθ)

end

## Jacobian utilities

function jacobian(ΔX::AbstractArray{Float32, 4}, Δθ::Array{Parameter, 1}, X::AbstractArray{Float32, 4}, FB::FluxBlock)
    throw(ArgumentError("Jacobian for Flux block not yet implemented"))
end

function adjointJacobian(ΔY::AbstractArray{Float32, 4}, X::AbstractArray{Float32, 4}, FB::FluxBlock)
    return backward(ΔY, X, FB; set_grad=false)
end

## Other utils

# Clear gradients
function clear_grad!(FB::FluxBlock)
    nparams = length(FB.params)
    for j=1:nparams
        FB.params[j].grad = nothing
    end
end

"""
    P = get_params(NL::NeuralNetLayer)

 Returns a cell array of all parameters in the network layer. Each cell
 entry contains a reference to the original parameter; i.e. modifying
 the paramters in `P`, modifies the parameters in `NL`.
"""
function get_params(FB::FluxBlock)
    return FB.params
end

function set_params!(FB::FluxBlock, θ::Array{Parameter, 1})
    model_params = Flux.params(FB.model)
    nparams = length(model_params.order)
    for j=1:nparams
        model_params.order[j] .= θ[j].data
        FB.params[j].data = θ[j].data
        FB.params[j].grad = θ[j].grad
    end
end


text:
# Top-level functions that dispatch to the individual constraints
function TO.cost!(J::Vector{<:Real}, conSet::ALConstraintSet)
    for i in eachindex(conSet.convals)
        TO.cost!(J, conSet.convals[i])
    end
end

function TO.cost_expansion!(E, conSet::ALConstraintSet)
    for i in eachindex(conSet.convals)
        TO.cost_expansion!(E, conSet.convals[i])
    end
end

# Cost function on a single constraint
function TO.cost!(J, conval::ALConVal)
    cone = TO.sense(conval)
    for (i,k) in enumerate(conval.inds)
        c = SVector(conval.vals[i])
        λ = SVector(conval.λ[i])
        μ = SVector(conval.μ[i])
        J[k] += TO.cost(cone, λ, c, μ)
    end
end

TO.cost(::TO.Equality, λ, c, μ) = λ'c + 0.5*c'Diagonal(μ)*c

function TO.cost(::TO.Inequality, λ, c, μ)
    a = @. (c >= 0) | (λ > 0)
	return λ'c + 0.5*c'Diagonal(μ .* a)*c
end

function TO.cost(cone::TO.SecondOrderCone, λ, c, μ)
	λbar = TO.projection(cone, λ - μ .* c)
	Iμ = Diagonal(1 ./ μ)
	return 0.5 * (λbar'Iμ*λbar - λ'Iμ*λ)
end

# Cost expansion on a single constraint across all time steps
function TO.cost_expansion!(E, conval)
    cone = TO.sense(conval)
    for (i,k) in enumerate(conval.inds)
        TO.cost_expansion!(cone, conval, i)
    end
    # add the values to the existing cost
    copy_expansion!(E, conval)
end

function TO.cost_expansion!(::Equality, conval, i)
    c = SVector(conval.vals[i])
    ∇c = conval.jac[i]
    λ = SVector(conval.λ[i])
    μ = conval.μ[i][1]

    λbar = λ + μ * c
    conval.const_hess[i] = true
    if prod(size(∇c)) < 24*24
        ∇c = SMatrix(∇c)
        conval.grad[i] .= ∇c'λbar
        conval.hess[i] .= μ*∇c'∇c
    else
        mul!(conval.grad[i].data, Transpose(∇c.data), λbar)
        mul!(conval.hess[i].data, Transpose(∇c), ∇c)
        conval.hess[i] .*= μ
    end
    return
end

function TO.cost_expansion!(::Inequality, conval, i)
    c = SVector(conval.vals[i])
    ∇c = conval.jac[i]
    tmp = conval.tmp
    λ = SVector(conval.λ[i])
    # μ = SVector(conval.μ[i])
    μ = conval.μ[i][1]
    a = @. (c >= 0) | (λ > 0)
    Iμ = Diagonal(a)

    λbar = λ + μ * (a .* c)
    conval.const_hess[i] = false

    if prod(size(∇c)) < 24*24
        ∇c = SMatrix(∇c)
        conval.grad[i] .= ∇c'λbar
        conval.hess[i] .= μ*∇c'Iμ*∇c
    else
        mul!(conval.grad[i].data, Transpose(∇c.data), λbar)
        mul!(tmp.data, Iμ, ∇c.data)  # TODO: do this directly on the SizedArrays
        mul!(conval.hess[i].data, Transpose(∇c), tmp)
        conval.hess[i] .*= μ
    end
end

function TO.cost_expansion!(cone::SecondOrderCone, conval, i)
    c = SVector(conval.vals[i])
    ∇c = conval.jac[i]
    λ = SVector(conval.λ[i])
    μ = SVector(conval.μ[i])
    Iμ = Diagonal(μ)
    ∇proj = conval.∇proj[i]  # pxp projection jacobian
    ∇²proj = conval.∇²proj[i]
    tmp = conval.tmp

    # The term inside the projection operator
    λbar = λ - μ .* c
    conval.const_hess[i] = false

    # Evaluate the projection and it's derivatives
    λp = TO.projection(cone, λbar)      # TODO: don't project more than you need to!
    TO.∇projection!(cone, ∇proj, λbar)        # evaluate the Jacobian of the projection operation
    TO.∇²projection!(cone, ∇²proj, λbar, λp)  # evalute the Jacobian of ∇Π(λ - μ*c)'Π(λ - μ*c)

    # Apply the chain rule
    μ = μ[1]
    mul!(tmp, ∇proj, ∇c)
    mul!(conval.grad[i], Transpose(tmp), λp)
    conval.grad[i] .*= -1                      # -∇cproj'λp
    mul!(conval.hess[i], Transpose(tmp), tmp)  # ∇cproj'∇cproj
    mul!(tmp, ∇²proj, ∇c)
    mul!(conval.hess[i], Transpose(∇c), tmp, 1.0, 1.0)
    conval.hess[i] .*= μ
end

function copy_expansion!(E::TO.CostExpansion{n,m}, conval::ALConVal{<:TO.StateConstraint}) where {n,m}
    ix,iu = SVector{n}(1:n), SVector{m}(1:m) .+ n
    for (i,k) in enumerate(conval.inds)
        E[k].q .+= conval.grad[i]
        E[k].Q .+= conval.hess[i]
        # E.const_hess[k] &= conval.is_const[i] & conval.const_hess[i]
    end
end

function copy_expansion!(E::TO.CostExpansion{n,m}, conval::ALConVal{<:TO.ControlConstraint}) where {n,m}
    ix,iu = SVector{n}(1:n), SVector{m}(1:m) .+ n
    for (i,k) in enumerate(conval.inds)
        E[k].r .+= conval.grad[i]
        E[k].R .+= conval.hess[i]
        # E.const_hess[k] &= conval.is_const[i] & conval.const_hess[i]
    end
end

function copy_expansion!(E::TO.CostExpansion{n,m}, conval::ALConVal{<:TO.StageConstraint}) where {n,m}
    ix,iu = SVector{n}(1:n), SVector{m}((1:m) .+ n)
    for (i,k) in enumerate(conval.inds)
        E[k].q .+= conval.grad[i].data[ix]
        E[k].r .+= conval.grad[i].data[iu]
        E[k].Q .+= conval.hess[i].data[ix,ix]
        E[k].R .+= conval.hess[i].data[iu,iu]
        E[k].H .+= conval.hess[i].data[iu,ix]
        # E.const_hess[k] &= conval.is_const[i] & conval.const_hess[i]
    end
end

function copy_expansion!(E::TO.CostExpansion{n,m}, conval::ALConVal{<:TO.CoupledConstraint}) where {n,m}
    ix,iu = SVector{n}(1:n), SVector{m}((1:m) .+ n)
    for (i,k) in enumerate(conval.inds)
        E[k].q .+= conval.grad[i][ix]
        E[k].r .+= conval.grad[i][iu]
        E[k].Q .+= conval.hess[i][ix,ix]
        E[k].R .+= conval.hess[i][iu,iu]
        E[k].H .+= conval.hess[i][iu,ix]
        E[k+1].q .+= conval.grad[i,2][ix]
        E[k+1].r .+= conval.grad[i,2][iu]
        E[k+1].Q .+= conval.hess[i,2][ix,ix]
        E[k+1].R .+= conval.hess[i,2][iu,iu]
        E[k+1].H .+= conval.hess[i,2][iu,ix]
        # E.const_hess[k] &= conval.is_const[i] & conval.const_hess[i]
    end
end


text:
# State tests

# Discrete state constructors
@test size(DiscreteState(1.0).x) == (1,)
@test size(DiscreteState(1.0, 0).x) == (1,)
@test size(UncertainDiscreteState(1.0, 2.0).x) == (1,)
@test size(UncertainDiscreteState(1.0, 2.0, 0).P) == (1, 1)

# Continuous state constructors
@test size(ContinuousState(1.0).x) == (1,)
@test size(ContinuousState(1.0, 0.0).x) == (1,)
@test size(UncertainContinuousState(1.0, 2.0).x) == (1,)
@test size(UncertainContinuousState(1.0, 2.0, 0.0).P) == (1, 1)

# Create testing states
discrete_state = DiscreteState(ones(2))
discrete_state = DiscreteState(ones(2), 0)
unc_discrete_state = UncertainDiscreteState(ones(2), eye(2))
unc_discrete_state = UncertainDiscreteState(ones(2), eye(2), 0)
continuous_state = ContinuousState(ones(2))
continuous_state = ContinuousState(ones(2), 0.0)
unc_continuous_state = UncertainContinuousState(ones(2), eye(2))
unc_continuous_state = UncertainContinuousState(ones(2), eye(2), 0.0)

# Test type conversions
@test StateEstimation.absolute_type(unc_discrete_state) <: DiscreteState
@test StateEstimation.absolute_type(unc_continuous_state) <: ContinuousState
@test StateEstimation.uncertain_type(discrete_state) <: UncertainDiscreteState
@test StateEstimation.uncertain_type(continuous_state) <: UncertainContinuousState
@test discrete_state == make_absolute(unc_discrete_state)
@test unc_discrete_state == make_uncertain(discrete_state, unc_discrete_state.P)
@test make_uncertain(discrete_state, zeros(2,2)) ==
    make_uncertain(discrete_state)
@test continuous_state == make_absolute(unc_continuous_state)
@test unc_continuous_state == make_uncertain(continuous_state,
                                             unc_continuous_state.P)
@test make_uncertain(continuous_state, zeros(2,2)) ==
    make_uncertain(continuous_state)

# Sample uncertain states
@test typeof(sample(unc_discrete_state)) <: DiscreteState
@test typeof(sample(unc_continuous_state)) <: ContinuousState

# Test distance metrics
@test distance(discrete_state, discrete_state) == 0
@test distance(discrete_state, unc_discrete_state) == 0
@test mahalanobis(discrete_state, unc_continuous_state) == 0
@test mahalanobis(DiscreteState([4.0, 1.0]), unc_discrete_state) == 3
@test mahalanobis(discrete_state, unc_discrete_state) ==
    mahalanobis(discrete_state, discrete_state, unc_discrete_state.P)

# Test absolute state expansion/reduction
test_state = ContinuousState(1.0)
new_state = expand_state(test_state, [2.0])
@test length(new_state.x) == 2
new_state = reduce_state(new_state, 1)
@test length(new_state.x) == 1

# Test uncertain state expansion/reduction
test_state = UncertainDiscreteState(ones(3), eye(3))
@test_throws DimensionMismatch expand_state(test_state, ones(2), eye(3))
new_state = expand_state(test_state, 2, 2)
@test length(new_state.x) == 4
@test size(new_state.P) == (4,4)
new_state = reduce_state(new_state, 3)
@test test_state.x == new_state.x
@test test_state.P == new_state.P

function test_state_addition_subtraction(state::AbstractState)
    for a = -1:0.5:1
        @test (state + a).x == (state.x+a)
        @test (a + state).x == (state.x+a)
        @test (state - a).x == (state.x-a)
        @test (a - state).x == (state.x-a)
        new_state = deepcopy(state)
        new_state .-= a
        new_state .+= a
        new_state .= a .+ new_state
        new_state .= a .- new_state
        @test new_state.x == state.x
    end
end
test_state_addition_subtraction(DiscreteState(ones(2)))
test_state_addition_subtraction(ContinuousState(ones(2)))
test_state_addition_subtraction(UncertainDiscreteState(ones(2), eye(2)))
test_state_addition_subtraction(UncertainContinuousState(ones(2), eye(2)))

function test_state_scalar_multiplication(state::AbstractState)
    for a = -1:0.5:1
        @test (state*a).x == (state.x*a)
        @test (a*state).x == (state.x*a)
        new_state = deepcopy(state)
        new_state .*= a
        @test new_state.x == (state.x*a)
        new_state .= a .* new_state
        @test new_state.x == (state.x*a^2)
        if typeof(state) <: AbstractUncertainState
            @test (state*a).P == (state.P*a^2)
            @test (a*state).P == (state.P*a^2)
            new_state = deepcopy(state)
            new_state .*= a
            @test new_state.P == (state.P*a^2)
            new_state .= a .* new_state
            @test new_state.P == (state.P*a^4)
        end
    end
end
test_state_scalar_multiplication(DiscreteState(ones(2)))
test_state_scalar_multiplication(ContinuousState(ones(2)))
test_state_scalar_multiplication(UncertainDiscreteState(ones(2), eye(2)))
test_state_scalar_multiplication(UncertainContinuousState(ones(2), eye(2)))

function test_state_matrix_multiplication(state::AbstractState)
    test_matrices = [eye(2), ones(2,2), -ones(2,2), zeros(2,2)]
    for idx = 1:length(test_matrices)
        A = test_matrices[idx]
        @test (A*state).x == A*state.x
        new_state = deepcopy(state)
        new_state .= A * state
        @test new_state.x == A*state.x
        if typeof(state) <: AbstractUncertainState
            @test (A*state).P == A*state.P*A'
            new_state = deepcopy(state)
            new_state .= A * state
            @test new_state.P == A*state.P*A'
        end
    end
end
test_state_scalar_multiplication(DiscreteState(ones(2)))
test_state_scalar_multiplication(ContinuousState(ones(2)))
test_state_scalar_multiplication(UncertainDiscreteState(ones(2), eye(2)))
test_state_scalar_multiplication(UncertainContinuousState(ones(2), eye(2)))

# Test in-place assignment operations
state1 = DiscreteState(ones(2))
state1 .= DiscreteState(zeros(2))
@test state1.x == zeros(2)
state1 = ContinuousState(ones(2))
state1 .= ContinuousState(zeros(2))
@test state1.x == zeros(2)
state1 = UncertainDiscreteState(ones(2), eye(2))
state1 .= UncertainDiscreteState(zeros(2), zeros(2,2))
@test (state1.x == zeros(2)) & (state1.P == zeros(2,2))
state1 = UncertainContinuousState(ones(2), eye(2))
state1 .= UncertainContinuousState(zeros(2), zeros(2,2))
@test (state1.x == zeros(2)) & (state1.P == zeros(2,2))

# State equality checks
@test DiscreteState(ones(2)) == DiscreteState([1.0, 1.0])
@test DiscreteState(ones(2)) != DiscreteState(zeros(2))
@test UncertainDiscreteState(ones(2), eye(2)) ==
    UncertainDiscreteState([1.0, 1.0], eye(2))
@test UncertainDiscreteState(ones(2), eye(2)) != DiscreteState(ones(2))
@test (DiscreteState(ones(2)) == ContinuousState(ones(2))) == false
@test ContinuousState(ones(2)) == ContinuousState([1.0, 1.0])
@test ContinuousState(ones(2)) != ContinuousState(zeros(2))
@test UncertainContinuousState(ones(2), eye(2)) ==
    UncertainContinuousState([1.0, 1.0], eye(2))
@test UncertainContinuousState(ones(2), eye(2)) != ContinuousState(ones(2))