text:
#
# TimestepVector and TimestepMatrix
#

#
# a.  General
#

# Get a timestep array of type T with N dimensions. Time labels will match those from the time dimension in md
function get_timestep_array(md::ModelDef, T, N, ti, value)
	if isuniform(md)
		first, stepsize = first_and_step(md)
		first === nothing && @warn "get_timestep_array: first === nothing"
        return TimestepArray{FixedTimestep{first, stepsize}, T, N, ti}(value)
    else
        TIMES = (time_labels(md)...,)
        return TimestepArray{VariableTimestep{TIMES}, T, N, ti}(value)
    end
end

# Return the index position of the time dimension in the datumdef or parameter. If there is no time dimension, return nothing
get_time_index_position(dims::Union{Nothing, Array{Symbol}}) = findfirst(isequal(:time), dims)
get_time_index_position(obj::Union{AbstractDatumDef, ArrayModelParameter}) = get_time_index_position(dim_names(obj))

function get_time_index_position(obj::AbstractCompositeComponentDef, comp_name::Symbol, datum_name::Symbol)
	get_time_index_position(dim_names(compdef(obj, comp_name), datum_name))
end

const AnyIndex = Union{Int, Vector{Int}, Tuple, Colon, OrdinalRange}
const AnyIndex_NonColon = Union{Int, Vector{Int}, Tuple, OrdinalRange}

# Helper function for getindex; throws a MissingException if data is missing, otherwise returns data
function _missing_data_check(data, t)
	if data === missing
		throw(MissingException("Cannot get index; data is missing. You may have tried to access a value in timestep $t that has not yet been computed."))
	else
		return data
	end
end

# Helper function for getindex; throws an error if the TimestepIndex index is out of range of the TimestepArray
function _index_bounds_check(data, dim, t)
	if size(data, dim) < t
		error("TimestepIndex index $t extends beyond bounds of TimestepArray dimension $dim")
	end
end

# Helper function for getindex; throws an error if you index into a N-dimensional TimestepArray with only one index
# if N > 1; Note that Julia does allow this and returns the column-major value, but this could produce unexpected
# behavior for users in this case so we do not allow it for now
function _single_index_check(data, idxs)
	num_idxs = length(idxs)
	num_dims = length(size(data))
	if num_idxs < num_dims
		error("Not enough indices provided to index into TimestepArray, $num_idxs provided, $num_dims required")
	end
end

# Helper functionfor getindex; throws an error if one indexes into a TimestepArray with an integer
function _throw_int_getindex_depwarning()
	msg = "Indexing with getindex into a TimestepArray with Integer(s) is deprecated, please index with a TimestepIndex(index::Int) instead ie. instead of t[2] use t[TimestepIndex(2)]"
	Base.depwarn("$msg, $(stacktrace())", :getindex)
end

# Helper function for setindex; throws an deprecation warning if one indexes into a TimestepArray with an integer
function _throw_int_setindex_depwarning()
	msg = "Indexing with setindex into a TimestepArray with Integer(s) is deprecated, please index with a TimestepIndex(index::Int) instead ie. instead of t[2] use t[TimestepIndec(2)]"
	Base.depwarn("$msg, $(stacktrace())", :setindex)
end

# Helper macro used by connector
macro allow_missing(expr)
	let e = gensym("e")
		retexpr = quote
			try
				$expr
			catch $e
				if $e isa MissingException
					missing
				else
					rethrow($e)
				end
			end
		end
		return esc(retexpr)
	end
end

# Helper functions for TimestepValue type
function _get_time_value_position(times::Union{Tuple, Array}, ts::TimestepValue{T}) where T
	t = findfirst(isequal.(ts.value, times))
	if t === nothing
		error("cannot use TimestepValue with value $(ts.value), value is not in the TimestepArray")
	end

	t_offset = t + ts.offset
	if t_offset > length(times)
		error("cannot get TimestepValue offset of $(ts.offset) from value $(ts.value), offset is after the end of the TimestepArray")
	end
	return t_offset
end

# Helper function to get the array of indices from an Array{TimestepIndex,1}
function _get_ts_indices(ts_array::Array{TimestepIndex, 1})
    return [ts.index for ts in ts_array]
end

function _get_ts_indices(ts_array::Array{TimestepValue{T}, 1}, times::Union{Tuple, Array}) where T
	return [_get_time_value_position(times, ts) for ts in ts_array]
end

#
# b. TimestepVector
#

function Base.getindex(v::TimestepVector{FixedTimestep{FIRST, STEP}, T}, ts::FixedTimestep{FIRST, STEP, LAST}) where {T, FIRST, STEP, LAST}
	data = v.data[ts.t]
	_missing_data_check(data, ts.t)
end

function Base.getindex(v::TimestepVector{VariableTimestep{TIMES}, T}, ts::VariableTimestep{TIMES}) where {T, TIMES}
	data = v.data[ts.t]
	_missing_data_check(data, ts.t)
end

function Base.getindex(v::TimestepVector{FixedTimestep{D_FIRST, STEP}, T}, ts::FixedTimestep{T_FIRST, STEP, LAST}) where {T, D_FIRST, T_FIRST, STEP, LAST}
	t = Int(ts.t + (T_FIRST - D_FIRST) / STEP)
	data = v.data[t]
	_missing_data_check(data, t)
end

function Base.getindex(v::TimestepVector{VariableTimestep{D_TIMES}, T}, ts::VariableTimestep{T_TIMES}) where {T, D_TIMES, T_TIMES}
	t = ts.t + findfirst(isequal(T_TIMES[1]), D_TIMES) - 1
	data = v.data[t]
	_missing_data_check(data, t)
end

function Base.getindex(v::TimestepVector{FixedTimestep{FIRST, STEP}, T_data}, ts::TimestepValue{T_time}) where {T_data, FIRST, STEP, T_time}
	LAST = FIRST + ((length(v.data)-1) * STEP)
	t = _get_time_value_position([FIRST:STEP:LAST...], ts)
	data = v.data[t]
	_missing_data_check(data, t)
end

function Base.getindex(v::TimestepVector{VariableTimestep{TIMES}, T_data}, ts::TimestepValue{T_time}) where {T_data, TIMES, T_time}
	t = _get_time_value_position(TIMES, ts)
	data = v.data[t]
	_missing_data_check(data, t)
end

function Base.getindex(v::TimestepVector, ts::TimestepIndex)
	t = ts.index
	_index_bounds_check(v.data, 1, t)
	data = v.data[t]
	_missing_data_check(data, t)
end

function Base.setindex!(v::TimestepVector{FixedTimestep{FIRST, STEP}, T}, val, ts::FixedTimestep{FIRST, STEP, LAST}) where {T, FIRST, STEP, LAST}
	setindex!(v.data, val, ts.t)
end

function Base.setindex!(v::TimestepVector{VariableTimestep{TIMES}, T}, val, ts::VariableTimestep{TIMES}) where {T, TIMES}
	setindex!(v.data, val, ts.t)
end

function Base.setindex!(v::TimestepVector{FixedTimestep{D_FIRST, STEP}, T}, val, ts::FixedTimestep{T_FIRST, STEP, LAST}) where {T, D_FIRST, T_FIRST, STEP, LAST}
	t = Int(ts.t + (T_FIRST - D_FIRST) / STEP)
	setindex!(v.data, val, t)
end

function Base.setindex!(v::TimestepVector{VariableTimestep{D_TIMES}, T}, val, ts::VariableTimestep{T_TIMES}) where {T, D_TIMES, T_TIMES}
	t = ts.t + findfirst(isequal(T_TIMES[1]), D_TIMES) - 1
	setindex!(v.data, val, t)
end

function Base.setindex!(v::TimestepVector{FixedTimestep{FIRST, STEP}, T_data}, val, ts::TimestepValue{T_time}) where {T_data, FIRST, STEP, T_time}
	LAST = FIRST + ((length(v.data)-1) * STEP)
	t = _get_time_value_position([FIRST:STEP:LAST...], ts)
	setindex!(v.data, val, t)
end

function Base.setindex!(v::TimestepVector{VariableTimestep{TIMES}, T_data}, val, ts::TimestepValue{T_time}) where {T_data, TIMES, T_time}
	t = _get_time_value_position(TIMES, ts)
	setindex!(v.data, val, t)
end

function Base.setindex!(v::TimestepVector, val, ts::TimestepIndex)
	setindex!(v.data, val, ts.index)
end

# int indexing version supports old-style components and internal functions, not
# part of the public API

function Base.getindex(v::TimestepVector{FixedTimestep{FIRST, STEP}, T}, i::AnyIndex_NonColon) where {T, FIRST, STEP}
	_throw_int_getindex_depwarning()
	return v.data[i]
end

function Base.getindex(v::TimestepVector{VariableTimestep{TIMES}, T}, i::AnyIndex_NonColon) where {T, TIMES}
	_throw_int_getindex_depwarning()
	return v.data[i]
end

function Base.setindex!(v::TimestepVector{FixedTimestep{Start, STEP}, T}, val, i::AnyIndex_NonColon) where {T, Start, STEP}
	_throw_int_setindex_depwarning()
	setindex!(v.data, val, i)
end

function Base.setindex!(v::TimestepVector{VariableTimestep{TIMES}, T}, val, i::AnyIndex_NonColon) where {T, TIMES}
	_throw_int_setindex_depwarning()
	setindex!(v.data, val, i)
end

function Base.length(v::TimestepVector)
	return length(v.data)
end

Base.lastindex(v::TimestepVector) = length(v)

#
# c. TimestepMatrix
#

function Base.getindex(mat::TimestepMatrix{FixedTimestep{FIRST, STEP}, T, 1}, ts::FixedTimestep{FIRST, STEP, LAST}, idx::AnyIndex) where {T, FIRST, STEP, LAST}
	data = mat.data[ts.t, idx]
	_missing_data_check(data, ts.t)
end

function Base.getindex(mat::TimestepMatrix{VariableTimestep{TIMES}, T, 1}, ts::VariableTimestep{TIMES}, idx::AnyIndex) where {T, TIMES}
	data = mat.data[ts.t, idx]
	_missing_data_check(data, ts.t)
end

function Base.getindex(mat::TimestepMatrix{FixedTimestep{D_FIRST, STEP}, T, 1}, ts::FixedTimestep{T_FIRST, STEP, LAST}, idx::AnyIndex) where {T, D_FIRST, T_FIRST, STEP, LAST}
	t = Int(ts.t + (T_FIRST - D_FIRST) / STEP)
	data = mat.data[t, idx]
	_missing_data_check(data, t)
end

function Base.getindex(mat::TimestepMatrix{VariableTimestep{D_TIMES}, T, 1}, ts::VariableTimestep{T_TIMES}, idx::AnyIndex) where {T, D_TIMES, T_TIMES}
	t = ts.t + findfirst(isequal(T_TIMES[1]), D_TIMES) - 1
	data = mat.data[t, idx]
	_missing_data_check(data, t)
end

function Base.getindex(mat::TimestepMatrix{FixedTimestep{FIRST, STEP}, T, 2}, idx::AnyIndex, ts::FixedTimestep{FIRST, STEP, LAST}) where {T, FIRST, STEP, LAST}
	data = mat.data[idx, ts.t]
	_missing_data_check(data, ts.t)
end

function Base.getindex(mat::TimestepMatrix{VariableTimestep{TIMES}, T, 2}, idx::AnyIndex, ts::VariableTimestep{TIMES}) where {T, TIMES}
	# WAS THIS: data = mat.data[ts.t, idx, ts.t]
	data = mat.data[idx, ts.t]
	_missing_data_check(data, ts.t)
end

function Base.getindex(mat::TimestepMatrix{FixedTimestep{D_FIRST, STEP}, T, 2}, idx::AnyIndex, ts::FixedTimestep{T_FIRST, STEP, LAST}) where {T, D_FIRST, T_FIRST, STEP, LAST}
	t = Int(ts.t + (T_FIRST - D_FIRST) / STEP)
	data = mat.data[idx, t]
	_missing_data_check(data, t)
end

function Base.getindex(mat::TimestepMatrix{VariableTimestep{D_TIMES}, T, 2}, idx::AnyIndex, ts::VariableTimestep{T_TIMES}) where {T, D_TIMES, T_TIMES}
	t = ts.t + findfirst(isequal(T_TIMES[1]), D_TIMES) - 1
	data = mat.data[idx, t]
	_missing_data_check(data, t)
end

function Base.getindex(mat::TimestepMatrix{FixedTimestep{FIRST, STEP}, T_data, 1}, ts::TimestepValue{T_time}, idx::AnyIndex) where {T_data, FIRST, STEP, T_time}
	LAST = FIRST + ((size(mat.data, 1) - 1) * STEP)
	t = _get_time_value_position([FIRST:STEP:LAST...], ts)
	data = mat.data[t, idx]
	_missing_data_check(data, t)
end

function Base.getindex(mat::TimestepMatrix{VariableTimestep{TIMES}, T_data, 1}, ts::TimestepValue{T_time}, idx::AnyIndex) where {T_data, TIMES, T_time}
	t = _get_time_value_position(TIMES, ts)
	data = mat.data[t, idx]
	_missing_data_check(data, t)
end

function Base.getindex(mat::TimestepMatrix, ts::TimestepIndex, idx::AnyIndex)
	t = ts.index
	_index_bounds_check(mat.data, 1, t)
	data = mat.data[t, idx]
	_missing_data_check(data, t)
end

function Base.getindex(mat::TimestepMatrix{FixedTimestep{FIRST, STEP}, T_data, 2}, idx::AnyIndex, ts::TimestepValue{T_time}) where {T_data, FIRST, STEP, T_time}
	LAST = FIRST + ((size(mat.data, 2) - 1) * STEP)
	t = _get_time_value_position([FIRST:STEP:LAST...], ts)
	data = mat.data[idx, t]
	_missing_data_check(data, t)
end

function Base.getindex(mat::TimestepMatrix{VariableTimestep{TIMES}, T_data, 2}, idx::AnyIndex, ts::TimestepValue{T_time}) where {T_data, TIMES, T_time}
	t = _get_time_value_position(TIMES, ts)
	data = mat.data[idx, t]
	_missing_data_check(data, t)
end

function Base.getindex(mat::TimestepMatrix, idx::AnyIndex, ts::TimestepIndex)
	t = ts.index
	_index_bounds_check(mat.data, 2, t)
	data = mat.data[idx, t]
	_missing_data_check(data, t)
end

function Base.setindex!(mat::TimestepMatrix{FixedTimestep{FIRST, STEP}, T, 1}, val, ts::FixedTimestep{FIRST, STEP, LAST}, idx::AnyIndex) where {T, FIRST, STEP, LAST}
	setindex!(mat.data, val, ts.t, idx)
end

function Base.setindex!(mat::TimestepMatrix{VariableTimestep{TIMES}, T, 1}, val, ts::VariableTimestep{TIMES}, idx::AnyIndex) where {T, TIMES}
	setindex!(mat.data, val, ts.t, idx)
end

function Base.setindex!(mat::TimestepMatrix{FixedTimestep{D_FIRST, STEP}, T, 1}, val, ts::FixedTimestep{T_FIRST, STEP, LAST}, idx::AnyIndex) where {T, D_FIRST, T_FIRST, STEP, LAST}
	t = Int(ts.t + (T_FIRST - D_FIRST) / STEP)
	setindex!(mat.data, val, t, idx)
end

function Base.setindex!(mat::TimestepMatrix{VariableTimestep{D_TIMES}, T, 1}, val, ts::VariableTimestep{T_TIMES}, idx::AnyIndex) where {T, D_TIMES, T_TIMES}
	t = ts.t + findfirst(isequal(T_TIMES[1]), D_TIMES) - 1
	setindex!(mat.data, val, t, idx)
end

function Base.setindex!(mat::TimestepMatrix{FixedTimestep{FIRST, STEP}, T, 2}, val, idx::AnyIndex, ts::FixedTimestep{FIRST, STEP, LAST}) where {T, FIRST, STEP, LAST}
	setindex!(mat.data, val, idx, ts.t)
end

function Base.setindex!(mat::TimestepMatrix{VariableTimestep{TIMES}, T, 2}, val, idx::AnyIndex, ts::VariableTimestep{TIMES}) where {T, TIMES}
	setindex!(mat.data, val, idx, ts.t)
end

function Base.setindex!(mat::TimestepMatrix{FixedTimestep{D_FIRST, STEP}, T, 2}, val, idx::AnyIndex, ts::FixedTimestep{T_FIRST, STEP, LAST}) where {T, D_FIRST, T_FIRST, STEP, LAST}
	t = Int(ts.t + (T_FIRST - D_FIRST) / STEP)
	setindex!(mat.data, val, idx, t)
end

function Base.setindex!(mat::TimestepMatrix{VariableTimestep{D_TIMES}, T, 2}, val, idx::AnyIndex, ts::VariableTimestep{T_TIMES}) where {T, D_TIMES, T_TIMES}
	t = ts.t + findfirst(isequal(T_TIMES[1]), D_TIMES) - 1
	setindex!(mat.data, val, idx, t)
end

function Base.setindex!(mat::TimestepMatrix{FixedTimestep{FIRST, STEP}, T_data, 1}, val, ts::TimestepValue{T_time}, idx::AnyIndex) where {T_data, FIRST, STEP, T_time}
	LAST = FIRST + ((size(mat.data, 1) - 1) * STEP)
	t = _get_time_value_position([FIRST:STEP:LAST...], ts)
	setindex!(mat.data, val, t, idx)
end

function Base.setindex!(mat::TimestepMatrix{VariableTimestep{TIMES}, T_data, 1}, val, ts::TimestepValue{T_time}, idx::AnyIndex) where {T_data, TIMES, T_time}
	t = _get_time_value_position(TIMES, ts)
	setindex!(mat.data, val, t, idx)
end

function Base.setindex!(mat::TimestepMatrix{FixedTimestep{FIRST, STEP}, T_data, 2}, val, idx::AnyIndex, ts::TimestepValue{T_time}) where {T_data, FIRST, STEP, T_time}
	LAST = FIRST + ((size(mat.data, 1) - 1) * STEP)
	t = _get_time_value_position([FIRST:STEP:LAST...], ts)
	setindex!(mat.data, val, idx, t)
end

function Base.setindex!(mat::TimestepMatrix{VariableTimestep{TIMES}, T_data, 2}, val, idx::AnyIndex, ts::TimestepValue{T_time}) where {T_data, TIMES, T_time}
	t = _get_time_value_position(TIMES, ts)
	setindex!(mat.data, val, idx, t)
end

function Base.setindex!(mat::TimestepMatrix, val, idx::AnyIndex, ts::TimestepIndex)
	setindex!(mat.data, val, idx, ts.index)
end

function Base.setindex!(mat::TimestepMatrix, val, ts::TimestepIndex, idx::AnyIndex)
	setindex!(mat.data, val, ts.index, idx)
end

# int indexing version supports old-style components and internal functions, not
# part of the public API

function Base.getindex(mat::TimestepMatrix{FixedTimestep{FIRST, STEP}, T, ti}, idx1::AnyIndex_NonColon, idx2::AnyIndex_NonColon) where {T, FIRST, STEP, ti}
	_throw_int_getindex_depwarning()
	return mat.data[idx1, idx2]
end

function Base.getindex(mat::TimestepMatrix{VariableTimestep{TIMES}, T, ti}, idx1::AnyIndex_NonColon, idx2::AnyIndex_NonColon) where {T, TIMES, ti}
	_throw_int_getindex_depwarning()
	return mat.data[idx1, idx2]
end

function Base.setindex!(mat::TimestepMatrix{FixedTimestep{FIRST, STEP}, T, ti}, val, idx1::Int, idx2::Int) where {T, FIRST, STEP, ti}
	_throw_int_setindex_depwarning()
	setindex!(mat.data, val, idx1, idx2)
end

function Base.setindex!(mat::TimestepMatrix{FixedTimestep{FIRST, STEP}, T, ti}, val, idx1::AnyIndex_NonColon, idx2::AnyIndex_NonColon) where {T, FIRST, STEP, ti}
	_throw_int_setindex_depwarning()
	mat.data[idx1, idx2] .= val
end

function Base.setindex!(mat::TimestepMatrix{VariableTimestep{TIMES}, T, ti}, val, idx1::Int, idx2::Int) where {T, TIMES, ti}
	_throw_int_setindex_depwarning()
	setindex!(mat.data, val, idx1, idx2)
end

function Base.setindex!(mat::TimestepMatrix{VariableTimestep{TIMES}, T, ti}, val, idx1::AnyIndex_NonColon, idx2::AnyIndex_NonColon) where {T, TIMES, ti}
	_throw_int_setindex_depwarning()
	mat.data[idx1, idx2] .= val
end

#
# TimestepArray methods
#
function Base.dotview(v::Mimi.TimestepArray, args...)
	# convert any timesteps to their underlying index
	args = map(arg -> (arg isa AbstractTimestep ? arg.t : arg), args)
	Base.dotview(v.data, args...)
end

Base.fill!(obj::TimestepArray, value) = fill!(obj.data, value)

Base.size(obj::TimestepArray) = size(obj.data)

Base.size(obj::TimestepArray, i::Int) = size(obj.data, i)

Base.ndims(obj::TimestepArray{T_ts, T, N, ti}) where {T_ts, T, N, ti} = N

Base.eltype(obj::TimestepArray{T_ts, T, N, ti}) where {T_ts, T, N, ti} = T

first_period(obj::TimestepArray{FixedTimestep{FIRST,STEP}, T, N, ti}) where {FIRST, STEP, T, N, ti} = FIRST
first_period(obj::TimestepArray{VariableTimestep{TIMES}, T, N, ti}) where {TIMES, T, N, ti} = TIMES[1]

last_period(obj::TimestepArray{FixedTimestep{FIRST, STEP}, T, N, ti}) where {FIRST, STEP, T, N, ti} = (FIRST + (size(obj, 1) - 1) * STEP)
last_period(obj::TimestepArray{VariableTimestep{TIMES}, T, N, ti}) where {TIMES, T, N, ti} = TIMES[end]

time_labels(obj::TimestepArray{FixedTimestep{FIRST, STEP}, T, N, ti}) where {FIRST, STEP, T, N, ti} = collect(FIRST:STEP:(FIRST + (size(obj, 1) - 1) * STEP))
time_labels(obj::TimestepArray{VariableTimestep{TIMES}, T, N, ti}) where {TIMES, T, N, ti} = collect(TIMES)

split_indices(idxs, ti) = idxs[1:ti - 1], idxs[ti], idxs[ti + 1:end]

function Base.getindex(arr::TimestepArray{FixedTimestep{FIRST, STEP}, T, N, ti}, idxs::Union{FixedTimestep{FIRST, STEP, LAST}, AnyIndex}...) where {T, N, ti, FIRST, STEP, LAST}
	idxs1, ts, idxs2 = split_indices(idxs, ti)
	return arr.data[idxs1..., ts.t, idxs2...]
end

function Base.getindex(arr::TimestepArray{VariableTimestep{TIMES}, T, N, ti}, idxs::Union{VariableTimestep{TIMES}, AnyIndex}...) where {T, N, ti, TIMES}
	idxs1, ts, idxs2 = split_indices(idxs, ti)
	return arr.data[idxs1..., ts.t, idxs2...]
end

function Base.getindex(arr::TimestepArray{FixedTimestep{D_FIRST, STEP}, T, N, ti}, idxs::Union{FixedTimestep{T_FIRST, STEP, LAST}, AnyIndex}...) where {T, N, ti, D_FIRST, T_FIRST, STEP, LAST}
	idxs1, ts, idxs2 = split_indices(idxs, ti)
	t = Int(ts.t + (FIRST - TIMES[1]) / STEP)
	return arr.data[idxs1..., t, idxs2...]
end

function Base.getindex(arr::TimestepArray{VariableTimestep{D_TIMES}, T, N, ti}, idxs::Union{VariableTimestep{T_TIMES}, AnyIndex}...) where {T, N, ti, D_TIMES, T_TIMES}
	idxs1, ts, idxs2 = split_indices(idxs, ti)
	t = ts.t + findfirst(isequal(T_TIMES[1]), D_TIMES) - 1
	return arr.data[idxs1..., t, idxs2...]
end

function Base.getindex(arr::TimestepArray{FixedTimestep{FIRST, STEP}, T_data, N, ti}, idxs::Union{TimestepValue{T_time}, AnyIndex}...) where {T_data, N, ti, FIRST, STEP, T_time}
	_single_index_check(arr.data, idxs)
	idxs1, ts, idxs2 = split_indices(idxs, ti)
	LAST = FIRST + ((size(arr.data, ti) - 1) * STEP)
	t = _get_time_value_position([FIRST:STEP:LAST...], ts)
	return arr.data[idxs1..., t, idxs2...]
end

function Base.getindex(arr::TimestepArray{VariableTimestep{TIMES}, T_data, N, ti}, idxs::Union{TimestepValue{T_time}, AnyIndex}...) where {T_data, N, ti, TIMES, T_time}
	_single_index_check(arr.data, idxs)
	idxs1, ts, idxs2 = split_indices(idxs, ti)
	t = _get_time_value_position(TIMES, ts)
	return arr.data[idxs1..., t, idxs2...]
end

function Base.getindex(arr::TimestepArray{FixedTimestep{FIRST, STEP}, T, N, ti}, idxs::Union{TimestepIndex, AnyIndex}...) where {T, N, ti, FIRST, STEP}
	_single_index_check(arr.data, idxs)
	idxs1, ts, idxs2 = split_indices(idxs, ti)
	t = ts.index
	_index_bounds_check(arr.data, ti, t)
	return arr.data[idxs1..., t, idxs2...]
end

function Base.getindex(arr::TimestepArray{VariableTimestep{TIMES}, T, N, ti}, idxs::Union{TimestepIndex, AnyIndex}...) where {T, N, ti, TIMES}
	_single_index_check(arr.data, idxs)
	idxs1, ts, idxs2 = split_indices(idxs, ti)
	t = ts.index
	_index_bounds_check(arr.data, ti, t)
	return arr.data[idxs1..., t, idxs2...]
end

function Base.setindex!(arr::TimestepArray{FixedTimestep{FIRST, STEP}, T, N, ti}, val, idxs::Union{FixedTimestep{FIRST, STEP, LAST}, AnyIndex}...) where {T, N, ti, FIRST, STEP, LAST}
	idxs1, ts, idxs2 = split_indices(idxs, ti)
	setindex!(arr.data, val, idxs1..., ts.t, idxs2...)
end

function Base.setindex!(arr::TimestepArray{VariableTimestep{TIMES}, T, N, ti}, val, idxs::Union{VariableTimestep{TIMES}, AnyIndex}...) where {T, N, ti, TIMES}
	idxs1, ts, idxs2 = split_indices(idxs, ti)
	setindex!(arr.data, val, idxs1..., ts.t, idxs2...)
end

function Base.setindex!(arr::TimestepArray{FixedTimestep{D_FIRST, STEP}, T, N, ti}, val, idxs::Union{FixedTimestep{T_FIRST, STEP, LAST}, AnyIndex}...) where {T, N, ti, D_FIRST, T_FIRST, STEP, LAST}
	idxs1, ts, idxs2 = split_indices(idxs, ti)
	t = ts.t + findfirst(isequal(T_FIRST[1]), D_FIRST) - 1
	setindex!(arr.data, val, idxs1..., t, idxs2...)
end

function Base.setindex!(arr::TimestepArray{VariableTimestep{D_TIMES}, T, N, ti}, val, idxs::Union{VariableTimestep{T_TIMES}, AnyIndex}...) where {T, N, ti, D_TIMES, T_TIMES}
	idxs1, ts, idxs2 = split_indices(idxs, ti)
	t = ts.t + findfirst(isequal(T_FIRST[1]), T_TIMES) - 1
	setindex!(arr.data, val, idxs1..., t, idxs2...)
end

function Base.setindex!(arr::TimestepArray{FixedTimestep{FIRST, STEP}, T_data, N, ti}, val, idxs::Union{TimestepValue{T_time}, AnyIndex}...) where {T_data, N, ti, FIRST, STEP, T_time}
	_single_index_check(arr.data, idxs)
	idxs1, ts, idxs2 = split_indices(idxs, ti)
	LAST = FIRST + ((size(arr.data, ti) - 1) * STEP)
	t = _get_time_value_position([FIRST:STEP:LAST...], ts)
	setindex!(arr.data, val, idxs1..., t, idxs2...)
end

function Base.setindex!(arr::TimestepArray{VariableTimestep{TIMES}, T_data, N, ti}, val, idxs::Union{TimestepValue{T_time}, AnyIndex}...) where {T_data, N, ti, TIMES, T_time}
	_single_index_check(arr.data, idxs)
	idxs1, ts, idxs2 = split_indices(idxs, ti)
	t = _get_time_value_position(TIMES, ts)
	setindex!(arr.data, val, idxs1..., t, idxs2...)
end

function Base.setindex!(arr::TimestepArray{FixedTimestep{FIRST, STEP}, T, N, ti}, val, idxs::Union{TimestepIndex, AnyIndex}...) where {T, N, ti, FIRST, STEP}
	idxs1, ts, idxs2 = split_indices(idxs, ti)
	t = ts.index
	setindex!(arr.data, val, idxs1..., t, idxs2...)
end

function Base.setindex!(arr::TimestepArray{VariableTimestep{TIMES}, T, N, ti}, val, idxs::Union{TimestepIndex, AnyIndex}...) where {T, N, ti, TIMES}
	idxs1, ts, idxs2 = split_indices(idxs, ti)
	t = ts.index
	setindex!(arr.data, val, idxs1..., t, idxs2...)
end

# Colon support - this allows the time dimension to be indexed with a colon, and
# the deprecation warning will become an error when integer indexing is fully deprecated

function Base.getindex(arr::TimestepArray{FixedTimestep{FIRST, STEP}, T, N, ti}, idxs::AnyIndex...) where {FIRST, STEP, T, N, ti}
    isa(idxs[ti], AnyIndex_NonColon) ? _throw_int_getindex_depwarning() : nothing
    return arr.data[idxs...]
end

function Base.getindex(arr::TimestepArray{VariableTimestep{TIMES}, T, N, ti}, idxs::AnyIndex...) where {TIMES, T, N, ti}
    isa(idxs[ti], AnyIndex_NonColon) ? _throw_int_getindex_depwarning() : nothing
    return arr.data[idxs...]
end

function Base.setindex!(arr::TimestepArray{FixedTimestep{FIRST, STEP}, T, N, ti}, val, idxs::AnyIndex...) where {FIRST, STEP, T, N, ti}
    isa(idxs[ti], AnyIndex_NonColon) ? _throw_int_getindex_depwarning() : nothing
    setindex!(arr.data, val, idxs...)
end

function Base.setindex!(arr::TimestepArray{VariableTimestep{TIMES}, T, N, ti}, val, idxs::AnyIndex...) where {TIMES, T, N, ti}
    isa(idxs[ti], AnyIndex_NonColon) ? _throw_int_getindex_depwarning() : nothing
    setindex!(arr.data, val, idxs...)
end

# Indexing with arrays of TimestepIndexes or TimestepValues
function Base.getindex(arr::TimestepArray{TS, T, N, ti}, idxs::Union{Array{TimestepIndex,1}, AnyIndex}...) where {TS, T, N, ti}
	idxs1, ts_array, idxs2 = split_indices(idxs, ti)
	ts_idxs = _get_ts_indices(ts_array)
	return arr.data[idxs1..., ts_idxs, idxs2...]
end

function Base.getindex(arr::TimestepArray{FixedTimestep{FIRST, STEP}, T_data, N, ti}, idxs::Union{Array{TimestepValue{T_time},1}, AnyIndex}...) where {T_data, N, ti, FIRST, STEP, T_time}
	idxs1, ts_array, idxs2 = split_indices(idxs, ti)
	LAST = FIRST + ((length(arr.data)-1) * STEP)
	ts_idxs = _get_ts_indices(ts_array, [FIRST:STEP:LAST...])
	return arr.data[idxs1..., ts_idxs, idxs2...]
end

function Base.getindex(arr::TimestepArray{VariableTimestep{TIMES}, T_data, N, ti}, idxs::Union{Array{TimestepValue{T_times},1}, AnyIndex}...) where {T_data, N, ti, TIMES, T_times}
	idxs1, ts_array, idxs2 = split_indices(idxs, ti)
	ts_idxs = _get_ts_indices(ts_array, TIMES)
	return arr.data[idxs1..., ts_idxs, idxs2...]
end

function Base.setindex!(arr::TimestepArray{TS, T, N, ti}, vals, idxs::Union{Array{TimestepIndex,1}, AnyIndex}...) where {TS, T, N, ti}
	idxs1, ts_array, idxs2 = split_indices(idxs, ti)
	ts_idxs = _get_ts_indices(ts_array)
	setindex!(arr.data, vals, idxs1..., ts_idxs, idxs2...)
end

function Base.setindex!(arr::TimestepArray{FixedTimestep{FIRST, STEP}, T_data, N, ti}, vals, idxs::Union{Array{TimestepValue{T_times},1}, AnyIndex}...) where {T_data, N, ti, FIRST, STEP, T_times}
	idxs1, ts_array, idxs2 = split_indices(idxs, ti)
	LAST = FIRST + ((length(arr.data)-1) * STEP)
	ts_idxs = _get_ts_indices(ts_array, [FIRST:STEP:LAST...])
	setindex!(arr.data, vals, idxs1..., ts_idxs, idxs2...)
end

function Base.setindex!(arr::TimestepArray{VariableTimestep{TIMES}, T_data, N, ti}, vals, idxs::Union{Array{TimestepValue{T_times},1}, AnyIndex}...) where {T_data, N, ti, TIMES, T_times}
	idxs1, ts_array, idxs2 = split_indices(idxs, ti)
	ts_idxs = _get_ts_indices(ts_array, TIMES)
	setindex!(arr.data, vals, idxs1..., ts_idxs, idxs2...)
end

"""
	hasvalue(arr::TimestepArray, ts::FixedTimestep)

Return `true` or `false`, `true` if the TimestepArray `arr` contains the Timestep `ts`.
"""
function hasvalue(arr::TimestepArray{FixedTimestep{FIRST, STEP}, T, N, ti}, ts::FixedTimestep{FIRST, STEP, LAST}) where {T, N, ti, FIRST, STEP, LAST}
	return 1 <= ts.t <= size(arr, 1)
end

"""
	hasvalue(arr::TimestepArray, ts::VariableTimestep)

Return `true` or `false`, `true` if the TimestepArray `arr` contains the Timestep `ts`.
"""
function hasvalue(arr::TimestepArray{VariableTimestep{TIMES}, T, N, ti}, ts::VariableTimestep{TIMES}) where {T, N, ti, TIMES}
	return 1 <= ts.t <= size(arr, 1)
end

function hasvalue(arr::TimestepArray{FixedTimestep{D_FIRST, STEP}, T, N, ti}, ts::FixedTimestep{T_FIRST, STEP, LAST}) where {T, N, ti, D_FIRST, T_FIRST, STEP, LAST}
	return D_FIRST <= gettime(ts) <= last_period(arr)
end

function hasvalue(arr::TimestepArray{VariableTimestep{D_FIRST}, T, N, ti}, ts::VariableTimestep{T_FIRST}) where {T, N, ti, T_FIRST, D_FIRST}
	return D_FIRST[1] <= gettime(ts) <= last_period(arr)
end

"""
	hasvalue(arr::TimestepArray, ts::FixedTimestep, idxs::Int...)

Return `true` or `false`, `true` if the TimestepArray `arr` contains the Timestep `ts` within
indices `idxs`. Used when Array and Timestep have different FIRST, validating all dimensions.
"""
function hasvalue(arr::TimestepArray{FixedTimestep{D_FIRST, STEP}, T, N, ti},
	ts::FixedTimestep{T_FIRST, STEP, LAST},
	idxs::Int...) where {T, N, ti, D_FIRST, T_FIRST, STEP, LAST}
	return D_FIRST <= gettime(ts) <= last_period(arr) && all([1 <= idx <= size(arr, i) for (i, idx) in enumerate(idxs)])
end

"""
	hasvalue(arr::TimestepArray, ts::VariableTimestep, idxs::Int...)

Return `true` or `false`, `true` if the TimestepArray `arr` contains the Timestep `ts` within
indices `idxs`. Used when Array and Timestep different TIMES, validating all dimensions.
"""
function hasvalue(arr::TimestepArray{VariableTimestep{D_FIRST}, T, N, ti},
	ts::VariableTimestep{T_FIRST},
	idxs::Int...) where {T, N, ti, D_FIRST, T_FIRST}

	return D_FIRST[1] <= gettime(ts) <= last_period(arr) && all([1 <= idx <= size(arr, i) for (i, idx) in enumerate(idxs)])
end


text:
struct Stack
    items :: Vector{Any}
    marks :: Vector{Int}
end

Stack() = Stack([], Int[])

function Base.length(s::Stack)
    length(s.items)
end

function Base.push!(s::Stack, x)
    push!(s.items, x)
    s
end

function Base.pop!(s::Stack)
    pop!(s.items)
end

function poptomark!(s::Stack)
    m = unmark!(s)
    ans = s.items[m+1:end]
    resize!(s.items, m)
    ans
end

function top(s::Stack)
    s.items[end]
end

function mark!(s::Stack)
    push!(s.marks, length(s.items))
    s
end

function unmark!(s::Stack)
    pop!(s.marks)
end

function topmark(s::Stack)
    s.marks[end]
end

ext: jl


text:
type RaggedMatrix{T} <: AbstractMatrix{T}
    data::Vector{T} # a Vector of non-zero entries
    cols::Vector{Int} # a Vector specifying the first index of each column
    m::Int #Number of rows
    function RaggedMatrix(data::Vector{T},cols::Vector{Int},m::Int)
        # make sure the cols are monitonically increasing
        @assert 1==cols[1]
        for j=1:length(cols)-1
            @assert cols[j] ≤ cols[j+1]
        end
        @assert cols[end] == length(data)+1

        # make sure we have less entries than the size of the matrix
        @assert length(data) ≤ m*(length(cols)-1)

        new(data,cols,m)
    end
end

RaggedMatrix(dat::Vector,cols::Vector{Int},m::Int) =
    RaggedMatrix{eltype(dat)}(dat,cols,m)

RaggedMatrix{T}(::Type{T},m::Int,colns::AbstractVector{Int}) =
    RaggedMatrix(Array(T,sum(colns)),[1;1+cumsum(colns)],m)

RaggedMatrix(m::Int,collengths::AbstractVector{Int}) = RaggedMatrix(Float64,m,collengths)

Base.size(A::RaggedMatrix) = (A.m,length(A.cols)-1)

colstop(A::RaggedMatrix,j::Integer) = min(A.cols[j+1]-A.cols[j],size(A,1))

Base.linearindexing{RM<:RaggedMatrix}(::Type{RM}) = Base.LinearSlow()

function getindex(A::RaggedMatrix,k::Int,j::Int)
    if k>size(A,1) || k < 1 || j>size(A,2) || j < 1
        throw(BoundsError(A,(k,j)))
    end

    if A.cols[j]+k-1 < A.cols[j+1]
        A.data[A.cols[j]+k-1]
    else
        zero(eltype(A))
    end
end

function Base.setindex!(A::RaggedMatrix,v,k::Int,j::Int)
    if k>size(A,1) || k < 1 || j>size(A,2) || j < 1
        throw(BoundsError(A,(k,j)))
    end

    if A.cols[j]+k-1 < A.cols[j+1]
        A.data[A.cols[j]+k-1]=v
    elseif v ≠ 0
        throw(BoundsError(A,(k,j)))
    end
    v
end

function Base.convert(::Type{Matrix},A::RaggedMatrix)
    ret = zeros(eltype(A),size(A,1),size(A,2))
    for j=1:size(A,2)
        ret[1:colstop(A,j),j] = view(A,1:colstop(A,j),j)
    end
    ret
end

Base.full(A::RaggedMatrix) = convert(Matrix,A)

function Base.convert(::Type{RaggedMatrix},B::BandedMatrix)
    l = bandwidth(B,1)
    ret = rzeros(eltype(B),size(B,1),l+(1:size(B,2)))
    for j=1:size(B,2),k=colrange(B,j)
        ret[k,j] = B[k,j]
    end
    ret
end

function Base.convert(::Type{RaggedMatrix},B::AbstractMatrix)
    l = bandwidth(B,1)
    ret = rzeros(eltype(B),size(B,1),Int[colstop(B,j) for j=1:size(B,2)])
    for j=1:size(B,2),k=colrange(B,j)
        ret[k,j] = B[k,j]
    end
    ret
end

for (op,bop) in ((:(Base.rand),:rrand),(:(Base.zeros),:rzeros),(:(Base.ones),:rones))
    @eval begin
        $bop{T}(::Type{T},m::Int,colns::AbstractVector{Int}) =
            RaggedMatrix($op(T,sum(colns)),[1;1+cumsum(colns)],m)
        $bop(m::Int,colns::AbstractVector{Int}) = $bop(Float64,m,colns)
    end
end

## BLAS

function Base.A_mul_B!(y::Vector,A::RaggedMatrix,b::Vector)
    m=size(A,2)

    if m ≠ length(b) || size(A,1) ≠ length(y)
        throw(BoundsError())
    end
    T=eltype(y)
    fill!(y,zero(T))
    for j=1:m
        kr=A.cols[j]:A.cols[j+1]-1
        BLAS.axpy!(b[j],view(A.data,kr),view(y,1:length(kr)))
    end
    y
end

function BLAS.axpy!(a,X::RaggedMatrix,Y::RaggedMatrix)
    if size(X) ≠ size(Y)
        throw(BoundsError())
    end

    if X.cols == Y.cols
        BLAS.axpy!(a,X.data,Y.data)
    else
        error("Not implemented.")
    end
end

colstop{T}(X::SubArray{T,2,RaggedMatrix{T},Tuple{UnitRange{Int},UnitRange{Int}}},
        j::Integer) = min(colstop(parent(X),j + first(parentindexes(X)[2])-1) -
                                            first(parentindexes(X)[1]) + 1,
                            size(X,1))

function BLAS.axpy!{T}(a,X::RaggedMatrix,
                       Y::SubArray{T,2,RaggedMatrix{T},Tuple{UnitRange{Int},UnitRange{Int}}})
    if size(X) ≠ size(Y)
        throw(BoundsError())
    end

    for j=1:size(X,2)
        @assert colstop(Y,j) ≤ colstop(X,j)
    end

    P = parent(Y)
    ksh = first(parentindexes(Y)[1]) - 1  # how much to shift
    jsh = first(parentindexes(Y)[2]) - 1  # how much to shift

    for j=1:size(X,2)
        kr = X.cols[j]:X.cols[j+1]-1
        BLAS.axpy!(a,view(X.data,kr),
                    view(P.data,(P.cols[j + jsh] + ksh-1) + (1:length(kr))))
    end

    Y
end

function *(A::RaggedMatrix,B::RaggedMatrix)
    cols = zeros(Int,size(B,2))
    T = promote_type(eltype(A),eltype(B))
    for j=1:size(B,2),k=1:colstop(B,j)
        cols[j] = max(cols[j],colstop(A,k))
    end

    unsafe_A_mul_B!(RaggedMatrix(T,size(A,1),cols),A,B)
end

function unsafe_A_mul_B!(Y::RaggedMatrix,A::RaggedMatrix,B::RaggedMatrix)
    fill!(Y.data,0)

    for j=1:size(B,2),k=1:colstop(B,j)
        BLAS.axpy!(B[k,j],view(A,1:colstop(A,k),k),view(Y.data,Y.cols[j]-1+(1:colstop(A,k))))
    end

    Y
end

function Base.A_mul_B!(Y::RaggedMatrix,A::RaggedMatrix,B::RaggedMatrix)
    for j=1:size(B,2)
        col = 0
        for k=1:colstop(B,j)
            col = max(col,colstop(A,k))
        end

        if col > colstop(Y,j)
            throw(BoundsError())
        end
    end

    unsafe_A_mul_B!(Y,A,B)
end


text:
import Knet.Ops20: bmm, bmm!
using Knet.KnetArrays: DevArray
using CUDA.CUBLAS: CUBLAS, cublasSgemmStridedBatched, cublasDgemmStridedBatched, cublasOperation_t #, handle

function bmm(A::R, B::R; transA::Bool = false, transB::Bool = false) where {R<:DevArray}
    sa, sb = size(A), size(B)
    m, k   = transA ? (sa[2],sa[1]) : (sa[1],sa[2])
    kb, n  = transB ? (sb[2],sb[1]) : (sb[1],sb[2])
    @assert kb == k && sa[3:end]==sb[3:end]
    a3, b3 = reshape(A,sa[1],sa[2],:), reshape(B,sb[1],sb[2],:)
    C = similar(A, m, n, size(a3,3))
    bmm!((transA ? 'T' : 'N'), (transB ? 'T' : 'N'), 1, a3, b3, 0, C)
    reshape(C,m,n,sb[3:end]...)
end

function bmm!(transA::AbstractChar, transB::AbstractChar, alpha::Number, A::R, B::R, beta::Number, C::R) where {R<:DevArray}
    cublasop(c::Char)=cublasOperation_t(if c==='N'; 0; elseif c==='T'; 1; elseif c==='C'; 2; else error("Unknown cublas op $c"); end)
    if ndims(A) != 3 || ndims(B) != 3
        throw(DimensionMismatch("$(map(size,(A,B,C)))"))
    end
    ma,ka,bsa = size(A)
    kb,nb,bsb = size(B)
    if bsa != bsb
        throw(DimensionMismatch("$(map(size,(A,B,C)))"))
    end
    bs = bsa
    if transA == 'N'
        m=ma; k=ka;
    else
        m=ka; k=ma;
    end
    if transB == 'N'
        k == kb || throw(DimensionMismatch("$(map(size,(A,B,C)))"))
        n=nb; k==kb;
    else
        k == nb || throw(DimensionMismatch("$(map(size,(A,B,C)))"))
        n=kb; k==nb;
    end
    (m == size(C,1) && n == size(C,2) && bs == size(C,3)) || throw(DimensionMismatch("$(map(size,(A,B,C)))"))
    lda,ldb,ldc= ma,kb,m
    transa = cublasop(transA); transb = cublasop(transB)
    T = eltype(A)
    alpha = T[alpha]; beta = T[beta]
    strideA, strideB, strideC = m*k, k*n, m*n
    if T<:Float64
        cublasDgemmStridedBatched(CUBLAS.handle(), transa, transb, m, n, k, alpha, A, lda, strideA, B, ldb, strideB, beta, C, ldc, strideC, bs)
    elseif T<:Float32
        cublasSgemmStridedBatched(CUBLAS.handle(), transa, transb, m, n, k, alpha, A, lda, strideA, B, ldb, strideB, beta, C, ldc, strideC, bs)
    else
        error("cublasXgemmStridedBatched does not support $T")
    end
    return C
end


text:
@testset "451.sort-characters-by-frequency.jl" begin
    @test frequency_sort("tree") == "eetr"
    @test frequency_sort("cccaaa") == "cccaaa"
    @test frequency_sort("Aabb") == "bbAa"
end


text:
# This is part of Calendars.jl. See the copyright note there.
# ====================== Gregorian dates ====================

# Day number of the start of the Gregorian calendar.
# CE-0001-01-01 ~ DN#1 = EpochGregorian.
# CE-1111-11-11 ~ DN#0405733
# CE-9999-12-31 ~ DN#3652059
const EpochGregorian = 1
const ValidYearsGregorian = (0, 9999) # by convention

# Is a divisible by b?
divisible(a, b) = rem(a, b) == 0

# True if year is a Gregorian leap year.
function isLeapYearGregorian(year::DPart)::Bool
    (divisible(year, 4) && ! divisible(year, 100)) || divisible(year, 400)
end

# Return the last day of the month for the Gregorian calendar.
function LastDayOfMonthGregorian(year::DPart, month::DPart)
    leap = isLeapYearGregorian(year)
    (month == 2 && leap) && return 29
    return MonthLength[month]
end

# Is the date a valid Gregorian date?
function isValidDateGregorian(cd::CDate, warn=true)
    cal, year, month, day = cd
    # val = (StringToSymbol(cal) == CE
    val = ( cal == CE
        && (ValidYearsGregorian[1] <= year <= ValidYearsGregorian[2])
        && (1 <= month && month <= 12)
        && (1 <= day <= LastDayOfMonthGregorian(year, month)))
    ! val && warn && @warn(Warning(cd))
    return val
end

# Represents the day number a valid Gregorian date?
function isValidDateGregorian(dn::DPart, warn=true)
    # val = EpochGregorian - 2 <= dn
    val = -1 <= dn
    ! val && warn && @warn(Warning(dn))
    return val
end

# Return the days this year so far.
function DayOfYearGregorian(year::DPart, month::DPart, day::DPart)
    day += div(367 * month - 362, 12)
    day += ((month <= 2) ? 0 : (isLeapYearGregorian(year) ? -1 : - 2))
    return day - 1
end

# Computes the day number from a valid Gregorian date.
function DayNumberGregorian(year::DPart, month::DPart, day::DPart)
    days = DayOfYearGregorian(year, month, day)

    dn = (EpochGregorian
        + days                 # days this year
        + 365 * (year - 1)     # days in previous years ignoring leap days
        + div(year - 1, 4)     # Julian leap days before this year...
        - div(year - 1, 100)   # ...minus prior century years...
        + div(year - 1, 400)   # ...plus prior years divisible by 400
    )
    return dn
end
DayNumberGregorian(d::CDate) = DayNumberGregorian(d[2], d[3], d[4])

# Return the Gregorian year from the day number.
function YearGregorian(dn::DPart)
    d = dn - EpochGregorian
    n400, d = divrem(d, 146097) # 400 years contain 146097 days
    n100, d = divrem(d, 36524)
    n4, d   = divrem(d, 1461)
    n1, d   = divrem(d, 365)
    year = 400 * n400 + 100 * n100 + 4 * n4 + n1
    (n100 == 4 || n1 == 4) && return year
    return year + 1
end

function DayOfYearGregorian(dn::DPart, year::DPart)
    mar = DayNumberGregorian(year, 3, 1)
    day = dn - YearStartGregorian(year)
    day += dn < mar ? 0 : (isLeapYearGregorian(year) ? 1 :  2)
    return day
end

# Computes the Gregorian date from a day number.
function DateGregorian(dn::DPart)
    ! isValidDateGregorian(dn) && return InvalidDate

    # Patch in the two allowed exceptions.
    dn == -1 && return (CE, 0, 12, 30)
    dn ==  0 && return (CE, 0, 12, 31)

    year = YearGregorian(dn)
    pdays = DayOfYearGregorian(dn, year)
    month = div(12 * pdays + 373, 367)
    day = dn - DayNumberGregorian(year, month, 1) + 1
    return (CE, year, month, day)::CDate
end

# Return the day number of the first day in the given Gregorian year.
function YearStartGregorian(year::DPart)
    return DayNumberGregorian(year, 1, 1)
end

# Return the day number of the last day in the given Gregorian year.
function YearEndGregorian(year::DPart)
    return DayNumberGregorian(year + 1, 1, 1) - 1
end

# Return the number of months in the given Gregorian year.
function MonthsInYearGregorian(year::DPart)
    return 12
end

# Return the number of days in the given Gregorian year.
function DaysInYearGregorian(year::DPart)
    return YearStartGregorian(year + 1) - YearStartGregorian(year)
end


text:
using
    Test,
    SeawaterPolynomials,
    SeawaterPolynomials.SecondOrderSeawaterPolynomials,
    SeawaterPolynomials.TEOS10

using SeawaterPolynomials: AbstractSeawaterPolynomial, BoussinesqEquationOfState

""" Test instantiation of a RoquetSeawaterPolynomial."""
function instantiate_roquet_polynomial(FT, coefficient_set)
    polynomial = RoquetSeawaterPolynomial(FT, coefficient_set)
    return typeof(polynomial) <: AbstractSeawaterPolynomial
end

""" Test instantiation of a RoquetEquationOfState."""
function instantiate_roquet_equation_of_state(FT, coefficient_set)
    eos = RoquetEquationOfState(FT, coefficient_set)
    return typeof(eos) <: BoussinesqEquationOfState
end

""" Test instantiation of TEOS10SeawaterPolynomial."""
function instantiate_teos10_polynomial(FT)
    polynomial = TEOS10SeawaterPolynomial(FT)
    return typeof(polynomial) <: AbstractSeawaterPolynomial
end

""" Test instantiation of a RoquetEquationOfState."""
function instantiate_teos10_equation_of_state(FT)
    eos = TEOS10EquationOfState(FT)
    return typeof(eos) <: BoussinesqEquationOfState
end

@testset "Second-order seawater polynomials" begin
    for coefficient_set in (
                            :Linear,
                            :Cabbeling,
                            :CabbelingThermobaricity,
                            :Freezing,
                            :SecondOrder
                           )

        for FT in (Float64, Float32)
            @test instantiate_roquet_polynomial(FT, coefficient_set)
            @test instantiate_roquet_equation_of_state(FT, coefficient_set)

            eos = RoquetEquationOfState(FT)

            @test SeawaterPolynomials.ρ′(0, 0, 0, eos) == 0
            @test SeawaterPolynomials.haline_contraction(0, 0, 0, eos) == eos.seawater_polynomial.R₁₀₀
            @test SeawaterPolynomials.thermal_expansion(0, 0, 0, eos) == eos.seawater_polynomial.R₀₁₀
        end
    end

end

@testset "TEOS-10 seawater polynomials" begin
    for FT in (Float64, Float32)
        @test instantiate_teos10_polynomial(FT)
        @test instantiate_teos10_equation_of_state(FT)

        eos = TEOS10EquationOfState(FT)
        # Test/check values from Roquet et al. (2014).
        Θ = 10   # [C]
        S = 30   # [g/kg]
        Z = 1e3  # [m]

        τ = SeawaterPolynomials.TEOS10.τ(Θ)
        s = SeawaterPolynomials.TEOS10.s(S)
        ζ = SeawaterPolynomials.TEOS10.ζ(Z)

        @test SeawaterPolynomials.TEOS10.r₀(ζ) ≈ 4.59763035
        @test SeawaterPolynomials.TEOS10.r′(τ, s, ζ) ≈ 1022.85377

        @test SeawaterPolynomials.TEOS10.ρ′(Θ, S, Z, eos) ≈ 1027.45140
        @test SeawaterPolynomials.TEOS10.thermal_expansion(Θ, S, Z, eos) ≈ 0.179646281
        @test SeawaterPolynomials.TEOS10.haline_contraction(Θ, S, Z, eos) ≈ 0.765555368
    end
end


text:
using Test
using Random: randn
using AutoGrad: @gcheck, Param
using Knet.Ops21: gelu, mmul

# Allow the user to override this:
if !isdefined(Main, :testparam)
    testparam(x...) = Param(randn(x...))
end

@testset "ops21" begin

    @testset "activation" begin
        x4 = testparam(5,4,3,2)
        @gcheck gelu.(x4)
    end

    @testset "mmul" begin
        x4 = testparam(5,4,3,2)
        w4 = testparam(6,5,4,3)
        @gcheck mmul(w4,x4,dims=3)
        x3 = testparam(3,4,5)
        w3 = testparam(5,4,3)
        @gcheck mmul(w3,x3)
    end
end


text:
module MixedBags

include("core.jl")
include("operations.jl")

end # module


text:
# Hash table acceleration index
struct UniqueHashIndex{D <: Dictionary} <: AbstractUniqueIndex
    dict::D
end

function UniqueHashIndex(a::AbstractArray)
    dict = Dictionary{eltype(a), SingleVector{eltype(keys(a))}}()

    @inbounds for i in keys(a)
        value = a[i]
        (hadindex, token) = gettoken!(dict, value)
        if hadindex
            error("Input not unique") # TODO Use appropriate Exception
        else # `value` is ready to be inserted into `dict` at slot `-index`
            @inbounds settokenvalue!(dict, token, SingleVector(i))
        end
    end
    return UniqueHashIndex(dict)
end

Base.summary(::UniqueHashIndex) = "UniqueHashIndex"

# Accelerations
Base.in(x, a::AcceleratedArray{<:Any, <:Any, <:Any, <:UniqueHashIndex}) = haskey(a.index.dict, x)

function Base.count(f::Fix2{typeof(isequal)}, a::AcceleratedArray{<:Any, <:Any, <:Any, <:UniqueHashIndex})
    if f.x in keys(a.index.dict)
        return 1
    else
        return 0
    end
end

function Base.findall(f::Fix2{typeof(isequal)}, a::AcceleratedArray{<:Any, <:Any, <:Any, <:UniqueHashIndex})
    (hasindex, token) = gettoken(a.index.dict, f.x)
    if hasindex
        return MaybeVector(@inbounds gettokenvalue(a.index.dict, token)[])
    else
        return MaybeVector{eltype(keys(a))}()
    end
end

function Base.filter(f::Fix2{typeof(isequal)}, a::AcceleratedArray{<:Any, <:Any, <:Any, <:UniqueHashIndex})
    if f.x in keys(a.index.dict)
        return MaybeVector{eltype(a)}(f.x)
    else
        return MaybeVector{eltype(a)}()
    end
end

function SplitApplyCombine.group(a::AcceleratedArray{<:Any, <:Any, <:Any, <:UniqueHashIndex}, b::AbstractArray)
    return map(inds -> SingleVector(@inbounds b[inds[]]), a.index.dict)
end

function SplitApplyCombine.groupreduce(::typeof(identity), f, op, a::AcceleratedArray{<:Any, <:Any, <:Any, <:UniqueHashIndex}; kw...)
    return map(inds -> @inbounds(reduce(op, a[inds[]]; kw...)), a.index.dict)
end

function SplitApplyCombine.groupfind(a::AcceleratedArray{<:Any, <:Any, <:Any, <:UniqueHashIndex})
    return a.index.dict
end

function SplitApplyCombine._innerjoin!(out, left::AbstractArray, right::AcceleratedArray{<:Any, <:Any, <:Any, <:UniqueHashIndex}, v::AbstractArray, ::typeof(isequal))
    @boundscheck if (axes(left)..., axes(right)...) != axes(v)
        throw(DimensionMismatch("innerjoin arrays do not have matching dimensions"))
    end

    dict = right.index.dict

    @inbounds for i_l ∈ keys(left)
        (hasindex, token) = gettoken(dict, left[i_l])
        if hasindex
            i_r = gettokenvalue(dict, token)[]
            push!(out, v[Tuple(i_l)..., Tuple(i_r)...])
        end
    end

    return out
end


text:
module SlurpTypes
import ...PatternStructure.Trees: SlurpHead
export SlurpHead, LazySlurp, GreedySlurp, GenericGreedySlurp, GenericLazySlurp,
       SimpleLastSlurp, SimpleGreedySlurpUntil, SimpleLazySlurpUntil, SimpleSlurp

abstract type LazySlurp   <: SlurpHead end
abstract type GreedySlurp <: SlurpHead end

struct GenericLazySlurp   <: LazySlurp end
struct GenericGreedySlurp <: GreedySlurp end

#-----------------------------------------------------------------------------------
# Simple: slurp is composed of one single binding name, e.g. *{a}
#-----------------------------------------------------------------------------------

abstract type SimpleGreedySlurp <: GreedySlurp end
abstract type SimpleLazySlurp   <: LazySlurp end

struct SimpleLastSlurp <: SimpleGreedySlurp
  post::Int
end

struct SimpleGreedySlurpUntil <: SimpleGreedySlurp
  until::Vector
  index::Int
end

struct SimpleLazySlurpUntil <: SimpleLazySlurp
  until::Vector
  index::Int
end

const SimpleSlurp = Union{SimpleGreedySlurp, SimpleLazySlurp}
end


text:
mutable struct DebugRect
    layoutobservables::GridLayoutBase.LayoutObservables
    width::Observable
    height::Observable
    tellwidth::Observable
    tellheight::Observable
    halign::Observable
    valign::Observable
    leftprot::Observable
    rightprot::Observable
    bottomprot::Observable
    topprot::Observable
end

observablify(x::Observable) = x
observablify(x, type=Any) = Observable{type}(x)

function DebugRect(; bbox = nothing, width=nothing, height=nothing,
    tellwidth = true, tellheight = true, halign=:center,
    valign=:center, topprot=0.0, leftprot=0.0, rightprot=0.0, bottomprot=0.0)

    width = observablify(width)
    height = observablify(height)
    tellwidth = observablify(tellwidth)
    tellheight = observablify(tellheight)
    halign = observablify(halign)
    valign = observablify(valign)
    topprot = observablify(topprot, Float32)
    leftprot = observablify(leftprot, Float32)
    rightprot = observablify(rightprot, Float32)
    bottomprot = observablify(bottomprot, Float32)

    protrusions::Observable{GridLayoutBase.RectSides{Float32}} = map(leftprot, rightprot, bottomprot, topprot) do l, r, b, t
        GridLayoutBase.RectSides{Float32}(l, r, b, t)
    end

    layoutobservables = GridLayoutBase.LayoutObservables(DebugRect, width,
        height, tellwidth, tellheight, halign, valign; suggestedbbox = bbox, protrusions = protrusions)

    DebugRect(layoutobservables, height, width, tellwidth, tellheight,
        halign, valign, leftprot, rightprot, bottomprot, topprot)
end


text:
export hiddenLayer

type hiddenLayer <:Layer
  bottom ::Array{Layer}
  top    ::Array{Layer}
  topRange ::Vector

  rᵢ     ::Array{Float32} #input rate
  W      ::Array{Float32}
  ∂W     ::Array{Float32}
  b      ::Array{Float32}
  ∂b     ::Array{Float32}
  I      ::Array{Float32}
  A      ::Function # activation function
  ∂A     ::Function # derivative of activation function
  rₒ     ::Array{Float32} # output rate

  ∂rᵢ_L  ::Array{Float32}
  ∂w_L   ::Array{Float32}
  ∂b_L   ::Array{Float32}
  ∂i_L   ::Array{Float32}
  ∂rₒ_L  ::Array{Float32}

  lock   ::Bool
end
hiddenLayer(inDim::Int,outDim::Int)=
  hiddenLayer(Layer[],
              Layer[],
              Vector{Array}(),

              map(Float32,rand(inDim,batchsize)),
              0.1*map(Float32,randn(outDim,inDim)),
              zeros(Float32,outDim,inDim),
              0.1*map(Float32,randn(outDim,1)),
              zeros(Float32,outDim,1),
              map(Float32,randn(outDim,batchsize)),
              recLinear,#sigmoid
              ∂recLinear,#∂sigmoid
              map(Float32,rand(outDim,batchsize)),

              map(Float32,randn(inDim,batchsize)),
              map(Float32,randn(outDim,inDim,batchsize)),
              map(Float32,randn(outDim,batchsize)),
              map(Float32,randn(outDim,batchsize)),
              map(Float32,randn(outDim,batchsize)),

              false)


text:
function init_sin_test!(problem, bl, state, aux, localgeo, t)
    (x, y, z) = localgeo.coord

    FT = eltype(state)
    param_set = parameter_set(bl)

    z = FT(z)
    _grav::FT = grav(param_set)
    _MSLP::FT = MSLP(param_set)

    # These constants are those used by Stevens et al. (2005)
    qref = FT(9.0e-3)
    q_pt_sfc = PhasePartition(qref)
    Rm_sfc = FT(gas_constant_air(param_set, q_pt_sfc))
    T_sfc = FT(292.5)
    P_sfc = _MSLP

    # Specify moisture profiles
    q_liq = FT(0)
    q_ice = FT(0)
    zb = FT(600)         # initial cloud bottom
    zi = FT(840)         # initial cloud top
    dz_cloud = zi - zb
    q_liq_peak = FT(0.00045)     # cloud mixing ratio at z_i

    if z > zb && z <= zi
        q_liq = (z - zb) * q_liq_peak / dz_cloud
    end

    if z <= zi
        θ_liq = FT(289.0)
        q_tot = qref
    else
        θ_liq = FT(297.5) + (z - zi)^(FT(1 / 3))
        q_tot = FT(1.5e-3)
    end

    w = FT(10 + 0.5 * sin(2 * π * ((x / 1500) + (y / 1500))))
    u = (5 + 2 * sin(2 * π * ((x / 1500) + (y / 1500))))
    v = FT(5 + 2 * sin(2 * π * ((x / 1500) + (y / 1500))))

    # Pressure
    H = Rm_sfc * T_sfc / _grav
    p = P_sfc * exp(-z / H)

    # Density, Temperature
    ts = PhaseEquil_pθq(param_set, p, θ_liq, q_tot)
    #ρ = air_density(ts)
    ρ = one(FT)

    e_kin = FT(1 / 2) * FT((u^2 + v^2 + w^2))
    e_pot = _grav * z
    E = ρ * total_energy(e_kin, e_pot, ts)

    state.ρ = ρ
    state.ρu = SVector(ρ * u, ρ * v, ρ * w)
    state.energy.ρe = E
    state.moisture.ρq_tot = ρ * q_tot
end


text:
function generate_test_data(filepath, n=10)
    # Generate a table with data
    tbl = table((id = ["0"*string(i) for i in 1:n],
                 IntField = [1000*i for i in 1:n],
                 FloatField = [i/10.0 for i in 1:n],
                 StringField = join.([collect(Iterators.take(string.('A':'E'), rand(0:10))) for _ in 1:n], " "),
                 StringField2 = join.([collect(Iterators.take(["abc", "abcd", "def", "ghi", "xyz", "αβγ", "ω"], rand(0:3))) for _ in 1:n], " "),
                 RandFloat = rand(n),
                 RandString = [randstring(rand(1:10)) for _ in 1:n],
                ),
                pkey=:id)

    # Write file
    open(filepath, "w") do io
        println(io, join(colnames(tbl), ','))
        for row in collect(tbl)
            println(io, join(row, ','))
        end
    end
end

function generate_test_embeddings(filepath, embtype=:WordVectors; dim=10, kind=:text)
    vocab = vcat(string.(collect(0:1000:10_000)),
                 string.('A':'E'),
                 ["abc", "abcd", "def", "ghi", "xyz", "αβγ", "ω"]
                )
    vectors = rand(dim, length(vocab))
    vectors./= permutedims(map(norm, (vectors[:, i] for i in 1:size(vectors,2))))
    local wv
    if embtype == :WordVectors
        wv = WordVectors(vocab, vectors)
    elseif embtype == :CompressedWordVectors
        wv = compress(WordVectors(vocab, vectors), k=10, m=2)
    elseif embtype == :GloveWordVectors
        wv = WordVectors(push!(vocab, "<unk>"), hcat(vectors, zeros(dim)))
    else
        @error "Unknown embeddings type. Use :WordVectors, :CompressedWordVectors or :GloveWordVectors."
    end
    write2disk(filepath, wv; kind=kind)
    return wv
end