text:
using LsqFit
using DataArrays, DataFrames

import Base.getindex
import Base.setindex!
import Base.show

module auxiliary end

type CannotFitException <: Exception
    msg::AbstractString
end

type BadDataException <: Exception
    msg::AbstractString
end

type NoResultsException <: Exception
    msg::AbstractString
end

"""
Describes a set of data and a model to be fit to.
"""
type Fitter

    "Model function to fit to."
    f::Function

    "Data to be fitted to."
    data::DataFrame

    "Fit parameters with initial guesses / best-fit values and uncertainties."
    _parameters::ModelParameters

    "Whether the most recent fitting attempt has converged."
    _converged::Bool

    "Datapoints to be ignored when fitting."
    _outliers::BitArray{1}

    "Studentized residuals of the fit."
    _residuals::Nullable{Vector{AbstractFloat}}

    "Covariance matrix of the best-fit parameters."
    _covariance::Union{Array{Float64, 2}, Void}

    "Various settings for the fitter."
    _settings::Dict{Symbol, Any}

    "Function to use for fitting, taking only the required two parameters."
    _f_fitting::Function

    "Figure number of the associated plot."
    _figure_number::Integer

end

"""
    Fitter(f::Function; kwargs...)

Creates a new Fitter object with model function `f`.
"""
function Fitter(f::Function; kwargs...)

    residuals = Nullable{Array}()

    arg_names = argument_names(f)
    variable_name = arg_names[1]
    param_names = arg_names[2:end]

    parameters = ModelParameters(
        [(p, FreeParameter(i, 1.0)) for (i, p) in enumerate(param_names)])

    settings = Dict{Symbol, Any}(
            :error_range     => 0.68,
            :autoplot        => true,
            :xscale          => "linear",
            :yscale          => "linear",
            :plot_curve      => true,
            :plot_guess      => true,
            :fpoints         => 1000,
            :xmin            => nothing,
            :xmax            => nothing,
            :xlabel          => "$(variable_name)",
            :ylabel          => "f($(variable_name))",
            :xvar            => variable_name,
            :yvar            => :y,
            :eyvar           => :y_err,
            :style_data      => Dict(:marker => "+", :color => "b",   :ls => ""),
            :style_outliers  => Dict(:marker => "+", :color => "0.5", :ls => ""),
            :style_fit       => Dict(:marker => "",  :color => "r",   :ls => "-"),
            :style_guess     => Dict(:marker => "",  :color => "k",   :ls => "--"))

    merge!(settings, Dict(kwargs))

    n_parameters = number_of_arguments(f) - 1 # remove one for the x

    if n_parameters == 0
        msg = "Cannot fit to a function with only one free parameter."
        throw(CannotFitException(msg))
    end

    f_fitting(x, p) = f(x, p...)

    data          = DataFrame()
    converged     = false
    covariance    = nothing
    figure_number = -1
    outliers      = BitArray{1}()

    Fitter(
        f, data, parameters, converged, outliers, residuals, covariance,
        settings, f_fitting, figure_number)
end

function setindex!(fitter::Fitter, val, key)
    if key ∈ keys(fitter._settings)
        fitter._settings[key] = val
    elseif key ∈ keys(fitter._parameters)
        fitter._parameters[key].value = val
    else
        throw(KeyError(key))
    end
end

function getindex(fitter::Fitter, key)
    val = nothing
    if key ∈ keys(fitter._settings)
        val = fitter._settings[key]
    elseif key ∈ keys(fitter._parameters)
        if isa(fitter._parameters[key], FixedParameter)
            val = fitter._parameters[key].value
        else
            if fitter._converged
                val = get(fitter._parameters[key].fit_value)
            else
                msg = "Fit results cannot be accessed until " *
                      "`fit!` has been called successfully."
                throw(NoResultsException(msg))
            end
        end
    else
        throw(KeyError(key))
    end
    val
end

function show(stream::IO, fitter::Fitter)
    description = "EasyPhys.Fitter with the following settings:\n\n"
    for (k, v) in fitter._settings
        description *= "\t$(rpad(k, 15)) => $(repr(v))\n"
    end
    description *= "\n"

    if n_free_parameters(fitter) < length(fitter._parameters)
        description *= "Constants:\n\n"
        for (key, val) in fitter._parameters
            if isa(val, FixedParameter)
                description *= "\t$(rpad(key, 15)) = $(val.value)\n"
            end
        end
        description *= "\n"
    end

    try
        χ²_guess = reduced_χ²(fitter, guesses(fitter))
        description *= "Guesses (χ² = $(χ²_guess)):\n\n"
    catch e
        if isa(e, BadDataException)
            description *= "Guesses:\n\n"
        else
            rethrow(e)
        end
    end

    for (key, val) in sort(
            [(k, v) for (k, v) in fitter._parameters if isa(v, FreeParameter)],
            by=(kv -> kv[2].position))
        description *= "\t$(rpad(key, 15)) = $(fitter._parameters[key].value)\n"
    end
    description *= "\n"

    if fitter._converged
        χ²_fit = reduced_χ²(fitter)
        confidence_percent = round(Integer, fitter[:error_range] * 100)
        description *= "Best-fit parameters with $(confidence_percent)\% "
        description *= "confidence interval (χ² = $(χ²_fit)):\n\n"
        for (key, val) in sort(
                [(k, v) for (k, v) in fitter._parameters if isa(v, FreeParameter)],
                by=(kv -> kv[2].position))
            description *= "\t$(rpad(key, 15)) = "
            description *= "$(get(val.fit_value)) ± "
            description *= "$(get(val.fit_uncertainty))\n"
        end
    else
        description *= "Fit results not yet present."
    end
    description *= "\n"

    write(stream, description)
end

"""
    null_results!(fitter::Fitter)

    null_results!(model_parameters::ModelParameters)

Sets all the `best_fit_value` and `fit_uncertainty` fields of `fitter`
or `model_parameters` to null.
"""
function null_results!(model_parameters::ModelParameters)
    for val in values(model_parameters)
        if isa(val, FreeParameter)
            val.fit_value       = Nullable{AbstractFloat}()
            val.fit_uncertainty = Nullable{AbstractFloat}()
        end
    end
end

null_results!(fitter::Fitter) = null_results!(fitter._parameters)

"""
    free!(fitter::Fitter, parameters...; parameter_value_pairs...)

    (fitter::Fitter) |> free!(parameters...; parameter_value_pairs...)

Frees `parameters` from any constraints when fitting. Guesses can
be provided using keyword arguments; otherwise, previous values are
used. Returns `fitter` so that similar calls can be chained together.
"""
@partially_applicable function free!(
        fitter::Fitter, parameters...; parameter_value_pairs...)

    for parameter in parameters
        p = fitter._parameters[parameter]
        fitter._parameters[parameter] = FreeParameter(p.position, p.value)
    end

    for (parameter, guess) in parameter_value_pairs
        p = fitter._parameters[parameter]
        fitter._parameters[parameter] = FreeParameter(p.position, guess)
    end

    null_results!(fitter)

    fitter
end

"""
    fix!(fitter::Fitter; parameter_value_pairs...)

    (fitter::Fitter) |> fix!(; parameter_value_pairs...)

For each `{parameter => value}` pair in `parameter_value_pairs`, fixes
`parameter` to `value` when fitting, treating it as a constant. Returns
`fitter` so that similar calls can be chained together.
"""
@partially_applicable function fix!(fitter::Fitter; parameter_value_pairs...)
    for (parameter, value) in parameter_value_pairs
        p = fitter._parameters[parameter]
        fitter._parameters[parameter] = FixedParameter(p.position, value)
    end

    null_results!(fitter)

    fitter
end

"""
    guess!(fitter::Fitter, values::Dict{Symbol, AbstractFloat})

    (fitter::Fitter) |> guess!(values::Dict{Symbol, AbstractFloat})

    guess!(fitter::Fitter, values)

    (fitter::Fitter) |> guess!(values)

    guess!(fitter::Fitter; parameter_value_pairs...)

    (fitter::Fitter) |> guess!(; parameter_value_pairs...)

Sets the initial guesses used to fit the free parameters of `fitter`.
Returns `fitter` so that similar calls can be chained together.
"""
@partially_applicable function guess!(
        fitter::Fitter, values::Dict{Symbol, AbstractFloat})

    for (k, v) in values
        parameter = fitter._parameters[k]
        if isa(parameter, FreeParameter)
            parameter.value = v
        else
            msg = "$(k) is a fixed parameter. Its value can be set using " *
                  "`fix!`, or it can be freed using `free!`"
            warn(msg)
        end
    end

    fitter
end

@partially_applicable function guess!(fitter::Fitter, values)
    free_parameters =
        [k for (k, v) in fitter._parameters if isa(v, FreeParameter)]

    if length(values) ≠ length(free_parameters)
        msg = "If supplied as an Array, `values` must contain exactly one " *
              "entry for each free parameter associated with `fitter`."
        warn(msg)
    end

    guess!(fitter, Dict{Symbol, AbstractFloat}(zip(free_parameters, values)))
end

@partially_applicable function guess!(fitter::Fitter; parameter_value_pairs...)
    guess!(fitter, Dict{Symbol, AbstractFloat}(parameter_value_pairs))
end

"""
    guesses(fitter::Fitter)

    guesses(model_parameters::ModelParameters)

Returns the initial guesses (or, for fixed parameters, the values) of
all fit parameters associated with `fitter` or `model_parameters`.
"""
function guesses(model_parameters::ModelParameters)
    guesses = ones(length(model_parameters))
    for val in values(model_parameters)
        guesses[val.position] = val.value
    end

    guesses
end

guesses(fitter::Fitter) = guesses(fitter._parameters)

"""
    n_free_parameters(fitter::Fitter)

    n_free_parameters(model_parameters::ModelParameters)

Returns the number of free (non-constant) parameters associated with
`fitter` or `model_parameters`.
"""
function n_free_parameters(model_parameters::ModelParameters)
    length([v for v in values(model_parameters) if isa(v, FreeParameter)])
end

n_free_parameters(fitter::Fitter) = n_free_parameters(fitter._parameters)

"""
    xdata(fitter::Fitter)

Gets the dependent data associated with `fitter`.
"""
xdata(fitter::Fitter) = fitter.data[fitter[:xvar]]

"""
    ydata(fitter::Fitter)

Gets the independent data associated with `fitter`.
"""
ydata(fitter::Fitter) = fitter.data[fitter[:yvar]]

"""
    eydata(fitter::Fitter)

Gets the errors in the independent data associated with `fitter`.
"""
eydata(fitter::Fitter) = fitter.data[fitter[:eyvar]]

"""
    set!(fitter::Fitter; kwargs...)

    (fitter::Fitter) |> set!(; kwargs...)

Updates the settings of `fitter` from the values given in `kwargs`.
Returns `fitter` so that similar calls can be chained together.
"""
@partially_applicable function set!(fitter::Fitter; kwargs...)
    for (k, v) in kwargs
        fitter[k] = v
    end

    fitter
end

"""
    set_data!(fitter::Fitter, xdata, ydata, eydata)

    set_data!(fitter::Fitter, dataframe)

    (fitter::Fitter) |> set_data!(xdata, ydata, eydata)

    (fitter::Fitter) |> set_data!(dataframe)

Updates the data associated with `fitter`. If a DataFrame is supplied, three
columns must be present in the order [xdata, ydata, eydata]. Returns `fitter`
so that similar calls can be chained together.
"""
@partially_applicable function set_data!(fitter::Fitter, xdata, ydata, eydata)
    n_data = length(xdata)
    if length(ydata) ≠ n_data
        msg = "xdata and ydata must have the same number of data"
        throw(BadDataException(msg))
    end

    if typeof(eydata) <: Number
        eydata = eydata * ones(xdata)
    elseif length(eydata) ≠ n_data
        msg = "eydata must be broadcastable to the size of xdata and ydata"
        throw(BadDataException(msg))
    end

    dataframe = DataFrame(Any[xdata, ydata, eydata], [:x, :y, :y_err])

    set_data!(fitter, dataframe)
end

@partially_applicable function set_data!(
        fitter::Fitter, dataframe::DataFrame,
        xvar=nothing, yvar=nothing, eyvar=nothing)
    if (size(dataframe, 2) ≠ 3)
        msg = "Exactly 3 columns of data must be supplied."
        throw(BadDataException(msg))
    end

    if  (xvar === nothing)  xvar = names(dataframe)[1] end
    if  (yvar === nothing)  yvar = names(dataframe)[2] end
    if (eyvar === nothing) eyvar = names(dataframe)[3] end

    fitter[:xvar] = xvar
    fitter[:yvar] = yvar
    fitter[:eyvar] = eyvar

    n_data = size(dataframe, 1)
    fitter._outliers = BitArray{1}(repeat([false], outer=n_data))

    fitter.data = dataframe

    if fitter[:autoplot]
        plot!(fitter)
    end

    fitter
end

"""
    apply_mask!(fitter::Fitter, mask::BitArray{1})

    (fitter::Fitter) |> apply_mask!(mask::BitArray{1})

Applies `mask` to the data associated with `fitter`. Returns `fitter`
so that similar calls can be chained together.
"""
@partially_applicable function apply_mask!(
        fitter::Fitter, mask::AbstractArray{Bool, 1})
    fitter._outliers = (@__dot__ ~mask)

    if fitter[:autoplot]
        plot!(fitter)
    end

    fitter
end

"""
    ignore_outliers!(fitter::Fitter, max_residual=10, params=[])

    (fitter::Fitter) |> ignore_outliers!(max_residual=10, params=[])

Ignores outlying datapoints associated with `fitter`, defined as those points
the studentized residuals of which are larger than `max_residual` when the
model function is evaluated with `params`. If no `params` are provided, uses
fit results.
"""
@partially_applicable function ignore_outliers!(
        fitter::Fitter, max_residual=10, params=[])
    outliers = abs.(studentized_residuals(fitter, params)) .> max_residual

    fitter |> apply_mask!(@__dot__ ~outliers)
end

"""
    xlims(fitter::Fitter)

Gets the active x-limits `xmin` and `xmax` of `fitter`.
"""
function xlims(fitter::Fitter)
    xmin = fitter[:xmin] === nothing ? minimum(xdata(fitter)) : fitter[:xmin]
    xmax = fitter[:xmax] === nothing ? maximum(xdata(fitter)) : fitter[:xmax]

    (xmin, xmax)
end

"""
    data_mask(fitter::Fitter)

Returns a mask indexing data which are within the fitting limits of `fitter`,
as determined by `fitter[:xmin]` and `fitter[:xmax]` as well as any outliers.
"""
function data_mask(fitter::Fitter)
    xmin, xmax = xlims(fitter)

    @compat (xmin .<= xdata(fitter) .<= xmax) .& (@__dot__ ~fitter._outliers)
end

"""
    fit!(fitter::Fitter; kwargs...)

    (fitter::Fitter) |> fit!(; kwargs...)

Fits the model function of `fitter` to the associated data. Updates the
settings of `fitter` from the values given in `kwargs` before fitting.
Returns `fitter` so that similar calls can be chained together.
"""
@partially_applicable function fit!(fitter::Fitter; kwargs...)
    if size(fitter.data, 1) == 0
        msg = "All xdata, ydata, and eydata must be set before calling `fit!`."
        throw(BadDataException(msg))
    end

    set!(fitter; kwargs...)

    fit_mask = data_mask(fitter)
    xdata_fit  = Vector{Float64}( xdata(fitter)[fit_mask])
    ydata_fit  = Vector{Float64}( ydata(fitter)[fit_mask])
    eydata_fit = Vector{Float64}(eydata(fitter)[fit_mask])

    all_params = sort(
        collect(keys(fitter._parameters)),
        by=(k -> fitter._parameters[k].position))

    fitting_params = filter(
        (k -> isa(fitter._parameters[k], FreeParameter)), all_params)

    if isempty(fitting_params)
        msg = "No free parameters! Cannot fit."
        throw(CannotFitException(msg))
    end

    fitting_constants = Dict{Symbol, AbstractFloat}(
        [(key, val.value) for (key, val) in fitter._parameters if isa(val, FixedParameter)])

    fitting_function_args = [:fitter, fitter[:xvar], fitting_params...]

    eval(
        auxiliary,
        quote
            fitting_function = $(Expr(:tuple, fitting_function_args...)) ->
                begin
                    $([:($k = $v) for (k, v) in fitting_constants]...)
                    fitter._f_fitting($(fitter[:xvar]), [$(all_params...)])
                end
        end)

    fitting_function = (x, p) ->
        eval(auxiliary, :(fitting_function($(fitter), $(x), $(p)...)))

    guesses = [fitter._parameters[k].value for k in fitting_params]
    weights = 1 ./ abs.(eydata_fit)

    fit_results = curve_fit(
        fitting_function, xdata_fit, ydata_fit, weights, guesses)

    if fit_results.converged
        fitter._converged = true

        params = fit_results.param
        param_errors = estimate_errors(fit_results, fitter[:error_range])

        for (i, (param, param_error)) in enumerate(zip(params, param_errors))
            fitter._parameters[fitting_params[i]].fit_value =
                Nullable{AbstractFloat}(param)
            fitter._parameters[fitting_params[i]].fit_uncertainty =
                Nullable{AbstractFloat}(param_error)
        end

        fitter._covariance = estimate_covar(fit_results)
        fitter._residuals = Nullable{Array}(-fit_results.resid)

        if fitter[:autoplot]
            plot!(fitter)
        end
    else
        fitter._converged = false

        null_results!(fitter)

        msg = "Fit did not converge."
        warn(msg)
    end

    fitter
end

"""
    parameter_covariance(fitter::Fitter)

Computes the covariance matrix of the best-fit parameters associated
with `fitter`.
"""
parameter_covariance(fitter::Fitter) = fitter._covariance

"""
    studentized_residuals(fitter::Fitter)

Computes the studentized residuals of the fit described by `fitter`.
"""
function studentized_residuals(fitter::Fitter, params=[])
    if size(fitter.data, 1) == 0
        msg = "All xdata, ydata, and eydata must be set " *
              "in order to calculate residuals."
        throw(BadDataException(msg))
    end

    resids = []
    if isempty(params)
        if fitter._converged
            resids = get(fitter._residuals)
        else
            msg = "Fit must converge before residuals can be accessed."
            throw(NoResultsException(msg))
        end
    else
        fit_mask = data_mask(fitter)
        weights = 1 ./ abs.(eydata(fitter)[fit_mask])

        resids = ydata(fitter)[fit_mask]
               - apply_f(fitter, xdata(fitter)[fit_mask], params)

        resids .*= weights
    end

    resids
end

"""
    reduced_χ²(fitter::Fitter)

Computes the reduced χ² of the fit described by `fitter`.
"""
function reduced_χ²(fitter::Fitter, params=[])
    squared_residuals = studentized_residuals(fitter, params).^2

    ndof = countnz(data_mask(fitter)) - n_free_parameters(fitter)

    sum(squared_residuals) / ndof
end

"""
    apply_f(fitter::Fitter, x)

    apply_f(fitter::Fitter, x, params)

    (fitter::Fitter) |> apply_f(x)

Applies the model function of `fitter` at points `x` using parameters `params`.
If no `params` is supplied and fit! has been called, defaults to using the
best-fit parameters.
"""
@partially_applicable function apply_f(fitter::Fitter, x, params=[])
    if isempty(params)
        if fitter._converged
            params = ones(length(fitter._parameters))
            for val in values(fitter._parameters)
                if isa(val, FixedParameter)
                    params[val.position] = val.value
                else
                    params[val.position] = get(val.fit_value)
                end
            end
        else
            msg = "Parameters must be explicitly supplied."
            throw(NoResultsException(msg))
        end
    end
    fitter._f_fitting(x, params)
end


text:
function plot_points()
    sdf = stacked_dataset()
    dropmissing!(sdf)
    xy = data(sdf) * mapping(:group, :score; row=:domain)
    draw(xy)
end

function plot_density()
    sdf = stacked_dataset()
    dropmissing!(sdf)

    dens = data(sdf) * mapping(:score, color=:group, row=:domain) * density()
    # Doesn't work for some reason.
    # axis = (; ylabel="")
    draw(dens)
end

ext: jl


text:
using ChiralWaveguide
using Test

include("WavePackets.jl")

include("WaveguideProblems.jl")

include("Systems.jl")

include("solve.jl")

include("ExactResults.jl")


text:
# This file is for g2o integration. For more information on the actual file
# formats, refer to: https://github.com/RainerKuemmerle/g2o/wiki/File-Format

export importG2o, exportG2o, findCommand, parseG2oInstruction!

## Common functions for g2o parsing

global commands = Dict(Pose2Pose2 => :EDGE_SE2,
                       Pose2Point2BearingRange => :LANDMARK)
#

"""
    $SIGNATURES

Hackish reverse lookup for dict[val] <= key.

Notes
- See DFG BiDictMap as better future solution (TODO)
"""
function findCommand(comm::Symbol)
  global commands
  for (key, val) in commands
    if val == comm
      return key
    end
  end
  error("$comm not found in known g2o commands, $(collect(values(commands)))")
  return nothing
end

"""
    $SIGNATURES

Read in every line from a g2o file.
"""
function importG2o(input_file::String)
    instructions = []
    # Read input file line by line.
    open(input_file) do file
        for ln in eachline(file)
            pieces = split(ln, ' ', keepempty=false)
            push!(instructions, pieces)
        end
    end
    return instructions
end

"""
    $SIGNATURES

Parses a single g2o instruction to add information to factor graph.
"""
function parseG2oInstruction!(fg::AbstractDFG,
                              instruction::Array{SubString{String},1})
    if instruction[1] == "EDGE_SE2"
        # Need to add a relative pose measurement between two variables.
        # Parse all of the variables starting with the symbols.
        from_pose = Symbol("x", instruction[2])
        to_pose = Symbol("x", instruction[3])

        # Parse the mean values for the relative pose measurements.
        x_state = parse(Float64, instruction[4])
        y_state = parse(Float64, instruction[5])
        theta_state = parse(Float64, instruction[6])

        # Parse the information matrix and take inverse to obtain covariance.
        info_mat = [parse(Float64, instruction[7]) parse(Float64, instruction[8])  parse(Float64, instruction[9]);
                    parse(Float64, instruction[8]) parse(Float64, instruction[10]) parse(Float64, instruction[11]);
                    parse(Float64, instruction[9]) parse(Float64, instruction[11]) parse(Float64, instruction[12])]
        cov_mat = inv(info_mat)
        # NOTE workaround to ensure cov_mat is Hermitian -- TODO use inf_mat directly for MvNormal
        cov_mat += cov_mat'
        cov_mat ./= 2

        # Make sure both variables are in the FG. Otherwise add them.
        if (from_pose in ls(fg)) == false
            addVariable!(fg, from_pose, Pose2)
        end
        if (to_pose in ls(fg)) == false
            addVariable!(fg, to_pose, Pose2)
        end

        # Add a factor between the two poses with the relative measurement.
        measurement = MvNormal([x_state; y_state; theta_state], cov_mat)
        rel_pose_factor = Pose2Pose2(measurement)
        addFactor!(fg, [from_pose, to_pose], rel_pose_factor)
    end
    return fg
end

## Export g2o functions

function getVariableListInts!(fct, uniqVarInt, varIntLabel)
  varlist = Int[]
  for vari in solverData(fct).fncargvID
    if !haskey(varIntLabel, vari)
      uniqVarInt[1] += 1
      varIntLabel[vari] = uniqVarInt[1]
    end
    push!(varlist, varIntLabel[vari])
  end
  varlist
end

function stringG2o!(dfg::AbstractDFG,
                    fc::Symbol,
                    fnc::Pose2Pose2,
                    varIntLabel::Dict,
                    uniqVarInt::Vector{Int};
                    overwriteMapping::Dict=Dict())::String
  #
  global commands
  # get variable numbers
  varlist = getVariableListInts!(getFactor(dfg,fc), uniqVarInt, varIntLabel)
  # get information matrix
  INF = 1 ./ sqrt(fnc.z.Σ.mat)
  INF[INF .== Inf] .= 0
  # get command
  comm = !haskey(overwriteMapping, Pose2Pose2) ? commands[Pose2Pose2] : overwriteMapping[Pose2Pose2]
  return "$comm $(varlist[1]) $(varlist[2]) $(fnc.z.μ[1]) $(fnc.z.μ[2]) $(fnc.z.μ[3]) $(INF[1,1]) $(INF[1,1]) $(INF[1,2]) $(INF[1,3]) $(INF[2,2]) $(INF[2,3]) $(INF[2,3]) $(INF[2,3])"
end

function stringG2o!(dfg::AbstractDFG,
                    fc::Symbol,
                    fnc::Pose2Point2BearingRange,
                    varIntLabel::Dict,
                    uniqVarInt::Vector{Int};
                    overwriteMapping::Dict=Dict{Symbol,Symbol}())::String
  #
  global commands
  # get variable numbers
  varlist = getVariableListInts!(getFactor(dfg,fc), uniqVarInt, varIntLabel)
  # get information matrix
  INF = [1/sqrt(fnc.bearing.σ); 0; 1/sqrt(fnc.range.σ)]
  # get command
  comm = !haskey(overwriteMapping, Pose2Point2BearingRange) ? commands[Pose2Point2BearingRange] : overwriteMapping[Pose2Point2BearingRange]
  return "$comm $(varlist[1]) $(varlist[2]) $(fnc.bearing.μ) $(fnc.range.μ) $(INF[1]) $(INF[2]) $(INF[3])"
end

function stringG2o!(dfg::AbstractDFG,
                    fc::Symbol,
                    fnc,
                    varIntLabel::Dict,
                    uniqVarInt::Vector{Int};
                    overwriteMapping::Dict=Dict{Symbol,Symbol}())
  #
  global commands
  error("unknown factor type $fnc")
end

"""
    $SIGNATURES

Export a factor graph to g2o file format.

Note:
- This funtion only supports Gaussian (i.e. Normal/MvNormal) factors, unpredictable witchcraft is used in other cases such as `AliasingScalarSampler` factor models.
"""
function exportG2o(dfg::AbstractDFG;
                   poseRegex::Regex=r"x\d",
                   solvable::Int=0,
                   ignorePriors::Bool=true,
                   filename::AbstractString="/tmp/test.txt",
                   overwriteMapping::Dict=Dict{Symbol, Symbol}())
  #
  uniqVarInt = Int[-1;]
  varIntLabel = Dict{Symbol, Int}()
  # all variables
  vars = ls(dfg, poseRegex, solvable=solvable) |> sortDFG
  # all factors
  fcts = lsf(dfg, solvable=solvable)
  # build text file based on factors, using pose variable order as guide
  io = open(filename,"w")
  for vs in vars
    # assign a unique number
    # all factors connected to variable
    vfcs = ls(dfg, vs, solvable=solvable)
    kvfcs = intersect(vfcs, fcts)
    for fc in kvfcs
      # ignore priors
      ignorePriors && isPrior(dfg, fc) ? (filter!(x->x!=fc, fcts); continue) : nothing
      # only add factors to g2o file once, remove if found
      !(fc in fcts) ? continue : filter!(x->x!=fc, fcts)
      # actually add the factor to the file
      fnc = getFactorType(dfg, fc)
      pstr = stringG2o!(dfg, fc, fnc, varIntLabel, uniqVarInt, overwriteMapping=overwriteMapping)
      println(io, pstr)
    end
  end
  close(io)

  return filename
end

##


text:
using Plots
using MieScattering
plotlyjs()
################################################################
# Fig 1.

# structure
m = 1.33 + .01im   # imag. component to demonstrate absorption

N = 1000
x = range(0., 13, length=N+1)[2:end]

Qsca = zeros(N)
Qabs = zeros(N)
for i = 1:N
    local qext, qsca, qabs, qback
    qext, qsca, qabs, qback = mie_scattering(m, x[i])
    Qsca[i] = qsca
    Qabs[i] = qabs
end

plt = plot(x, Qsca,
    xlabel="ka",
    ylabel="Scattering efficiency",
    legend=false,
    label="LABEL",
    title="TITLE",
    yscale=:log10,
    ylims=(0.005, maximum(Qsca))
)
plot!(x, Qabs)
display(plt)

println("max:" , maximum(Qsca), " ind: ",  argmax(Qsca))
println("min:" , minimum(Qsca), " ind: ",  argmin(Qsca))

################################################################
# Fig. 2

# structure
m = 2 + 1im   # imag. component to demonstrate absorption

N=1000
x = range(0, 13, length=N+1)[2:end]

Qsca = zeros(N)
Qabs = zeros(N)
for i = 1:N
    local qext, qsca, qabs, qback
    qext, qsca, qabs, qback = mie_scattering(m, x[i])
    Qsca[i] = qsca
    Qabs[i] = qabs
end

plt = plot(x, Qsca,
    xlabel="ka",
    ylabel="Scattering efficiency",
    legend=false,
    label="LABEL",
    title="TITLE",
    yscale=:log10,
    ylims=(0.005, maximum(Qabs))
)
plot!(x, Qabs)
display(plt)

println("max:" , maximum(Qsca), " ind: ",  argmax(Qsca))
println("min:" , minimum(Qsca), " ind: ",  argmin(Qsca))

################################################################
# Fig. 3.

m = 1.33 + .01im
x = 3

theta = range(0, π, length=2*180+1)
SL, SR, SU = scattering_function(theta, m, x)
SU /= π * x^2

angles = theta*180/π
# TODO: check why sqrt
plot(
    angles, SU,
    label="Unpolarized",
    ylabel="Intensity (\$|S|^2}\$)",
    xlabel="phase angle ϴ",
    yaxis=:log10
)

qext, qsca, qabs, qback = mie_scattering(m, x)
plot!(angles, qsca/4π * ones(length(angles)))

# plot!(theta*180/π, SL, label="left-polarized (perpendicular)")
# plot!(theta*180/π, SR, label="right-polarized (parallel)")

println("max:" , maximum(SU), " ind: ", angles[argmax(SU)])
println("min:" , minimum(SU), " ind: ", angles[argmin(SU)])

m = 1.33+0.01im
x = 7

theta = range(0, π, length=2*180+1)
SL, SR, SU = scattering_function(theta, m, x)
SU /= π * x^2

angles = theta*180/π
# TODO: check why sqrt
plot(
    angles, SU,
    label="Unpolarized",
    ylabel="Intensity (\$|S|^2}\$)",
    xlabel="phase angle ϴ",
    yaxis=:log10
)

qext, qsca, qabs, qback = mie_scattering(m, x)
plot!(angles, qsca/4π * ones(length(angles)))

# plot!(theta*180/π, SL, label="left-polarized (perpendicular)")
# plot!(theta*180/π, SR, label="right-polarized (parallel)")

println("max:" , maximum(SU), " ind: ", angles[argmax(SU)])
println("min:" , minimum(SU), " ind: ", angles[argmin(SU)])

m = 2+1im
x = 7

theta = range(0, π, length=2*180+1)
SL, SR, SU = scattering_function(theta, m, x)
SU /= π * x^2

angles = theta*180/π
# TODO: check why sqrt
plot(
    angles, SU,
    label="Unpolarized",
    ylabel="Intensity (\$|S|^2}\$)",
    xlabel="phase angle ϴ",
    yaxis=:log10
)

qext, qsca, qabs, qback = mie_scattering(m, x)
plot!(angles, qsca/4π * ones(length(angles)))

# plot!(theta*180/π, SL, label="left-polarized (perpendicular)")
# plot!(theta*180/π, SR, label="right-polarized (parallel)")

println("max:" , maximum(SU), " ind: ", angles[argmax(SU)])
println("min:" , minimum(SU), " ind: ", angles[argmin(SU)])

max_stars_repo_name: fekad/jlmie

max_issues_repo_name: fekad/jlmie

max_forks_repo_name: fekad/jlmie


text:
input = joinpath(@__DIR__, "input")
foods = readlines(input)

function parse_contents(foods)
    acc = []
    for food in foods
        (recipes, allergens) = split(food, " (contains ")
        recipes = map(s -> strip(s), split(recipes, " "))
        allergens = map(s -> strip(s), split(allergens[1:end - 1], ","))
        push!(acc, (recipes, allergens))
    end
    return acc
end

function tainted(foods)
    acc = Dict{String, Vector{String}}()
    for (recipe, allergens) in foods
        for allergen in allergens
            if haskey(acc, allergen)
                as = acc[allergen]
                acc[allergen] = intersect(as, recipe)
            else
                acc[allergen] = recipe
            end
        end
    end
    return acc
end

function all_ingredients(foods)
    acc = Set{String}()
    for (recipe, _) in foods
        for ingredient in recipe
            push!(acc, ingredient)
        end
    end
    return acc
end

function all_tainted_ingredients(tainted)
    acc = Set{String}()
    for (_, ingredients) in tainted
        for ingredient in ingredients
            push!(acc, ingredient)
        end
    end
    return acc
end

function count_okay(dishes, okay)
    count = 0
    for (ingredients, _) in dishes
        for ingredient in ingredients
            if ingredient in okay
                count += 1
            end
        end
    end
    return count
end

function solve(tainted)
    result = Dict{String, String}()
    seen = []

    while length(result) < length(tainted)
        for (k, v) in tainted
            left = setdiff(v, seen)
            if length(left) == 1
                append!(seen, left)
                result[k] = left[1]
            end
        end
    end

    return result
end

dishes = parse_contents(foods)
t = tainted(dishes)
a = all_ingredients(dishes)
at = all_tainted_ingredients(t)
okay = setdiff(a, at)

solved = solve(t)
solved = collect(solved)
sort!(solved, by=(kv -> first(kv)))

p1 = count_okay(dishes, okay)
p2 = join([v for (_, v) in solved], ",")

println("-----------------------------------------------------------------------")
println("allergen assessment -- part one :: $p1")
println("allergen assessment -- part two :: $p2")
println("-----------------------------------------------------------------------")

@assert(p1 == 2203)
@assert(p2 == "fqfm,kxjttzg,ldm,mnzbc,zjmdst,ndvrq,fkjmz,kjkrm")


text: Base.getproperty(ele::HTMLElement, name::String) = getattr(ele, name)

max_stars_repo_name: oxinabox/oxinabox

max_issues_repo_name: oxinabox/oxinabox

max_forks_repo_name: oxinabox/oxinabox


text:
#generate some preliminary graphs of beta

#-----> primary entry point for preliminary beta
function formpreliminarybeta()

  local beta::DataFrame
  local crsp::DataFrame

  #form the data to grph
  beta = readprocessbetadata()
  validpermnos = unique(beta.permno)
  crsp = prepcrsp(
    crsptype = :prelimcrspm,
    refreshcrsp=PARAM[:refreshpreliminarybeta],
    crsppath=PARAM[:preliminarypath],
    validpermnos = validpermnos,
    csvextension = CSV_EXTENSION,
    incsvstream = IN_CSV_STREAM,
    binprefix = "prelimbeta",
    crspcolumns = [:permno, :date, :price, :mktcap])::DataFrame

  beta = innerjoin(beta, crsp, on=[:date, :permno])

  series = aggregatebetadata(beta)
  PRELIM_STYLE_IDX::Dict = Dict(:textsize => 12.0,
      :linesize => 0.5,
      :shapesize => 0.25,
      #:xdomain => [1926, 2018],
      :width => 400,
      :height => 300,
      :xaxisname=>"Year",
      :legendtextsize => 9.0,
      :ydomain => [-1.0,2.0]#,
      #=:seriescolor => :black=#)

  #println(describe(series))
  graphbetadata(series,
    styleidx=PRELIM_STYLE_IDX,
    graphname="loadings",
    yaxisname = "Beta")

  graphbetadata(series,
    Fy=:valueabs,
    styleidx=PRELIM_STYLE_IDX,
    graphname="absloadings",
    yaxisname = "abs(Beta)")

end

function graphbetadata(
  series::AbstractDataFrame;
    styleidx::Dict = error("styleidx is required"),
  Fx=:date,
  Fy=:value,
  Fgroup=:variable,
  graphname::String = error("graph name is required"),
  yaxisname::String = error("y-axis name is required"))

  outputfolder::String = PARAM[:figurepath]

  panel = @vlplot() +
    @vlplot( #first part contains the line and circle dots
      data = series,
      width=styleidx[:width],
      height=styleidx[:height],
      title=nothing,
      mark={:line,
        clip=true,
        size=styleidx[:linesize]},
      x={Fx,
        axis={title=styleidx[:xaxisname],
          #formatType = "date",
          #format= "yyy",
          grid=false,
          titleFontSize=styleidx[:textsize],
          tickCount=11,
          titleFontWeight=:normal,
          labels=true},
          #=labelOpacity = {
            condition = {test = "datum.$(Fx) % 10 == 0",
              value = 0.0},
            value = 1.0}},=#
        type=:temporal,
        #scale={domain=styleidx[:xdomain]}
        },
      y={Fy,
        axis={title=yaxisname,
            grid=false,
            tickCount=4,
            titleFontSize=styleidx[:textsize],
            titleFontWeight=:normal},
          type=:quantitative,
          scale={zero=false,
            domain=styleidx[:ydomain]}
        },
      color={Fgroup,
        type=:nominal}
      )

  save("$outputfolder\\$graphname.pdf", panel)
end

#read in the appropriate file
function readprocessbetadata(;
  inbinstream::Function = IN_BIN_STREAM,
  outbinstream::Function = OUT_BIN_STREAM)

  prelimbetapath::String = "$(PARAM[:preliminarypath])\\$(PARAM[:prelimbetadata])"
  wrdsdateformat::DateFormat = PARAM[:wrdsdateformat]

  local beta::DataFrame
  local prelimbeta

  if PARAM[:refreshpreliminarybeta]
    beta = CSV.read("$(prelimbetapath).csv") |> DataFrame
    rename!(beta, (cleanname).(names(beta)))
    beta.date = (i::Int->Date("$i", wrdsdateformat)).(beta.date)
    outbinstream("$(prelimbetapath).$BIN_EXTENSION", beta)
  else
    beta = inbinstream("$(prelimbetapath).$BIN_EXTENSION")
  end

  return beta
end

function aggregatebetadata(beta::AbstractDataFrame;
    aggfields::Vector{Symbol} = [:bmkt, :bsmb, :bhml, :bumd])
  local series::DataFrame
  local seriesabs::DataFrame

  #create the aggregated dataframe by hand
  sbeta::SubDataFrame = view(beta, :, [:date; aggfields; :mktcap])
  series = DataFrame(date = unique(beta.date))
  series.bmkt = Vector{MFloat64}(undef, size(series,1))
  series.bsmb = similar(series.bmkt)
  series.bhml = similar(series.bmkt)
  series.bumd = similar(series.bmkt)
  series.absbmkt = similar(series.bmkt)
  series.absbsmb = similar(series.bmkt)
  series.absbhml = similar(series.bmkt)
  series.absbumd = similar(series.bmkt)

  seriesrows::Dict = Dict(r.date=>r for r ∈ eachrow(series))

  ssbetas = groupby(sbeta, :date)
  Threads.@threads for i::Int ∈ 1:length(ssbetas)
    ssbeta::SubDataFrame = ssbetas[i]
    r::DataFrameRow = seriesrows[ssbeta.date[1]]

    #aggregate each of the betas
    for f ∈ aggfields
      s3beta = view(ssbeta, completecases(ssbeta[!, [f, :mktcap]]), [f, :mktcap])
      (size(s3beta,1) == 0) && continue
      mktcap::Float64 = sum(s3beta.mktcap)
      r[f] = sum(s3beta.mktcap .* s3beta[!,f]) / mktcap
      r[Symbol(:abs, string(f))] = sum(s3beta.mktcap .* (abs).(s3beta[!,f])) / mktcap
    end
  end

  println(describe(series))
  #create a long series for each of the two graphs
  templong1 = stack(series[!,[:date; aggfields]], Not([:date]))

  #now create a long df for  the absolute values
  select!(series, Not(aggfields))
  rename!(series, (f::Symbol->Symbol(:abs, string(f))=>f).(aggfields))
  templong2 = stack(series[!,names(series)], Not([:date]), value_name=:valueabs)

  #join the two long dfs
  series = innerjoin(templong1, templong2, on=[:date, :variable])
  series = series[completecases(series),:]

  return series
end

#averages teh betas across permnos
function aggregatebetadataequalweighted(beta::AbstractDataFrame)
  local series::DataFrame
  local seriesabs::DataFrame

  sbeta::SubDataFrame = view(beta, :, [:date, :bmkt, :bsmb, :bhml, :bumd])
  println(describe(sbeta))
  aggfunc(v::AbstractVector{<:MFloat64}) = mean(skipmissing(v))
  series = aggregate(groupby(sbeta, :date), aggfunc)

  #fix the names
  rename!(series, (s::Symbol-> s=>nofunction(s, "_aggfunc")).(names(series)))

  #melt to graph
  series = stack(series, Not([:date]))

  #make the absolue value version
  aggfuncabs(v::AbstractVector{<:MFloat64}) = mean(skipmissing((abs).(v)))
  seriesabs = aggregate(groupby(sbeta, :date), aggfuncabs)
  rename!(seriesabs, (s::Symbol-> s=>nofunction(s, "_aggfuncabs")).(names(seriesabs)))
  seriesabs = stack(seriesabs, Not([:date]), value_name=:valueabs)

  series = innerjoin(series, seriesabs, on=[:date, :variable])

  return series
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
export SuperpositionMeasure

"""
    struct SuperpositionMeasure{X,NT} <: AbstractMeasure
        components :: NT
    end

Superposition of measures is analogous to mixture distributions, but (because
measures need not be normalized) requires no scaling.

The superposition of two measures μ and ν can be more concisely written as μ + ν.

Superposition measures satisfy

    basemeasure(μ + ν) == basemeasure(μ) + basemeasure(ν)
"""
struct SuperpositionMeasure{NT} <: AbstractMeasure
    components :: NT
end

# SuperpositionMeasure(ms :: AbstractMeasure...) = SuperpositionMeasure{X,length(ms)}(ms)

# SuperpositionMeasure(m::NTuple{N, Measure{X}}) where {N,X} = SuperpositionMeasure(m...)

# Base.length(::SuperpositionMeasure{X,N}) where {X,N} = N

# function Base.:+(μ::SuperpositionMeasure{N1}, ν::SuperpositionMeasure{N2}) where {N1,N2}
#     components = (μ.components..., ν.components...)
#     SuperpositionMeasure{X, N1+N2}(components)
# end

# function Base.:+(μ::AbstractMeasure, ν::SuperpositionMeasure{X,N}) where {X,N}
#     components = (μ, ν.components...)
#     SuperpositionMeasure{X,N+1}(components)
# end

# function Base.:+(μ::SuperpositionMeasure{X,N}, ν::AbstractMeasure) where {X,N}
#     components = (μ.components..., ν)
#     SuperpositionMeasure{X,N+1}(components)
# end

function Base.:+(μ::AbstractMeasure, ν::AbstractMeasure)
    components = (μ, ν)
    SuperpositionMeasure(components)
end

logdensity(μ::SuperpositionMeasure, x) = logsumexp((logdensity(m,x) for m in μ.components))

basemeasure(μ::SuperpositionMeasure) = SuperpositionMeasure(map(basemeasure, μ.components))

# TODO: Fix `rand` method (this one is wrong)
# function Base.rand(μ::SuperpositionMeasure{X,N}) where {X,N}
#     return rand(rand(μ.components))
# end


text:
#=
Created on Sunday 19 April 2020
Last update: Sunday 26 April 2020

@author: Michiel Stock
michielfmstock@gmail.com

Illustration of the algorithms for solving TSP.

Currently implements:
    - nearest neighbors
    - greedy
    - random insertion
    - hill climbing
    - simulated annealing
    - tabu search
=#

using STMO
using STMO.TSP

using Plots

tsp = totoro_tsp()
name = "totoro"

n = length(tsp)

# NEAREST NEIGHBORS

"""
    nearestneighbors(tsp::TravelingSalesmanProblem; start::Int)

Solves the TSP using the nearest neighbors. Provide a starting city using `start`.
If none is provide, a random one is chosen.
"""
function nearestneighbors(tsp::TravelingSalesmanProblem; start::Int)
    cities_to_add = Set(cities(tsp))
    delete!(cities_to_add, start)
    tour = [start]
    current = start
    cost = 0.0
    while length(cities_to_add) > 0
        # find closest city not in tour
        c, next = minimum((dist(tsp, current, n), n) for n in cities_to_add)
        push!(tour, next)
        delete!(cities_to_add, next)
        cost += c
        current = next
    end
    return tour, cost
end

nearestneighbors(tsp::TravelingSalesmanProblem) = nearestneighbors(tsp, start=rand(cities(tsp)))

"""
    bestnearestneighbors(tsp::TravelingSalesmanProblem;
                        ntry::Union{Nothing,Int}=nothing)

Chooses the best nearest neighbor solution over all cities. If `ntry` is provided,
a random number of cities is tried. This is done if searching all cities is too
involved.
"""
function bestnearestneighbors(tsp::TravelingSalesmanProblem;
                        ntry::Union{Nothing,Int}=nothing)
    if ntry isa Nothing
        cities_to_try = cities(tsp)
    else
        cities_to_try = rand(cities(tsp), ntry)
    end
    best_cost = Inf64
    best_tour = [1, 2]
    for start in cities_to_try
        tour, cost = nearestneighbors(tsp, start=start)
        if cost < best_cost
            best_cost = cost
            best_tour = tour
        end
    end
    return best_tour, best_cost
end

tour_nn, cost_nn = nearestneighbors(tsp)
tour_nnbest, cost_nnbest = bestnearestneighbors(tsp, ntry=100)  # try 100 starts

p_nn = plot_cities(tsp, markersize=1)
plot_tour!(tsp, tour_nn, label="random")
plot_tour!(tsp, tour_nnbest, label="best", color=mygreen)
title!("Nearest Neighbor\n cost = $cost_nnbest")

# GREEDY

using DataStructures

"""
    greedy(tsp::TravelingSalesmanProblem)

Uses the greedy algorithm to solve the TSP.
"""
function greedy(tsp::TravelingSalesmanProblem)
    usf = DisjointSets(cities(tsp))
    times_added = Dict(c=>0 for c in cities(tsp))
    n = length(tsp)
    edges = [(dist(tsp, i, j), i, j) for i in 1:n-1 for j in (i+1):n]
    cost = 0.0
    selected_edges = Tuple{Int,Int}[]
    for (c, i, j) in sort!(edges)
        if !in_same_set(usf, i, j) && times_added[i] < 2 && times_added[j] < 2
            union!(usf, i, j)
            push!(selected_edges, (i, j))
            times_added[i] += 1
            times_added[j] += 1
            cost += c
            length(selected_edges) == n  && break
        end
    end
    # close the loop
    current = findfirst(c->times_added[c]==1, cities(tsp))
    # knit edges in a tour
    tour = [current]
    sizehint!(tour, n)
    used = Set([current])
    while length(tour) < n
        for (i, j) in selected_edges
            if i==current && j ∉ used
                push!(tour, j)
                push!(used, j)
                current = j
                break
            elseif j==current && i ∉ used
                push!(tour, i)
                push!(used, i)
                current = i
                break
            end
        end
    end
    return tour, cost
end

tour_greedy, cost_greedy = greedy(tsp)

p_greedy = plot_cities(tsp, markersize=1)
plot_tour!(tsp, tour_greedy)
title!("Greedy\n cost = $cost_greedy")

"""
    insertion(tsp::TravelingSalesmanProblem)

Randomly inserts cities in a tour at a place where it has the lowest cost.
"""
function insertion(tsp::TravelingSalesmanProblem)
    n = length(tsp)
    start = rand(cities(tsp))
    tour = [start]
    sizehint!(tour, n)
    cities_to_try = Set(cities(tsp))
    delete!(cities_to_try, start)
    cost = 0.0
    while length(tour) < n
        city = rand(cities_to_try)
        connection_cost = Inf64
        best_pos = 0
        # find cheapest place to insert city
        for pos in 1:length(tour)
            ci = (pos > 1) ? tour[pos-1] : tour[end]
            cj = tour[pos]
            # cost of connecting
            Δc = dist(tsp, ci, city) + dist(tsp, city, cj) - dist(tsp, ci, cj)
            if Δc < connection_cost
                connection_cost = Δc
                best_pos = pos
            end
        end
        insert!(tour, best_pos, city)
        delete!(cities_to_try, city)
        cost += connection_cost
    end
    return tour, computecost(tsp, tour)
end

tour_insertion, cost_insertion = insertion(tsp)

p_insertion = plot_cities(tsp, markersize=1)
plot_tour!(tsp, tour_insertion)
title!("Insertion (random)\ncost=$cost_insertion")

# 2-opt

"""
    hillclimbing!(tsp, tour; verbose=false, maxitter=Inf)

Uses hill climbing to improve a tour by finding each iteration the best path
between two cities to flip.
"""
function hillclimbing!(tsp, tour; verbose=false, maxitter=Inf)
    n = length(tsp)
    improved = true
    cost = computecost(tsp, tour)
    iter = 0
    costs = [cost]
    while improved
        best_Δc = 0.0  # any change is sufficient to continue
        improved = false
        best_i, best_j = 0, 0
        for i in 1:n-2
            for j in (i+1):n
                Δc = deltaflipcost(tsp, tour, i, j)
                if Δc < best_Δc
                    best_Δc = Δc
                    best_i, best_j = i, j
                    improved = true
                end
            end
        end
        !improved && break
        flip!(tour, best_i, best_j)
        iter += 1
        cost += best_Δc
        push!(costs, cost)
        verbose && println("Iteration $iter: Δcost = $best_Δc (flipped $best_i, $best_j)")
        iter > maxitter && break
    end
    println("converged in $iter steps")
    return tour, cost, costs
end

"""
    hillclimbing(tsp; verbose=false, maxitter=Inf)

Uses hill climbing to improve a tour by finding each iteration the best path
between two cities to flip. Starts from the given order of the cities.
"""
hillclimbing(tsp; kwargs...) = hillclimbing!(tsp, collect(1:length(tsp)); kwargs...)

tour_hc, cost_hc, costs_hc = hillclimbing(tsp, verbose=true)

p_hc = plot_cities(tsp, markersize=1)
plot_tour!(tsp, tour_hc)
title!("Hill Climbing\ncost=$cost_hc")

# SIMULATED ANNEALING

"""
    simulatedannealing!tsp, tour;
                    Tmax, Tmin, r, kT::Int, verbose=false)

Uses simulated annealing to improve a tour by finding each iteration the best path
between two cities to flip.
"""
function simulatedannealing!(tsp, tour;
                Tmax, Tmin, r, kT::Int, verbose=false)
    n = length(tsp)
    improved = true
    cost = computecost(tsp, tour)
    T = Tmax
    n_steps = (log(Tmin) - log(Tmax)) / log(r) |> ceil |> Int
    costs = Vector{typeof(cost)}(undef, n_steps)
    iter = 0
    while T > Tmin
        iter += 1
        acc = 0
        for k in 1:kT
            # choose two cities to swap
            i, j = rand(1:n, 2)
            # improvement
            Δc = deltaflipcost(tsp, tour, i, j)
            if Δc < 0.0 || rand() < exp(- Δc / T)
                cost += Δc
                flip!(tour, i, j)
                acc += 1
            end
        end
        verbose && println("T = $T : cost = $cost (acc % = $(acc / kT))")
        costs[iter] = cost
        T *= r
    end
    return tour, cost, costs
end

simulatedannealing(tsp; kwargs...) = simulatedannealing!(tsp,
                                            collect(1:length(tsp)); kwargs...)

tour_sa = collect(1:n)

tour_sa, cost_sa, costs_sa = simulatedannealing!(tsp, tour_sa;
                Tmax=1e7, Tmin=1e-4, r=0.975, kT=100000, verbose=true)

p_sa = plot_cities(tsp, markersize=1)
plot_tour!(tsp, tour_sa)
title!("Simulated Annealing\ncost=$cost_sa")

# TABU SEARCH

"""
    tabusearch!(tsp::TravelingSalesmanProblem, tour; ntabu::Int, niter::Int,
                            verbose=false)

Improves a tour by iteratively performing the best local improvement, similarly
to `hillclimbing`. In tabu search however, after a position of the tour is
modified, it is blocked for `ntabu` steps. This behaviour is meant to escape
local minima.
"""
function tabusearch!(tsp::TravelingSalesmanProblem, tour; ntabu::Int, niter::Int,
                            verbose=false)
    n = length(tsp)
    cost = computecost(tsp, tour)
    # positions in the tour that are tabood
    tabood = zeros(Int, n)
    costs = Vector(undef, niter)
    for iter in 1:niter
        best_Δc = Inf64
        best_i, best_j = 0, 0
        for i in 1:n-2
            # if i is tabood skip it
            tabood[i] > iter && continue
            for j in (i+1):n
                # if j is tabood skip it
                tabood[j] > iter && continue
                Δc = deltaflipcost(tsp, tour, i, j)
                if Δc < best_Δc
                    best_Δc = Δc
                    best_i, best_j = i, j
                end
            end
        end
        flip!(tour, best_i, best_j)
        cost += best_Δc
        # tabu these for ntabu steps
        tabood[best_i] = iter + ntabu
        tabood[best_j] = iter + ntabu
        costs[iter] = cost
        verbose && println("Iteration $iter: Δcost = $best_Δc (flipped $best_i, $best_j)")
    end
    return tour, cost, costs
end

tabusearch(tsp; kwargs...) = tabusearch!(tsp, collect(1:length(tsp)); kwargs...)

tour_tabu, cost_tabu, costs_tabu = tabusearch(tsp, ntabu=20,
                niter=2_000, verbose=true)

p_tabu = plot_cities(tsp, markersize=1)
plot_tour!(tsp, tour_tabu)
title!("Tabu Search\ncost=$cost_tabu")

# SUMMARIZE
# ---------

plot(p_nn, p_greedy, p_insertion, p_hc, p_sa, p_tabu, size=(1000, 1000))
savefig("figures/tsp_algorithms_$name.png")

plot(costs_hc, label="hill climbing")
plot!(costs_sa, label="simulated annealing")
plot!(costs_tabu, label="tabu search")
xlabel!("Iteration")
ylabel!("Tour cost")
savefig("figures/tsp_convergence_$name.png")


text:
"""
    ITouchable(nexid)

An interface for touchable widgets.
"""
struct ITouchable <: AbstractINextion
    _nid::NexID

    function ITouchable(nid::NexID)
        new(nid)
    end
end


text:
"""
    InhomogPoisson(λ)

Inhomogenous Poisson process with intensity function `λ(t)`.
See also `ThinningAlg`.
"""
struct InhomogPoisson{T<:Function}
    λ::T
end

abstract type PoissonSampler
end

"""
    ThinningAlg(λmax)

Sampling method for `InhomogPoisson` by the 'thinning' algorithm.

#### Examples:

```
sample(ThinningAlg(λmax), T, InhomogPoisson(λ))
```
"""
struct ThinningAlg <: PoissonSampler
    λmax::Float64
end

"""
    sample(::Thinning, T, P::InhomogPoisson) -> tt
"""
function sample(method::ThinningAlg, T, P::InhomogPoisson)
    λmax = method.λmax
    t = 0.0
    tt = zeros(0)
    while t <= T
        t = t - log(rand())/λmax
        if rand() ≤ P.λ(t)/λmax
            push!(tt, t)
        end
    end
    tt
end


text:
"""
	write_newick([file::String,] tree::Tree)

Write `tree` as a newick string in `file`.
  If `file` is not provided, return the newick string.
"""
write_newick(file::String, tree::Tree, mode="w") = write_newick(file, tree.root, mode)
write_newick(tree::Tree) = write_newick(tree.root)

"""
	write_newick([file::String,] root::TreeNode)
"""
function write_newick(file::String, root::TreeNode, mode = "w")
	out = write_newick(root)
	open(file, mode) do f
		write(f, out)
	end

	return nothing
end
write_newick(root::TreeNode) = write_newick!("", root)*";"

"""
"""
function write_newick!(s::String, root::TreeNode)
	if !isempty(root.child)
		s *= '('
		temp = sort(root.child, by=x->length(POTleaves(x)))
		for c in temp
			s = write_newick!(s, c)
			s *= ','
		end
		s = s[1:end-1] # Removing trailing ','
		s *= ')'
	end
	s *= root.label
	if !ismissing(root.tau)
		s *= ':'
		s *= string(root.tau)
	end
	if root.isroot && ismissing(root.tau)
		s *= ":0"
	end

	return s
end

max_stars_repo_name: PierreBarrat/TreeTools

max_issues_repo_name: PierreBarrat/TreeTools

max_forks_repo_name: PierreBarrat/TreeTools


text:
module TestNullableArray

using Base.Test
using CategoricalArrays
using NullableArrays
using CategoricalArrays: DefaultRefType
using Compat

typealias String Compat.ASCIIString

# == currently throws an error for Nullables
(==) = isequal

for isordered in (false, true)
    for R in (CategoricalArrays.DefaultRefType, UInt8, UInt, Int8, Int)
        # Vector with no null values
        for a in (["b", "a", "b"],
                  Nullable{String}["b", "a", "b"],
                  NullableArray(["b", "a", "b"]))
            x = NullableCategoricalVector{String, R}(a, ordered=isordered)
            na = eltype(a) <: Nullable ? a : convert(Array{Nullable{String}}, a)

            @test x == na
            @test ordered(x) === isordered
            @test levels(x) == sort(map(get, unique(na)))
            @test size(x) === (3,)
            @test length(x) === 3

            @test convert(NullableCategoricalArray, x) === x
            @test convert(NullableCategoricalArray{String}, x) === x
            @test convert(NullableCategoricalArray{String, 1}, x) === x
            @test convert(NullableCategoricalArray{String, 1, R}, x) === x
            @test convert(NullableCategoricalArray{String, 1, DefaultRefType}, x) == x
            @test convert(NullableCategoricalArray{String, 1, UInt8}, x) == x

            @test convert(NullableCategoricalVector, x) === x
            @test convert(NullableCategoricalVector{String}, x) === x
            @test convert(NullableCategoricalVector{String, R}, x) === x
            @test convert(NullableCategoricalVector{String, DefaultRefType}, x) == x
            @test convert(NullableCategoricalVector{String, UInt8}, x) == x

            for y in (NullableCategoricalArray(x, ordered=isordered),
                      NullableCategoricalArray{String}(x, ordered=isordered),
                      NullableCategoricalArray{String, 1}(x, ordered=isordered),
                      NullableCategoricalArray{String, 1, R}(x, ordered=isordered),
                      NullableCategoricalArray{String, 1, DefaultRefType}(x, ordered=isordered),
                      NullableCategoricalArray{String, 1, UInt8}(x, ordered=isordered),
                      NullableCategoricalVector(x, ordered=isordered),
                      NullableCategoricalVector{String}(x, ordered=isordered),
                      NullableCategoricalVector{String, R}(x, ordered=isordered),
                      NullableCategoricalVector{String, DefaultRefType}(x, ordered=isordered),
                      NullableCategoricalVector{String, UInt8}(x, ordered=isordered),
                      categorical(x, ordered=isordered),
                      categorical(x, false, ordered=isordered),
                      categorical(x, true, ordered=isordered))
                @test isa(y, NullableCategoricalVector{String})
                @test ordered(y) === isordered
                @test y == x
                @test y !== x
                @test y.refs !== x.refs
                @test y.pool !== x.pool
            end

            @test convert(NullableCategoricalArray, a) == x
            @test convert(NullableCategoricalArray{String}, a) == x
            @test convert(NullableCategoricalArray{String, 1}, a) == x
            @test convert(NullableCategoricalArray{String, 1, R}, a) == x
            @test convert(NullableCategoricalArray{String, 1, DefaultRefType}, a) == x
            @test convert(NullableCategoricalArray{String, 1, UInt8}, a) == x

            @test convert(NullableCategoricalVector, a) == x
            @test convert(NullableCategoricalVector{String}, a) == x
            @test convert(NullableCategoricalVector{String, R}, a) == x
            @test convert(NullableCategoricalVector{String, DefaultRefType}, a) == x
            @test convert(NullableCategoricalVector{String, UInt8}, a) == x

            @test NullableCategoricalArray{String}(a, ordered=isordered) == x
            @test NullableCategoricalArray{String, 1}(a, ordered=isordered) == x
            @test NullableCategoricalArray{String, 1, R}(a, ordered=isordered) == x
            @test NullableCategoricalArray{String, 1, DefaultRefType}(a, ordered=isordered) == x
            @test NullableCategoricalArray{String, 1, UInt8}(a, ordered=isordered) == x

            @test NullableCategoricalVector(a, ordered=isordered) == x
            @test NullableCategoricalVector{String}(a, ordered=isordered) == x
            @test NullableCategoricalVector{String, R}(a, ordered=isordered) == x
            @test NullableCategoricalVector{String, DefaultRefType}(a, ordered=isordered) == x
            @test NullableCategoricalVector{String, UInt8}(a, ordered=isordered) == x

            for (y, R1, R2, compact) in ((a, DefaultRefType, UInt8, true),
                                         (a, DefaultRefType, DefaultRefType, false),
                                         (x, R, UInt8, true),
                                         (x, R, R, false))
                x2 = categorical(y, ordered=isordered)
                @test x2 == y
                if eltype(y) <: Nullable
                    @test isa(x2, NullableCategoricalVector{String, R1})
                else
                    @test isa(x2, CategoricalVector{String, R1})
                end
                @test ordered(x2) === isordered

                x2 = categorical(y, compact, ordered=isordered)
                @test x2 == y
                if eltype(y) <: Nullable
                    @test isa(x2, NullableCategoricalVector{String, R2})
                else
                    @test isa(x2, CategoricalVector{String, R2})
                end
                @test ordered(x2) === isordered
            end

            x2 = compact(x)
            @test x2 == x
            @test isa(x2, NullableCategoricalVector{String, UInt8})
            @test ordered(x2) === ordered(x)
            @test levels(x2) == levels(x)

            x2 = copy(x)
            @test x2 == x
            @test ordered(x2) === ordered(x)
            @test typeof(x2) === typeof(x)
            @test levels(x2) == levels(x)

            if !ordered(x)
                @test ordered!(x, true) === true
            end
            @test get(x[1] > x[2])
            @test get(x[3] > x[2])

            @test ordered!(x, false) === false
            @test ordered(x) === false
            @test_throws Exception x[1] > x[2]
            @test_throws Exception x[3] > x[2]

            @test x[1] === Nullable(x.pool.valindex[1])
            @test x[2] === Nullable(x.pool.valindex[2])
            @test x[3] === Nullable(x.pool.valindex[1])
            @test_throws BoundsError x[4]

            @test x[1:2] == Nullable{String}["b", "a"]
            @test typeof(x[1:2]) === typeof(x)

            x[1] = x[2]
            @test x[1] === Nullable(x.pool.valindex[2])
            @test x[2] === Nullable(x.pool.valindex[2])
            @test x[3] === Nullable(x.pool.valindex[1])

            x[3] = "c"
            @test x[1] === Nullable(x.pool.valindex[2])
            @test x[2] === Nullable(x.pool.valindex[2])
            @test x[3] === Nullable(x.pool.valindex[3])
            @test levels(x) == ["a", "b", "c"]

            x[2:3] = "b"
            @test x[1] === Nullable(x.pool.valindex[2])
            @test x[2] === Nullable(x.pool.valindex[1])
            @test x[3] === Nullable(x.pool.valindex[1])
            @test levels(x) == ["a", "b", "c"]

            droplevels!(x) == ["a", "b"]
            @test x[1] === Nullable(x.pool.valindex[1])
            @test x[2] === Nullable(x.pool.valindex[2])
            @test x[3] === Nullable(x.pool.valindex[2])
            @test levels(x) == ["a", "b"]

            levels!(x, ["b", "a"]) == ["b", "a"]
            @test x[1] === Nullable(x.pool.valindex[1])
            @test x[2] === Nullable(x.pool.valindex[2])
            @test x[3] === Nullable(x.pool.valindex[2])
            @test levels(x) == ["b", "a"]

            @test_throws ArgumentError levels!(x, ["a"])
            @test_throws ArgumentError levels!(x, ["e", "b"])
            @test_throws ArgumentError levels!(x, ["e", "a", "b", "a"])

            @test levels!(x, ["e", "a", "b"]) == ["e", "a", "b"]
            @test x[1] === Nullable(x.pool.valindex[1])
            @test x[2] === Nullable(x.pool.valindex[2])
            @test x[3] === Nullable(x.pool.valindex[2])
            @test levels(x) == ["e", "a", "b"]

            x[1] = "c"
            @test x[1] === Nullable(x.pool.valindex[4])
            @test x[2] === Nullable(x.pool.valindex[2])
            @test x[3] === Nullable(x.pool.valindex[2])
            @test levels(x) == ["e", "a", "b", "c"]

            @test_throws ArgumentError levels!(x, ["e", "c"])
            @test levels!(x, ["e", "c"], nullok=true) == ["e", "c"]
            @test x[1] === Nullable(x.pool.valindex[2])
            @test x[2] === eltype(x)()
            @test x[3] === eltype(x)()
            @test levels(x) == ["e", "c"]

            push!(x, "e")
            @test length(x) == 4
            @test isequal(x, NullableArray(["c", "", "", "e"], [false, true, true, false]))
            @test levels(x) == ["e", "c"]

            push!(x, "zz")
            @test length(x) == 5
            @test isequal(x, NullableArray(["c", "", "", "e", "zz"], [false, true, true, false, false]))
            @test levels(x) == ["e", "c", "zz"]

            push!(x, x[1])
            @test length(x) == 6
            @test isequal(x, NullableArray(["c", "", "", "e", "zz", "c"], [false, true, true, false, false, false]))
            @test levels(x) == ["e", "c", "zz"]

            push!(x, eltype(x)())
            @test length(x) == 7
            @test isequal(x, NullableArray(["c", "", "", "e", "zz", "c", ""], [false, true, true, false, false, false, true]))
            @test isnull(x[end])
            @test levels(x) == ["e", "c", "zz"]

            append!(x, x)
            @test isequal(x, NullableArray(["c", "", "", "e", "zz", "c", "", "c", "", "", "e", "zz", "c", ""], [false, true, true, false, false, false, true, false, true, true, false, false, false, true]))
            @test levels(x) == ["e", "c", "zz"]
            @test ordered(x) === false
            @test length(x) == 14

            b = ["z","y","x"]
            y = NullableCategoricalVector{String, R}(b)
            append!(x, y)
            @test length(x) == 17
            @test ordered(x) === false
            @test levels(x) == ["e", "c", "zz", "x", "y", "z"]
            @test isequal(x, NullableArray(["c", "", "", "e", "zz", "c", "", "c", "", "", "e", "zz", "c", "", "z", "y", "x"], [false, true, true, false, false, false, true, false, true, true, false, false, false, true, false, false, false]))

            empty!(x)
            @test ordered(x) === false
            @test length(x) == 0
            @test levels(x) == ["e", "c", "zz", "x", "y", "z"]
        end

        # Vector with null values
        for a in (Nullable{String}["a", "b", Nullable()],
                  NullableArray(Nullable{String}["a", "b", Nullable()]))
            x = NullableCategoricalVector{String, R}(a, ordered=isordered)

            @test x == a
            @test levels(x) == map(get, filter(x->!isnull(x), unique(a)))
            @test size(x) === (3,)
            @test length(x) === 3

            @test convert(NullableCategoricalArray, x) === x
            @test convert(NullableCategoricalArray{String}, x) === x
            @test convert(NullableCategoricalArray{String, 1}, x) === x
            @test convert(NullableCategoricalArray{String, 1, R}, x) === x
            @test convert(NullableCategoricalArray{String, 1, DefaultRefType}, x) == x
            @test convert(NullableCategoricalArray{String, 1, UInt8}, x) == x

            @test convert(NullableCategoricalVector, x) === x
            @test convert(NullableCategoricalVector{String}, x) === x
            @test convert(NullableCategoricalVector{String, R}, x) === x
            @test convert(NullableCategoricalVector{String, DefaultRefType}, x) == x
            @test convert(NullableCategoricalVector{String, UInt8}, x) == x

            for y in (NullableCategoricalArray(x, ordered=isordered),
                      NullableCategoricalArray{String}(x, ordered=isordered),
                      NullableCategoricalArray{String, 1}(x, ordered=isordered),
                      NullableCategoricalArray{String, 1, R}(x, ordered=isordered),
                      NullableCategoricalArray{String, 1, DefaultRefType}(x, ordered=isordered),
                      NullableCategoricalArray{String, 1, UInt8}(x, ordered=isordered),
                      NullableCategoricalVector(x, ordered=isordered),
                      NullableCategoricalVector{String}(x, ordered=isordered),
                      NullableCategoricalVector{String, R}(x, ordered=isordered),
                      NullableCategoricalVector{String, DefaultRefType}(x, ordered=isordered),
                      NullableCategoricalVector{String, UInt8}(x, ordered=isordered),
                      categorical(x, ordered=isordered),
                      categorical(x, false, ordered=isordered),
                      categorical(x, true, ordered=isordered))
                @test isa(y, NullableCategoricalVector{String})
                @test ordered(y) === isordered
                @test y == x
                @test y !== x
                @test y.refs !== x.refs
                @test y.pool !== x.pool
            end

            @test convert(NullableCategoricalArray, a) == x
            @test convert(NullableCategoricalArray{String}, a) == x
            @test convert(NullableCategoricalArray{String, 1}, a) == x
            @test convert(NullableCategoricalArray{String, 1, R}, a) == x
            @test convert(NullableCategoricalArray{String, 1, DefaultRefType}, a) == x
            @test convert(NullableCategoricalArray{String, 1, UInt8}, a) == x

            @test convert(NullableCategoricalVector, a) == x
            @test convert(NullableCategoricalVector{String}, a) == x
            @test convert(NullableCategoricalVector{String, R}, a) == x
            @test convert(NullableCategoricalVector{String, DefaultRefType}, a) == x
            @test convert(NullableCategoricalVector{String, UInt8}, a) == x

            @test NullableCategoricalArray{String}(a, ordered=isordered) == x
            @test NullableCategoricalArray{String, 1}(a, ordered=isordered) == x
            @test NullableCategoricalArray{String, 1, R}(a, ordered=isordered) == x
            @test NullableCategoricalArray{String, 1, DefaultRefType}(a, ordered=isordered) == x
            @test NullableCategoricalArray{String, 1, UInt8}(a, ordered=isordered) == x

            @test NullableCategoricalVector(a, ordered=isordered) == x
            @test NullableCategoricalVector{String}(a, ordered=isordered) == x
            @test NullableCategoricalVector{String, R}(a, ordered=isordered) == x
            @test NullableCategoricalVector{String, DefaultRefType}(a, ordered=isordered) == x
            @test NullableCategoricalVector{String, UInt8}(a, ordered=isordered) == x

            for (y, R1, R2, compact) in ((a, DefaultRefType, UInt8, true),
                                         (a, DefaultRefType, DefaultRefType, false),
                                         (x, R, UInt8, true),
                                         (x, R, R, false))
                x2 = categorical(y, ordered=isordered)
                @test x2 == y
                @test isa(x2, NullableCategoricalVector{String, R1})
                @test ordered(x2) === isordered

                x2 = categorical(y, compact, ordered=isordered)
                @test x2 == y
                @test isa(x2, NullableCategoricalVector{String, R2})
                @test ordered(x2) === isordered
            end

            x2 = compact(x)
            @test x2 == x
            @test isa(x2, NullableCategoricalVector{String, UInt8})
            @test levels(x2) == levels(x)

            x2 = copy(x)
            @test x2 == x
            @test typeof(x2) === typeof(x)
            @test levels(x2) == levels(x)

            @test x[1] === Nullable(x.pool.valindex[1])
            @test x[2] === Nullable(x.pool.valindex[2])
            @test x[3] === eltype(x)()
            @test_throws BoundsError x[4]

            @test x[1:2] == Nullable{String}["a", "b"]
            @test typeof(x[1:2]) === typeof(x)

            @test x[2:3] == Nullable{String}["b", Nullable()]
            @test typeof(x[2:3]) === typeof(x)

            x[1] = "b"
            @test x[1] === Nullable(x.pool.valindex[2])
            @test x[2] === Nullable(x.pool.valindex[2])
            @test x[3] === eltype(x)()

            x[3] = "c"
            @test x[1] === Nullable(x.pool.valindex[2])
            @test x[2] === Nullable(x.pool.valindex[2])
            @test x[3] === Nullable(x.pool.valindex[3])
            @test levels(x) == ["a", "b", "c"]

            x[1] = Nullable()
            @test x[1] === eltype(x)()
            @test x[2] === Nullable(x.pool.valindex[2])
            @test x[3] === Nullable(x.pool.valindex[3])
            @test levels(x) == ["a", "b", "c"]

            x[2:3] = Nullable()
            @test x[1] === eltype(x)()
            @test x[2] === eltype(x)()
            @test x[3] === eltype(x)()
            @test levels(x) == ["a", "b", "c"]
        end

        # Vector created from range (i.e. non-Array AbstractArray),
        # direct conversion to a vector with different eltype
        a = 0.0:0.5:1.5
        x = NullableCategoricalVector{Float64, R}(a, ordered=isordered)

        @test x == map(Nullable, a)
        @test ordered(x) === isordered
        @test levels(x) == unique(a)
        @test size(x) === (4,)
        @test length(x) === 4

        @test convert(NullableCategoricalArray, x) === x
        @test convert(NullableCategoricalArray{Float64}, x) === x
        @test convert(NullableCategoricalArray{Float64, 1}, x) === x
        @test convert(NullableCategoricalArray{Float64, 1, R}, x) === x
        @test convert(NullableCategoricalArray{Float64, 1, DefaultRefType}, x) == x
        @test convert(NullableCategoricalArray{Float64, 1, UInt8}, x) == x

        @test convert(NullableCategoricalVector, x) === x
        @test convert(NullableCategoricalVector{Float64}, x) === x
        @test convert(NullableCategoricalVector{Float64, R}, x) === x
        @test convert(NullableCategoricalVector{Float64, DefaultRefType}, x) == x
        @test convert(NullableCategoricalVector{Float64, UInt8}, x) == x

        for y in (NullableCategoricalArray(x, ordered=isordered),
                  NullableCategoricalArray{Float64}(x, ordered=isordered),
                  NullableCategoricalArray{Float64, 1}(x, ordered=isordered),
                  NullableCategoricalArray{Float64, 1, R}(x, ordered=isordered),
                  NullableCategoricalArray{Float64, 1, DefaultRefType}(x, ordered=isordered),
                  NullableCategoricalArray{Float64, 1, UInt8}(x, ordered=isordered),
                  NullableCategoricalVector(x, ordered=isordered),
                  NullableCategoricalVector{Float64}(x, ordered=isordered),
                  NullableCategoricalVector{Float64, R}(x, ordered=isordered),
                  NullableCategoricalVector{Float64, DefaultRefType}(x, ordered=isordered),
                  NullableCategoricalVector{Float64, UInt8}(x, ordered=isordered),
                  categorical(x, ordered=isordered),
                  categorical(x, false, ordered=isordered),
                  categorical(x, true, ordered=isordered))
            @test isa(y, NullableCategoricalVector{Float64})
            @test ordered(y) === isordered
            @test y == x
            @test y !== x
            @test y.refs !== x.refs
            @test y.pool !== x.pool
        end

        @test convert(NullableCategoricalArray, a) == x
        @test convert(NullableCategoricalArray{Float64}, a) == x
        @test convert(NullableCategoricalArray{Float32}, a) == x
        @test convert(NullableCategoricalArray{Float64, 1}, a) == x
        @test convert(NullableCategoricalArray{Float32, 1}, a) == x
        @test convert(NullableCategoricalArray{Float64, 1, R}, a) == x
        @test convert(NullableCategoricalArray{Float32, 1, R}, a) == x
        @test convert(NullableCategoricalArray{Float64, 1, DefaultRefType}, a) == x
        @test convert(NullableCategoricalArray{Float32, 1, DefaultRefType}, a) == x
        @test convert(NullableCategoricalArray{Float64, 1, UInt8}, a) == x
        @test convert(NullableCategoricalArray{Float32, 1, UInt8}, a) == x

        @test convert(NullableCategoricalVector, a) == x
        @test convert(NullableCategoricalVector{Float64}, a) == x
        @test convert(NullableCategoricalVector{Float32}, a) == x
        @test convert(NullableCategoricalVector{Float64, R}, a) == x
        @test convert(NullableCategoricalVector{Float32, R}, a) == x
        @test convert(NullableCategoricalVector{Float64, DefaultRefType}, a) == x
        @test convert(NullableCategoricalVector{Float32, DefaultRefType}, a) == x
        @test convert(NullableCategoricalVector{Float64, UInt8}, a) == x
        @test convert(NullableCategoricalVector{Float32, UInt8}, a) == x

        @test NullableCategoricalArray{Float64}(a, ordered=isordered) == x
        @test NullableCategoricalArray{Float32}(a, ordered=isordered) == x
        @test NullableCategoricalArray{Float64, 1}(a, ordered=isordered) == x
        @test NullableCategoricalArray{Float32, 1}(a, ordered=isordered) == x
        @test NullableCategoricalArray{Float64, 1, R}(a, ordered=isordered) == x
        @test NullableCategoricalArray{Float32, 1, R}(a, ordered=isordered) == x
        @test NullableCategoricalArray{Float64, 1, DefaultRefType}(a, ordered=isordered) == x
        @test NullableCategoricalArray{Float32, 1, DefaultRefType}(a, ordered=isordered) == x

        @test NullableCategoricalVector(a, ordered=isordered) == x
        @test NullableCategoricalVector{Float64}(a, ordered=isordered) == x
        @test NullableCategoricalVector{Float32}(a, ordered=isordered) == x
        @test NullableCategoricalVector{Float64, R}(a, ordered=isordered) == x
        @test NullableCategoricalVector{Float32, R}(a, ordered=isordered) == x
        @test NullableCategoricalVector{Float64, DefaultRefType}(a, ordered=isordered) == x
        @test NullableCategoricalVector{Float32, DefaultRefType}(a, ordered=isordered) == x

        for (y, R1, R2, compact) in ((a, DefaultRefType, UInt8, true),
                                     (a, DefaultRefType, DefaultRefType, false),
                                     (x, R, UInt8, true),
                                     (x, R, R, false))
            x2 = categorical(y, ordered=isordered)
            @test x2 == collect(y)
            if eltype(y) <: Nullable
                @test isa(x2, NullableCategoricalVector{Float64, R1})
            else
                @test isa(x2, CategoricalVector{Float64, R1})
            end
            @test ordered(x2) === isordered

            x2 = categorical(y, compact, ordered=isordered)
            @test x2 == collect(y)
            if eltype(y) <: Nullable
                @test isa(x2, NullableCategoricalVector{Float64, R2})
            else
                @test isa(x2, CategoricalVector{Float64, R2})
            end
            @test ordered(x2) === isordered
        end

        x2 = compact(x)
        @test x2 == x
        @test ordered(x2) === ordered(x)
        @test isa(x2, NullableCategoricalVector{Float64, UInt8})
        @test levels(x2) == levels(x)

        x2 = copy(x)
        @test x2 == x
        @test ordered(x2) === ordered(x)
        @test typeof(x2) === typeof(x)
        @test levels(x2) == levels(x)

        @test x[1] === Nullable(x.pool.valindex[1])
        @test x[2] === Nullable(x.pool.valindex[2])
        @test x[3] === Nullable(x.pool.valindex[3])
        @test x[4] === Nullable(x.pool.valindex[4])
        @test_throws BoundsError x[5]

        @test x[1:2] == Nullable{Float64}[0.0, 0.5]
        @test typeof(x[1:2]) === typeof(x)

        x[2] = 1
        @test x[1] === Nullable(x.pool.valindex[1])
        @test x[2] === Nullable(x.pool.valindex[3])
        @test x[3] === Nullable(x.pool.valindex[3])
        @test x[4] === Nullable(x.pool.valindex[4])
        @test levels(x) == unique(a)

        x[1:2] = -1
        @test x[1] === Nullable(x.pool.valindex[5])
        @test x[2] === Nullable(x.pool.valindex[5])
        @test x[3] === Nullable(x.pool.valindex[3])
        @test x[4] === Nullable(x.pool.valindex[4])
        @test levels(x) == vcat(unique(a), -1)

        push!(x, 2.0)
        @test length(x) == 5
        @test isequal(x, NullableArray([-1.0, -1.0, 1.0, 1.5, 2.0]))
        @test ordered(x) === isordered
        @test levels(x) == [0.0,  0.5,  1.0,  1.5, -1.0,  2.0]

        push!(x, x[1])
        @test length(x) == 6
        @test isequal(x, NullableArray([-1.0, -1.0, 1.0, 1.5, 2.0, -1.0]))
        @test ordered(x) === isordered
        @test levels(x) == [0.0,  0.5,  1.0,  1.5, -1.0,  2.0]

        append!(x, x)
        @test length(x) == 12
        @test isequal(x, NullableArray([-1.0, -1.0, 1.0, 1.5, 2.0, -1.0, -1.0, -1.0, 1.0, 1.5, 2.0, -1.0]))
        @test ordered(x) === isordered
        @test levels(x) == [0.0,  0.5,  1.0,  1.5, -1.0,  2.0]

        b = [2.5, 3.0, -3.5]
        y = NullableCategoricalVector{Float64, R}(b)
        append!(x, y)
        @test length(x) == 15
        @test isequal(x, NullableArray([-1.0, -1.0, 1.0, 1.5, 2.0, -1.0, -1.0, -1.0, 1.0, 1.5, 2.0, -1.0, 2.5, 3.0, -3.5]))
        @test ordered(x) === isordered
        @test levels(x) == [0.0, 0.5, 1.0, 1.5, -1.0, 2.0, -3.5, 2.5, 3.0]

        empty!(x)
        @test length(x) == 0
        @test ordered(x) === isordered
        @test levels(x) == [0.0, 0.5, 1.0, 1.5, -1.0, 2.0, -3.5, 2.5, 3.0]

        # Matrix with no null values
        for a in (["a" "b" "c"; "b" "a" "c"],
                  Nullable{String}["a" "b" "c"; "b" "a" "c"],
                  NullableArray(["a" "b" "c"; "b" "a" "c"]))
            na = eltype(a) <: Nullable ? a : convert(Array{Nullable{String}}, a)
            x = NullableCategoricalMatrix{String, R}(a, ordered=isordered)

            @test x == na
            @test ordered(x) === isordered
            @test levels(x) == map(get, unique(na))
            @test size(x) === (2, 3)
            @test length(x) === 6

            @test convert(NullableCategoricalArray, x) === x
            @test convert(NullableCategoricalArray{String}, x) === x
            @test convert(NullableCategoricalArray{String, 2}, x) === x
            @test convert(NullableCategoricalArray{String, 2, R}, x) === x
            @test convert(NullableCategoricalArray{String, 2, DefaultRefType}, x) == x
            @test convert(NullableCategoricalArray{String, 2, UInt8}, x) == x

            @test convert(NullableCategoricalMatrix, x) === x
            @test convert(NullableCategoricalMatrix{String}, x) === x
            @test convert(NullableCategoricalMatrix{String, R}, x) === x
            @test convert(NullableCategoricalMatrix{String, DefaultRefType}, x) == x
            @test convert(NullableCategoricalMatrix{String, UInt8}, x) == x

            for y in (NullableCategoricalArray(x, ordered=isordered),
                      NullableCategoricalArray{String}(x, ordered=isordered),
                      NullableCategoricalArray{String, 2}(x, ordered=isordered),
                      NullableCategoricalArray{String, 2, R}(x, ordered=isordered),
                      NullableCategoricalArray{String, 2, DefaultRefType}(x, ordered=isordered),
                      NullableCategoricalArray{String, 2, UInt8}(x, ordered=isordered),
                      NullableCategoricalMatrix(x, ordered=isordered),
                      NullableCategoricalMatrix{String}(x, ordered=isordered),
                      NullableCategoricalMatrix{String, R}(x, ordered=isordered),
                      NullableCategoricalMatrix{String, DefaultRefType}(x, ordered=isordered),
                      NullableCategoricalMatrix{String, UInt8}(x, ordered=isordered),
                      categorical(x, ordered=isordered),
                      categorical(x, false, ordered=isordered),
                      categorical(x, true, ordered=isordered))
                @test isa(y, NullableCategoricalMatrix{String})
                @test ordered(y) === isordered
                @test y == x
                @test y !== x
                @test y.refs !== x.refs
                @test y.pool !== x.pool
            end

            @test convert(NullableCategoricalArray, a) == x
            @test convert(NullableCategoricalArray{String}, a) == x
            @test convert(NullableCategoricalArray{String, 2, R}, a) == x
            @test convert(NullableCategoricalArray{String, 2, DefaultRefType}, a) == x
            @test convert(NullableCategoricalArray{String, 2, UInt8}, a) == x

            @test convert(NullableCategoricalMatrix, a) == x
            @test convert(NullableCategoricalMatrix{String}, a) == x
            @test convert(NullableCategoricalMatrix{String, R}, a) == x
            @test convert(NullableCategoricalMatrix{String, DefaultRefType}, a) == x
            @test convert(NullableCategoricalMatrix{String, UInt8}, a) == x

            @test NullableCategoricalArray{String}(a, ordered=isordered) == x
            @test NullableCategoricalArray{String, 2}(a, ordered=isordered) == x
            @test NullableCategoricalArray{String, 2}(a, ordered=isordered) == x
            @test NullableCategoricalArray{String, 2, R}(a, ordered=isordered) == x
            @test NullableCategoricalArray{String, 2, DefaultRefType}(a, ordered=isordered) == x
            @test NullableCategoricalArray{String, 2, UInt8}(a, ordered=isordered) == x

            @test NullableCategoricalMatrix(a, ordered=isordered) == x
            @test NullableCategoricalMatrix{String}(a, ordered=isordered) == x
            @test NullableCategoricalMatrix{String, R}(a, ordered=isordered) == x
            @test NullableCategoricalMatrix{String, DefaultRefType}(a, ordered=isordered) == x
            @test NullableCategoricalMatrix{String, UInt8}(a, ordered=isordered) == x

            for (y, R1, R2, compact) in ((a, DefaultRefType, UInt8, true),
                                         (a, DefaultRefType, DefaultRefType, false),
                                         (x, R, UInt8, true),
                                         (x, R, R, false))
            x2 = categorical(y, ordered=isordered)
            @test x2 == y
            if eltype(y) <: Nullable
                @test isa(x2, NullableCategoricalMatrix{String, R1})
            else
                @test isa(x2, CategoricalMatrix{String, R1})
            end
            @test ordered(x2) === isordered

            x2 = categorical(y, compact, ordered=isordered)
            @test x2 == y
            if eltype(y) <: Nullable
                @test isa(x2, NullableCategoricalMatrix{String, R2})
            else
                @test isa(x2, CategoricalMatrix{String, R2})
            end
            @test ordered(x2) === isordered
            end

            x2 = compact(x)
            @test x2 == x
            @test isa(x2, NullableCategoricalMatrix{String, UInt8})
            @test ordered(x2) === ordered(x)
            @test levels(x2) == levels(x)

            x2 = copy(x)
            @test x2 == x
            @test typeof(x2) === typeof(x)
            @test ordered(x2) === ordered(x)
            @test levels(x2) == levels(x)

            @test x[1] === Nullable(x.pool.valindex[1])
            @test x[2] === Nullable(x.pool.valindex[2])
            @test x[3] === Nullable(x.pool.valindex[2])
            @test x[4] === Nullable(x.pool.valindex[1])
            @test x[5] === Nullable(x.pool.valindex[3])
            @test x[6] === Nullable(x.pool.valindex[3])
            @test_throws BoundsError x[7]

            @test x[1,1] === Nullable(x.pool.valindex[1])
            @test x[2,1] === Nullable(x.pool.valindex[2])
            @test x[1,2] === Nullable(x.pool.valindex[2])
            @test x[2,2] === Nullable(x.pool.valindex[1])
            @test x[1,3] === Nullable(x.pool.valindex[3])
            @test x[2,3] === Nullable(x.pool.valindex[3])
            @test_throws BoundsError x[1,4]
            @test_throws BoundsError x[4,1]
            @test_throws BoundsError x[4,4]

            @test x[1:2,:] == x
            @test typeof(x[1:2,:]) === typeof(x)
            @test x[1:2,1] == Nullable{String}["a", "b"]
            @test typeof(x[1:2,1]) === NullableCategoricalVector{String, R}

            x[1] = "z"
            @test x[1] === Nullable(x.pool.valindex[4])
            @test x[2] === Nullable(x.pool.valindex[2])
            @test x[3] === Nullable(x.pool.valindex[2])
            @test x[4] === Nullable(x.pool.valindex[1])
            @test x[5] === Nullable(x.pool.valindex[3])
            @test x[6] === Nullable(x.pool.valindex[3])
            @test levels(x) == ["a", "b", "c", "z"]

            x[1,:] = "a"
            @test x[1] === Nullable(x.pool.valindex[1])
            @test x[2] === Nullable(x.pool.valindex[2])
            @test x[3] === Nullable(x.pool.valindex[1])
            @test x[4] === Nullable(x.pool.valindex[1])
            @test x[5] === Nullable(x.pool.valindex[1])
            @test x[6] === Nullable(x.pool.valindex[3])
            @test levels(x) == ["a", "b", "c", "z"]

            x[1,1:2] = "z"
            @test x[1] === Nullable(x.pool.valindex[4])
            @test x[2] === Nullable(x.pool.valindex[2])
            @test x[3] === Nullable(x.pool.valindex[4])
            @test x[4] === Nullable(x.pool.valindex[1])
            @test x[5] === Nullable(x.pool.valindex[1])
            @test x[6] === Nullable(x.pool.valindex[3])
            @test levels(x) == ["a", "b", "c", "z"]
        end

        # Matrix with null values
        for a in (Nullable{String}["a" Nullable() "c"; "b" "a" Nullable()],
                  NullableArray(Nullable{String}["a" Nullable() "c"; "b" "a" Nullable()]))
            x = NullableCategoricalMatrix{String, R}(a, ordered=isordered)

            @test x == a
            @test ordered(x) === isordered
            @test levels(x) == map(get, filter(x->!isnull(x), unique(a)))
            @test size(x) === (2, 3)
            @test length(x) === 6

            @test convert(NullableCategoricalArray, x) === x
            @test convert(NullableCategoricalArray{String}, x) === x
            @test convert(NullableCategoricalArray{String, 2}, x) === x
            @test convert(NullableCategoricalArray{String, 2, R}, x) === x
            @test convert(NullableCategoricalArray{String, 2, DefaultRefType}, x) == x
            @test convert(NullableCategoricalArray{String, 2, UInt8}, x) == x

            @test convert(NullableCategoricalMatrix, x) === x
            @test convert(NullableCategoricalMatrix{String}, x) === x
            @test convert(NullableCategoricalMatrix{String, R}, x) === x
            @test convert(NullableCategoricalMatrix{String, DefaultRefType}, x) == x
            @test convert(NullableCategoricalMatrix{String, UInt8}, x) == x

            for y in (NullableCategoricalArray(x, ordered=isordered),
                      NullableCategoricalArray{String}(x, ordered=isordered),
                      NullableCategoricalArray{String, 2}(x, ordered=isordered),
                      NullableCategoricalArray{String, 2, R}(x, ordered=isordered),
                      NullableCategoricalArray{String, 2, DefaultRefType}(x, ordered=isordered),
                      NullableCategoricalArray{String, 2, UInt8}(x, ordered=isordered),
                      NullableCategoricalMatrix(x, ordered=isordered),
                      NullableCategoricalMatrix{String}(x, ordered=isordered),
                      NullableCategoricalMatrix{String, R}(x, ordered=isordered),
                      NullableCategoricalMatrix{String, DefaultRefType}(x, ordered=isordered),
                      NullableCategoricalMatrix{String, UInt8}(x, ordered=isordered),
                      categorical(x, ordered=isordered),
                      categorical(x, false, ordered=isordered),
                      categorical(x, true, ordered=isordered))
                @test isa(y, NullableCategoricalMatrix{String})
                @test ordered(y) === isordered
                @test y == x
                @test y !== x
                @test y.refs !== x.refs
                @test y.pool !== x.pool
            end

            @test convert(NullableCategoricalArray, a) == x
            @test convert(NullableCategoricalArray{String}, a) == x
            @test convert(NullableCategoricalArray{String, 2, R}, a) == x
            @test convert(NullableCategoricalArray{String, 2, DefaultRefType}, a) == x
            @test convert(NullableCategoricalArray{String, 2, UInt8}, a) == x

            @test convert(NullableCategoricalMatrix, a) == x
            @test convert(NullableCategoricalMatrix{String}, a) == x
            @test convert(NullableCategoricalMatrix{String, R}, a) == x
            @test convert(NullableCategoricalMatrix{String, DefaultRefType}, a) == x
            @test convert(NullableCategoricalMatrix{String, UInt8}, a) == x

            @test NullableCategoricalArray{String}(a, ordered=isordered) == x
            @test NullableCategoricalArray{String, 2}(a, ordered=isordered) == x
            @test NullableCategoricalArray{String, 2}(a, ordered=isordered) == x
            @test NullableCategoricalArray{String, 2, R}(a, ordered=isordered) == x
            @test NullableCategoricalArray{String, 2, DefaultRefType}(a, ordered=isordered) == x
            @test NullableCategoricalArray{String, 2, UInt8}(a, ordered=isordered) == x

            @test NullableCategoricalMatrix(a, ordered=isordered) == x
            @test NullableCategoricalMatrix{String}(a, ordered=isordered) == x
            @test NullableCategoricalMatrix{String, R}(a, ordered=isordered) == x
            @test NullableCategoricalMatrix{String, DefaultRefType}(a, ordered=isordered) == x
            @test NullableCategoricalMatrix{String, UInt8}(a, ordered=isordered) == x

            for (y, R1, R2, compact) in ((a, DefaultRefType, UInt8, true),
                                         (a, DefaultRefType, DefaultRefType, false),
                                         (x, R, UInt8, true),
                                         (x, R, R, false))
                x2 = categorical(y, ordered=isordered)
                @test x2 == y
                @test isa(x2, NullableCategoricalMatrix{String, R1})
                @test ordered(x2) === isordered

                x2 = categorical(y, compact, ordered=isordered)
                @test x2 == y
                @test isa(x2, NullableCategoricalMatrix{String, R2})
                @test ordered(x2) === isordered
            end

            x2 = compact(x)
            @test x2 == x
            @test isa(x2, NullableCategoricalMatrix{String, UInt8})
            @test ordered(x2) === ordered(x)
            @test levels(x2) == levels(x)

            x2 = copy(x)
            @test x2 == x
            @test typeof(x2) === typeof(x)
            @test ordered(x2) === ordered(x)
            @test levels(x2) == levels(x)

            @test x[1] === Nullable(x.pool.valindex[1])
            @test x[2] === Nullable(x.pool.valindex[2])
            @test x[3] === eltype(x)()
            @test x[4] === Nullable(x.pool.valindex[1])
            @test x[5] === Nullable(x.pool.valindex[3])
            @test x[6] === eltype(x)()
            @test_throws BoundsError x[7]

            @test x[1,1] === Nullable(x.pool.valindex[1])
            @test x[2,1] === Nullable(x.pool.valindex[2])
            @test x[1,2] === eltype(x)()
            @test x[2,2] === Nullable(x.pool.valindex[1])
            @test x[1,3] === Nullable(x.pool.valindex[3])
            @test x[2,3] === eltype(x)()
            @test_throws BoundsError x[1,4]
            @test_throws BoundsError x[4,1]
            @test_throws BoundsError x[4,4]

            @test x[1:2,:] == x
            @test typeof(x[1:2,:]) === typeof(x)
            @test x[1:2,1] == Nullable{String}["a", "b"]
            @test typeof(x[1:2,1]) === NullableCategoricalVector{String, R}

            x[1] = "z"
            @test x[1] === Nullable(x.pool.valindex[4])
            @test x[2] === Nullable(x.pool.valindex[2])
            @test x[3] === eltype(x)()
            @test x[4] === Nullable(x.pool.valindex[1])
            @test x[5] === Nullable(x.pool.valindex[3])
            @test x[6] === eltype(x)()
            @test levels(x) == ["a", "b", "c", "z"]

            x[1,:] = "a"
            @test x[1] === Nullable(x.pool.valindex[1])
            @test x[2] === Nullable(x.pool.valindex[2])
            @test x[3] === Nullable(x.pool.valindex[1])
            @test x[4] === Nullable(x.pool.valindex[1])
            @test x[5] === Nullable(x.pool.valindex[1])
            @test x[6] === eltype(x)()
            @test levels(x) == ["a", "b", "c", "z"]

            x[1,1:2] = Nullable("z")
            @test x[1] === Nullable(x.pool.valindex[4])
            @test x[2] === Nullable(x.pool.valindex[2])
            @test x[3] === Nullable(x.pool.valindex[4])
            @test x[4] === Nullable(x.pool.valindex[1])
            @test x[5] === Nullable(x.pool.valindex[1])
            @test x[6] === eltype(x)()
            @test levels(x) == ["a", "b", "c", "z"]

            x[1] = Nullable()
            @test x[1] === eltype(x)()
            @test x[2] === Nullable(x.pool.valindex[2])
            @test x[3] === Nullable(x.pool.valindex[4])
            @test x[4] === Nullable(x.pool.valindex[1])
            @test x[5] === Nullable(x.pool.valindex[1])
            @test x[6] === eltype(x)()
            @test levels(x) == ["a", "b", "c", "z"]

            x[1,1:2] = Nullable()
            @test x[1] === eltype(x)()
            @test x[2] === Nullable(x.pool.valindex[2])
            @test x[3] === eltype(x)()
            @test x[4] === Nullable(x.pool.valindex[1])
            @test x[5] === Nullable(x.pool.valindex[1])
            @test x[6] === eltype(x)()
            @test levels(x) == ["a", "b", "c", "z"]

            x[:,2] = Nullable()
            @test x[1] === eltype(x)()
            @test x[2] === Nullable(x.pool.valindex[2])
            @test x[3] === eltype(x)()
            @test x[4] === eltype(x)()
            @test x[5] === Nullable(x.pool.valindex[1])
            @test x[6] === eltype(x)()
            @test levels(x) == ["a", "b", "c", "z"]

            # Constructor with values plus missingness array
            x = NullableCategoricalArray(1:3, [true, false, true], ordered=isordered)
            @test x == Nullable{Int}[Nullable(), 2, Nullable()]
            @test ordered(x) == isordered
            @test levels(x) == [2]

            if VERSION >= v"0.5.0-dev"
                x = NullableCategoricalVector(1:3, [true, false, true], ordered=isordered)
                @test x == Nullable{Int}[Nullable(), 2, Nullable()]
                @test ordered(x) === isordered
                @test levels(x) == [2]

                x = NullableCategoricalMatrix([1 2; 3 4], [true false; false true],
                                              ordered=isordered)
                @test x == Nullable{Int}[Nullable() 2; 3 Nullable()]
                @test ordered(x) === isordered
                @test levels(x) == [2, 3]
            end
        end

        # Uninitialized array
        v = Any[NullableCategoricalArray(2, ordered=isordered),
                NullableCategoricalArray(String, 2, ordered=isordered),
                NullableCategoricalArray{String}(2, ordered=isordered),
                NullableCategoricalArray{String, 1}(2, ordered=isordered),
                NullableCategoricalArray{String, 1, R}(2, ordered=isordered),
                NullableCategoricalVector{String}(2, ordered=isordered),
                NullableCategoricalVector{String, R}(2, ordered=isordered),
                NullableCategoricalArray(2, 3, ordered=isordered),
                NullableCategoricalArray(String, 2, 3, ordered=isordered),
                NullableCategoricalArray{String}(2, 3, ordered=isordered),
                NullableCategoricalArray{String, 2}(2, 3, ordered=isordered),
                NullableCategoricalArray{String, 2, R}(2, 3, ordered=isordered),
                NullableCategoricalMatrix{String}(2, 3, ordered=isordered),
                NullableCategoricalMatrix{String, R}(2, 3, ordered=isordered)]

        # See conditional definition of constructors in array.jl and nullablearray.jl
        if VERSION >= v"0.5.0-dev"
            push!(v, NullableCategoricalVector(2, ordered=isordered),
                     NullableCategoricalVector(String, 2, ordered=isordered),
                     NullableCategoricalMatrix(2, 3, ordered=isordered),
                     NullableCategoricalMatrix(String, 2, 3, ordered=isordered))
        end

        for x in v
            @test ordered(x) === isordered
            @test isnull(x[1])
            @test isnull(x[2])
            @test levels(x) == []

            x2 = compact(x)
            @test x2 == x
            if VERSION >= v"0.5.0-dev"
                @test isa(x2, NullableCategoricalArray{String, ndims(x), UInt8})
            else
                @test isa(x2, NullableCategoricalArray{typeof(x).parameters[1], ndims(x), UInt8})
            end
            @test ordered(x2) === ordered(x)
            @test levels(x2) == []

            x2 = copy(x)
            @test x2 == x
            @test typeof(x2) === typeof(x)
            @test ordered(x2) === ordered(x)
            @test levels(x2) == []

            x[1] = "c"
            @test x[1] === Nullable(x.pool.valindex[1])
            @test isnull(x[2])
            @test levels(x) == ["c"]

            x[1] = "a"
            @test x[1] === Nullable(x.pool.valindex[2])
            @test isnull(x[2])
            @test levels(x) == ["c", "a"]

            x[2] = Nullable()
            @test x[1] === Nullable(x.pool.valindex[2])
            @test x[2] === eltype(x)()
            @test levels(x) == ["c", "a"]

            x[1] = Nullable("b")
            @test x[1] === Nullable(x.pool.valindex[3])
            @test x[2] === eltype(x)()
            @test levels(x) == ["c", "a", "b"]
        end
    end
end

end


text:
using Tetris
using Base.Test

@test isdefined(Main, Symbol("@tetris"))
@test isdefined(Tetris, :setup)