text:
import Base: sum, prod, diff, mean, var, std, quantile, minimum, maximum, median, middle, cov, cor, cumprod, cumsum, cummin, cummax

# define several helper functions such as sum or mean over nullable arrays (AbstractArrayWrapper{Nullable{T}}).

"""

`dropnaiter(arr)`

Generate an iterator from a nullable array `arr`, which iterates over only non-null elements.

##### Examples

```julia
julia> for x in dropnaiter(@nalift([1,2,NA,4,5]))
         println(x)
       end
1
2
4
5
```

"""
dropnaiter{T,N}(arr::AbstractArray{Nullable{T},N}) = NonNAIterator{T,N,typeof(arr)}(arr)

"""

`enum_dropnaiter(arr)`

Generate an iterator from a nullable array `arr`, which yields (index, elem) for an integer `index` for non-null element positions of `arr` and a non-null element `elem`.

##### Examples

```julia
julia> for x in enum_dropnaiter(@nalift([:A,:B,NA,NA,:C]))
         println(x)
       end
(1,:A)
(2,:B)
(5,:C)
```

"""
enum_dropnaiter{T,N}(arr::AbstractArray{Nullable{T},N}) = EnumerateNonNAIterator{T,N,typeof(arr)}(arr)

immutable NonNAIterator{T,N,A<:AbstractArray}
  array::A #AbstractArray{Nullable{T},N}
end
Base.eltype{T,N,A}(::Type{NonNAIterator{T,N,A}}) = T

Base.start(iter::NonNAIterator) = start(iter.array)
Base.next{T,N,U,A}(iter::NonNAIterator{T,N,A}, state::U) = begin
  (next_item::Nullable{T}, next_state::U) = next(iter.array, state)
  while next_item.isnull
    (next_item, next_state) = next(iter.array, next_state)
  end
  (next_item.value, next_state)::Tuple{T,U}
  # recursion does not increase/decrease performance in some case.
  #if next_item.isnull
  #  # I need to check if TCO is applied here.
  #  next(iter, next_state)::Tuple{T,U}
  #else
  #  (next_item.value, next_state)::Tuple{T,U}
  #end
end
Base.done{T,N,U,A}(iter::NonNAIterator{T,N,A}, state::U) = begin
  isdonenow::Bool = done(iter.array, state)
  if isdonenow
    true
  else
    (next_item, next_state) = next(iter.array, state)
    if next_item.isnull
      done(iter, next_state)::Bool
    else
      false
    end
  end
end

immutable EnumerateNonNAIterator{T,N,A<:AbstractArray}
  array::A #::AbstractArray{Nullable{T},N}
end

Base.eltype{T,N,A}(::Type{EnumerateNonNAIterator{T,N,A}}) = Tuple{Int,T}
Base.start{T,N}(iter::EnumerateNonNAIterator{T,N}) = (1, start(iter.array))
Base.next{T,N}(iter::EnumerateNonNAIterator{T,N}, state) = begin
  (next_item, next_state) = next(iter.array, state[2])
  if next_item.isnull
    next(iter, (state[1]+1, next_state)) #::Tuple{Tuple{Int,T},Tuple{Int,typeof(state)}}
  else
    ((state[1], next_item.value), (state[1]+1, next_state)) #::Tuple{Tuple{Int,T},Tuple{Int,typeof(state)}}
  end
end
Base.done{T,N}(iter::EnumerateNonNAIterator{T,N}, state) = begin
  isdonenow = done(iter.array, state[2])
  if isdonenow
    true
  else
    (next_item, next_state) = next(iter.array, state[2])
    if next_item.isnull
      done(iter, (state[1]+1, next_state))::Bool
    else
      false
    end
  end
end

immutable ZipNonNAIterator{T,N,M,A}
  arrays::A
end

"""

`zip_dropnaiter(arrs...)`

Generate a zipped iterator from nullable arrays `arrs...`. If any element in `arrs...` is null, the iterator will skip it and move to the next element tuple.

##### Examples

```julia
julia> for x in zip_dropnaiter(@nalift([11,12,NA,NA,15]),
                               @nalift([:X,NA,:Z,NA,:V]),
                               @nalift([71,72,73,NA,75]))
         println(x)
       end
(11,:X,71)
(15,:V,75)
```

"""
zip_dropnaiter{N}(arrs::AbstractArray{TypeVar(:T),N}...) = begin
  T = Tuple{map(x->eltype(eltype(x)), arrs)...}
  ZipNonNAIterator{T,N,length(arrs),typeof(arrs)}(arrs)
end
Base.eltype{T,N,M,A}(::Type{ZipNonNAIterator{T,N,M,A}}) = T
Base.start{T,N,M}(iter::ZipNonNAIterator{T,N,M}) = map(start, iter.arrays)
Base.next{T,N,M}(iter::ZipNonNAIterator{T,N,M}, state) = begin
  nexts = map(iter.arrays, state) do array,s; next(array,s) end
  for (n,s) in nexts
    if n.isnull
      return next(iter, map(x->x[2],nexts))
    end
  end
  (map(x->x[1].value,nexts), map(x->x[2],nexts))
end
Base.done{T,N,M}(iter::ZipNonNAIterator{T,N,M}, state) = begin
  isdonenow = reduce(|, map(iter.arrays, state) do array,s; done(array,s) end)
  if isdonenow
    return true
  end
  nexts = map(iter.arrays, state) do array,s; next(array,s) end
  for (n,s) in nexts
    if n.isnull
      return done(iter, map(x->x[2],nexts))::Bool
    end
  end
  return false
end

# for some reason, I cannot make dropnaiter type stable.
# for now, just use the plain for loop.
Base.sum{T,N,A}(arr::AbstractArrayWrapper{Nullable{T},N,A}) = begin
  acc = zero(T)
  for x in arr
    if !x.isnull
      acc += x.value
    end
  end
  Nullable(acc)
end
Base.sum{T<:AbstractFloat,N,A}(arr::AbstractArrayWrapper{Nullable{T},N,FloatNAArray{T,N,A}}) = begin
  arrdata = arr.a.data
  acc = zero(T)
  for x in arrdata
    if !isnan(x)
      acc += x
    end
  end
  Nullable(acc)
end

Base.prod{T,N,A}(arr::AbstractArrayWrapper{Nullable{T},N,A}) = begin
  acc = one(T)
  for x in arr
    if !x.isnull
      acc *= x.value
    end
  end
  Nullable(acc)
end
Base.prod{T<:AbstractFloat,N,A}(arr::AbstractArrayWrapper{Nullable{T},N,FloatNAArray{T,N,A}}) = begin
  arrdata = arr.a.data
  acc = one(T)
  for x in arrdata
    if !isnan(x)
      acc *= x
    end
  end
  Nullable(acc)
end

Base.mean{T}(arr::AbstractArrayWrapper{Nullable{T}}) = (iter=dropnaiter(arr);isempty(iter) ? Nullable{T}() : Nullable(mean(iter)))
Base.mean{T<:AbstractFloat,N,A}(arr::AbstractArrayWrapper{Nullable{T},N,FloatNAArray{T,N,A}}) = begin
  acc = zero(T)
  count = 0
  for elem in arr.a.data
    if !isnan(elem)
      count += 1
      acc += elem
    end
  end
  if count == 0
    Nullable{T}()
  else
    Nullable(acc / convert(T,count))
  end
end

Base.var{T}(arr::AbstractArrayWrapper{Nullable{T}}; corrected::Bool=true, mean=nothing) = begin
  iter = dropnaiter(arr)
  isempty(iter) ? Nullable{T}() : Nullable(var(iter; corrected=corrected, mean=mean))
end
Base.std{T}(arr::AbstractArrayWrapper{Nullable{T}}; corrected::Bool=true, mean=nothing) = begin
  iter = dropnaiter(arr)
  isempty(iter) ? Nullable{T}() : Nullable(std(iter; corrected=corrected, mean=mean))
end

# TODO need to optimize this without converting to another arrays first.
Base.cov{T,U}(arr1::AbstractArrayWrapper{Nullable{T}}, arr2::AbstractArrayWrapper{Nullable{U}}) = begin
  zipped = zip_dropnaiter(arr1, arr2)
  v1 = Array(T,length(arr1))
  v2 = Array(U,length(arr2))
  n = 1
  for z in zipped
    v1[n] = z[1]
    v2[n] = z[2]
    n += 1
  end
  resize!(v1, n-1)
  resize!(v2, n-1)
  r = cov(v1, v2)
  if isnan(r)
    Nullable{typeof(r)}()
  else
    Nullable(r)
  end
end

# TODO need to optimize this without converting to another arrays first.
Base.cor{T,U}(arr1::AbstractArrayWrapper{Nullable{T}}, arr2::AbstractArrayWrapper{Nullable{U}}) = begin
  zipped = zip_dropnaiter(arr1, arr2)
  v1 = Array(T,length(arr1))
  v2 = Array(U,length(arr2))
  n = 1
  for z in zipped
    v1[n] = z[1]
    v2[n] = z[2]
    n += 1
  end
  if n == 1
    # well, it's not really correct. If T and U are Int64, the result might be a Float64.
    Nullable{Base.return_types(/, (promote_type(T,U),Int))[1]}()
  else
    resize!(v1, n-1)
    resize!(v2, n-1)
    r = cor(v1, v2)
    if isnan(r)
      Nullable{typeof(r)}()
    else
      Nullable(cor(v1, v2))
    end
  end
end

Base.minimum{T}(arr::AbstractArrayWrapper{Nullable{T}}) = begin
  r=type_array(collect(dropnaiter(arr)))
  isempty(r) ? Nullable{T}() : Nullable(minimum(r))
end
Base.minimum{T<:AbstractFloat,N,A}(arr::AbstractArrayWrapper{Nullable{T},N,FloatNAArray{T,N,A}}) = (r=minimum(arr.a.data);isnan(r) ? Nullable{T}() : Nullable(r))

Base.maximum{T}(arr::AbstractArrayWrapper{Nullable{T}}) = begin
  r=type_array(collect(dropnaiter(arr)))
  isempty(r) ? Nullable{T}() : Nullable(maximum(r))
end
Base.maximum{T<:AbstractFloat,N,A}(arr::AbstractArrayWrapper{Nullable{T},N,FloatNAArray{T,N,A}}) = (r=maximum(arr.a.data);isnan(r) ? Nullable{T}() : Nullable(r))

Base.median{T}(arr::AbstractArrayWrapper{Nullable{T}}) = median_helper(typeof(one(T) / 2), arr)
median_helper{T,DIVTYPE}(::Type{DIVTYPE}, arr::AbstractArrayWrapper{Nullable{T}}) = begin
  r=collect_nonnas(arr)
  isempty(r) ? Nullable{DIVTYPE}() : Nullable(median(r))
end

collect_nonnas{T}(arr::AbstractArrayWrapper{Nullable{T}}) = type_array(collect(dropnaiter(arr)))
collect_nonnas{T<:AbstractFloat,N,A}(arr::AbstractArrayWrapper{Nullable{T},N,FloatNAArray{T,N,A}}) = begin
  nonnas = similar(arr.a.data, T)
  count = 0
  for elem in arr.a.data
    if !isnan(elem)
      count += 1
      nonnas[count] = elem
    end
  end
  resize!(nonnas, count)
  nonnas
end

Base.middle{T}(arr::AbstractArrayWrapper{Nullable{T}}) = middle_helper(typeof(one(T) / 2), arr)
middle_helper{T,DIVTYPE}(::Type{DIVTYPE}, arr::AbstractArrayWrapper{Nullable{T}}) = begin
  r=type_array(collect(dropnaiter(arr)))
  isempty(r) ? Nullable{DIVTYPE}() : Nullable(middle(r))
end
Base.middle{T<:AbstractFloat,N,A}(arr::AbstractArrayWrapper{Nullable{T},N,FloatNAArray{T,N,A}}) = (r=middle(arr.a.data);isnan(r) ? Nullable{T}() : Nullable(r))

Base.quantile{T}(arr::AbstractArrayWrapper{Nullable{T}}, q::AbstractVector) = quantile_helper(typeof(one(T) / 2), arr, q)
quantile_helper{T,DIVTYPE}(::Type{DIVTYPE}, arr::AbstractArrayWrapper{Nullable{T}}, q::AbstractVector) = begin
  r=collect_nonnas(arr)
  if isempty(r)
    Nullable{DIVTYPE}[]
  else
    map(Nullable, quantile(r, q))
  end
end
Base.quantile{T}(arr::AbstractArrayWrapper{Nullable{T}}, q::Number) = quantile_helper(typeof(one(T) / 2), arr, q)
quantile_helper{T,DIVTYPE}(::Type{DIVTYPE}, arr::AbstractArrayWrapper{Nullable{T}}, q::Number) = begin
  r=collect_nonnas(arr)
  if isempty(r)
    Nullable{DIVTYPE}()
  else
    Nullable(quantile(r, q))
  end
end

for op = [:sum, :prod, :mean, :var, :std, :median, :middle, :minimum, :maximum, :quantile]
  @eval begin
    $(op)(arr::DictArray) = mapvalues($(op), arr.data)
    $(op)(arr::LabeledArray) = $(op)(arr.data)
  end
end
Base.quantile(arr::DictArray, q::Number) = mapvalues(v->quantile(v, q), arr.data)
Base.quantile(arr::LabeledArray, q::Number) = quantile(arr.data, q)

for op = [:cov, :cor]
  @eval begin
    $(op)(arr1::DictArray, arr2::DictArray) = begin
      axis1 = keys(arr1)
      axis2 = keys(arr2)
      # the covariance matrix should not be large.
      len1 = length(axis1)
      len2 = length(axis2)
      inds = reshape(collect(zip(repeat(collect(1:len1),inner=[1],outer=[len2]),
                                 repeat(collect(1:len2),inner=[len1],outer=[1]))),len1,len2)
      data = map(inds) do xy
        $(op)(arr1.data.values[xy[1]], arr2.data.values[xy[2]])
      end
      larr(data, axis1=axis1, axis2=axis2)
    end
    $(op)(arr1::LabeledArray, arr2::LabeledArray) = begin
      assert(arr1.axes == arr2.axes)
      $(op)(arr1.data, arr2.data)
    end
    $(op)(arr::DictArray{TypeVar(:T),1}) = $(op)(arr, arr)
    $(op)(arr::DictArray{TypeVar(:T),2}) = $(op)(arr, arr)
    $(op)(arr::DictArray) = $(op)(arr, arr)
    $(op)(arr::LabeledArray{TypeVar(:T),1}) = $(op)(arr, arr)
    $(op)(arr::LabeledArray{TypeVar(:T),2}) = $(op)(arr, arr)
    $(op)(arr::LabeledArray) = $(op)(arr, arr)
  end
end

for opspec in [:cumsum, :cumprod, :cummin, :cummax, :diff]
  @eval begin
    $(opspec){T}(arr::AbstractArrayWrapper{Nullable{T}};rev=false) = $(opspec)(arr, 1;rev=rev)
    $(opspec){T,N,A}(arr::AbstractArrayWrapper{Nullable{T},N,FloatNAArray{T,N,A}};rev=false) = $(opspec)(arr, 1;rev=rev)
    $(opspec)(arr::DictArray, dims::Integer...;rev=false) = DictArray(mapvalues(v->$(opspec)(v, dims...;rev=rev), arr.data))
    $(opspec)(arr::LabeledArray, dims::Integer...;rev=false) = LabeledArray($(opspec)(arr.data, dims...;rev=rev), arr.axes)
  end
end

for opspec in [:msum, :mprod, :mmean, :mminimum, :mmaximum, :mmedian, :mmiddle, :nafill]
  @eval begin
    $(opspec){T}(arr::AbstractArrayWrapper{Nullable{T}};rev=false,window=0) = $(opspec)(arr, 1;rev=rev,window=window)
    $(opspec){T,N,A}(arr::AbstractArrayWrapper{Nullable{T},N,FloatNAArray{T,N,A}};rev=false,window=0) = $(opspec)(arr, 1;rev=rev,window=window)
    $(opspec)(arr::DictArray, dims::Integer...;rev=false,window=0) = DictArray(mapvalues(v->$(opspec)(v, dims...;rev=rev,window=window), arr.data))
    $(opspec)(arr::LabeledArray, dims::Integer...;rev=false,window=0) = LabeledArray($(opspec)(arr.data, dims...;rev=rev,window=window), arr.axes)
  end
end

mquantile{T}(arr::AbstractArrayWrapper{Nullable{T}}, quantile::Number;rev=false,window=0) = mquantile(arr, quantile, 1; rev=rev,window=window)
mquantile{T,N,A}(arr::AbstractArrayWrapper{Nullable{T},N,FloatNAArray{T,N,A}}, quantile::Number;rev=false,window=0) = mquantile(arr, quantile, 1; rev=rev,window=window)
mquantile(arr::DictArray, quantile::Number, dims::Integer...;rev=false,window=0) =
  DictArray(mapvalues(v->mquantile(v, quantile, dims...;rev=rev,window=window), arr.data))
mquantile(arr::LabeledArray, quantile::Number, dims::Integer...;rev=false,window=0) =
  LabeledArray(mquantile(arr.data, quantile, dims...;rev=rev,window=window), arr.axes)

Base.cumsum{T}(arr::AbstractArrayWrapper{Nullable{T}}, dims::Integer...;rev=false) = begin
  function update!(tgt, src)
    acc = zero(T)
    for i in eachindex(src)
      if !src[i].isnull
        acc += src[i].value
      end
      tgt[i] = Nullable(acc)
    end
  end
  result = similar(arr)
  map_array_preserve_shape!(update!, result.a, arr.a, dims...;rev=rev)
  result
end
Base.cumsum{T,N,A}(arr::AbstractArrayWrapper{Nullable{T},N,FloatNAArray{T,N,A}}, dims::Integer...;rev=false) = begin
  function update!(tgt, src)
    acc = zero(T)
    for i in eachindex(src)
      if !isnan(src[i])
        acc += src[i]
      end
      tgt[i] = acc
    end
  end
  result = similar(arr)
  map_array_preserve_shape!(update!, result.a.data, arr.a.data, dims...;rev=rev)
  result
end

Base.cumprod{T}(arr::AbstractArrayWrapper{Nullable{T}}, dims::Integer...;rev=false) = begin
  function update!(tgt, src)
    acc = one(T)
    for i in eachindex(src)
      if !src[i].isnull
        acc *= src[i].value
      end
      tgt[i] = Nullable(acc)
    end
  end
  result = similar(arr)
  map_array_preserve_shape!(update!, result.a, arr.a, dims...;rev=rev)
  result
end
Base.cumprod{T,N,A}(arr::AbstractArrayWrapper{Nullable{T},N,FloatNAArray{T,N,A}},
                    dims::Integer...
                    ;rev=false) = begin
  function update!(tgt, src)
    acc = one(T)
    for i in eachindex(src)
      if !isnan(src[i])
        acc *= src[i]
      end
      tgt[i] = acc
    end
  end
  result = similar(arr)
  map_array_preserve_shape!(update!, result.a.data, arr.a.data, dims...;rev=rev)
  result
end

Base.cummin{T}(arr::AbstractArrayWrapper{Nullable{T}}, dims::Integer...;rev=false) = begin
  function update!(tgt, src)
    acc = Nullable{T}()
    for i in eachindex(src)
      if !src[i].isnull
        acc = acc.isnull ? src[i] : Nullable(min(acc.value, src[i].value))
      end
      tgt[i] = acc
    end
  end
  result = similar(arr)
  map_array_preserve_shape!(update!, result.a, arr.a, dims...;rev=rev)
  result
end
Base.cummin{T,N,A}(arr::AbstractArrayWrapper{Nullable{T},N,FloatNAArray{T,N,A}},
                   dims::Integer...
                   ;rev=false) = begin
  function update!(tgt, src)
    acc = convert(T, NaN)
    for i in eachindex(src)
      if !isnan(src[i])
        acc = isnan(acc) ? src[i] : min(acc, src[i])
      end
      tgt[i] = acc
    end
  end
  result = similar(arr)
  map_array_preserve_shape!(update!, result.a.data, arr.a.data, dims...;rev=rev)
  result
end

Base.cummax{T}(arr::AbstractArrayWrapper{Nullable{T}}, dims::Integer...;rev=false) = begin
  function update!(tgt, src)
    acc = Nullable{T}()
    for i in eachindex(src)
      if !src[i].isnull
        acc = acc.isnull ? src[i] : Nullable(max(acc.value, src[i].value))
      end
      tgt[i] = acc
    end
  end
  result = similar(arr)
  map_array_preserve_shape!(update!, result.a, arr.a, dims...;rev=rev)
  result
end
Base.cummax{T,N,A}(arr::AbstractArrayWrapper{Nullable{T},N,FloatNAArray{T,N,A}},
                   dims::Integer...
                   ;rev=false) = begin
  function update!(tgt, src)
    acc = convert(T, NaN)
    for i in eachindex(src)
      if !isnan(src[i])
        acc = isnan(acc) ? src[i] : max(acc, src[i])
      end
      tgt[i] = acc
    end
  end
  result = similar(arr)
  map_array_preserve_shape!(update!, result.a.data, arr.a.data, dims...;rev=rev)
  result
end

cummiddle{T}(arr::AbstractArrayWrapper{Nullable{T}}, dims::Integer...;rev=false) = begin
  DIVTYPE = typeof(one(T) / 2)
  function update!(tgt, src)
    accmin = Nullable{T}()
    accmax = Nullable{T}()
    for i in eachindex(src)
      if !src[i].isnull
        accmin = accmin.isnull ? src[i] : Nullable(min(accmin.value, src[i].value))
        accmax = accmax.isnull ? src[i] : Nullable(max(accmax.value, src[i].value))
      end
      tgt[i] = (accmin.isnull || accmax.isnull ? Nullable{DIVTYPE}() : Nullable((accmin.value + accmax.value) / 2)) :: Nullable{DIVTYPE}
    end
  end
  result = similar(arr, Nullable{DIVTYPE})
  map_array_preserve_shape!(update!, result.a, arr.a, dims...;rev=rev)
  result
end
cummiddle{T,N,A}(arr::AbstractArrayWrapper{Nullable{T},N,FloatNAArray{T,N,A}}, dims::Integer...;rev=false) = begin
  DIVTYPE = typeof(one(T) / 2)
  function update!(tgt, src)
    accmin = convert(T, NaN)
    accmax = convert(T, NaN)
    for i in eachindex(src)
      if !isnan(src[i])
        accmin = isnan(accmin) ? src[i] : min(accmin, src[i])
        accmax = isnan(accmax) ? src[i] : max(accmax, src[i])
      end
      tgt[i] = isnan(accmin) || isnan(accmax) ? convert(DIVTYPE, NaN) : (accmin + accmax) / 2
    end
  end
  result = similar(arr, Nullable{DIVTYPE})
  map_array_preserve_shape!(update!, result.a.data, arr.a.data, dims...;rev=rev)
  result
end

nafill0(arr::AbstractArrayWrapper, dims::Integer...;rev=false) = AbstractArrayWrapper(nafill0(arr.a, dims...; rev=rev))
nafill0{T}(arr::AbstractArray{Nullable{T}}, dims::Integer...;rev=false) = begin
  function update!(tgt, src)
    lastelem = Nullable{T}()
    for i in eachindex(src)
      if !src[i].isnull
        lastelem = src[i]
      end
      tgt[i] = lastelem
    end
  end
  result = similar(arr)
  map_array_preserve_shape!(update!, result, arr, dims...;rev=rev)
  result
end
nafill0{T}(arr::FloatNAArray{T}, dims::Integer...;rev=false) = begin
  function update!(tgt, src)
    lastelem = convert(T, NaN)
    for i in eachindex(src)
      if !isnan(src[i])
        lastelem = src[i]
      end
      tgt[i] = lastelem
    end
  end
  result = similar(arr)
  map_array_preserve_shape!(update!, result.data, arr.data, dims...;rev=rev)
  result
end

"""

`nafill(arr, dims... [; rev=false, window=0])` for `arr` of type `AbstractArrayWrapper`/`LabeledArray`/`DictArray`.

Fill forward (backward if `rev=true`) `arr` using non-null values from the last `window` elements, or latest non-null value from the beginning if `window=0`.

##### Arguments

* `arr`: `AbstractArrayWrapper`/`LabeledArray`/`DictArray`, the input array. When applied to `DictArray`, `nafill` is applied to each field. When applied to `LabeledArray`, `nafill` is applied to the base.
* `dims`: by default `dims=(1,)`. That is, the fill forward is performed along the first direction. If `dims=(n1, n2,...)`, for each slice spanned along the directions `n1`, `n2`, ..., the fill forward is taken along the leading dimension in `dims` first (i.e. `sum(dims)`), and then the next dimension, and so on.
* `rev`: If `rev=true`, the backward filling is calculated instead, starting for the last elements. By default, `rev=false`.
* `window`: If `window>0`, only the last `window` elements, including the one in consideration, will be used to fill forward. If `window=0`, `nafill` fills forward `arr` using all the elements so far. `NA` will be ignored.

##### Examples

```julia
julia> t = @nalift([1 NA;NA 4;NA NA])
3x2 DataCubes.AbstractArrayWrapper{Nullable{Int64},2,Array{Nullable{Int64},2}}:
 Nullable(1)        Nullable{Int64}()
 Nullable{Int64}()  Nullable(4)
 Nullable{Int64}()  Nullable{Int64}()

julia> nafill(t)
3x2 DataCubes.AbstractArrayWrapper{Nullable{Int64},2,Array{Nullable{Int64},2}}:
 Nullable(1)  Nullable{Int64}()
 Nullable(1)  Nullable(4)
 Nullable(1)  Nullable(4)

julia> nafill(t,2)
3x2 DataCubes.AbstractArrayWrapper{Nullable{Int64},2,Array{Nullable{Int64},2}}:
 Nullable(1)        Nullable(1)
 Nullable{Int64}()  Nullable(4)
 Nullable{Int64}()  Nullable{Int64}()

julia> nafill(t,2,1)
3x2 DataCubes.AbstractArrayWrapper{Nullable{Int64},2,Array{Nullable{Int64},2}}:
 Nullable(1)  Nullable(1)
 Nullable(1)  Nullable(4)
 Nullable(1)  Nullable(4)

julia> nafill(t, rev=true)
3x2 DataCubes.AbstractArrayWrapper{Nullable{Int64},2,Array{Nullable{Int64},2}}:
 Nullable(1)        Nullable(4)
 Nullable{Int64}()  Nullable(4)
 Nullable{Int64}()  Nullable{Int64}()

julia> nafill(t, window=2)
3x2 DataCubes.AbstractArrayWrapper{Nullable{Int64},2,Array{Nullable{Int64},2}}:
 Nullable(1)        Nullable{Int64}()
 Nullable(1)        Nullable(4)
 Nullable{Int64}()  Nullable(4)
```

"""
function nafill end

nafill(arr::AbstractArrayWrapper, dims::Integer...;rev=false, window=0) = AbstractArrayWrapper(nafill(arr.a, dims...; rev=rev, window=window))
nafill{T}(arr::AbstractArray{Nullable{T}}, dims::Integer...;rev=false, window=0) = begin
  if window == 0
    return nafill0(arr, dims...;rev=rev)
  end
  function update!(tgt, src)
    lastelem = Nullable{T}()
    age = 0
    for i in eachindex(src)
      if src[i].isnull
        age += 1
        if age >= window
          lastelem = Nullable{T}()
        end
      else
        lastelem = src[i]
        age = 0
      end
      tgt[i] = lastelem
    end
  end
  result = similar(arr)
  map_array_preserve_shape!(update!, result, arr, dims...;rev=rev)
  result
end
nafill{T}(arr::FloatNAArray{T}, dims::Integer...;rev=false, window=0) = begin
  if window == 0
    return nafill0(arr, dims...;rev=rev)
  end
  function update!(tgt, src)
    lastelem = convert(T, NaN)
    age = 0
    for i in eachindex(src)
      if isnan(src[i])
        age += 1
        if age >= window
          lastelem = convert(T, NaN)
        end
      else
        lastelem = src[i]
        age = 0
      end
      tgt[i] = lastelem
    end
  end
  result = similar(arr)
  map_array_preserve_shape!(update!, result.data, arr.data, dims...;rev=rev)
  result
end

"""

`diff(arr, dims... [; rev=false])` for `arr` of type `AbstractArrayWrapper`/`LabeledArray`/`DictArray`.

Take the difference between adjacent elements of `arr` along the directions belonging to the integers `dims`.
Note that `diff` applied to `AbstractArrayWrapper` (or to `LabeledArray` or `DictArray` by extension) will have the same shape as the original array. The first elements will be the first elements of the input array. This will ensure cumsum(diff(arr)) == diff(cumsum(arr)) == arr if there is no `Nullable` element.

##### Arguments

* `arr`: `AbstractArrayWrapper`/`LabeledArray`/`DictArray`, the input array. When applied to `DictArray`, `diff` is applied to each field. When applied to `LabeledArray`, `diff` is applied to the base.
* `dims`: by default `dims=(1,)`. That is, difference is calculated along the first direction. If `dims=(n1, n2,...)`, for each slice spanned along the directions `n1`, `n2`, ..., difference is taken along the leading dimension in `dims` first (i.e. `sum(dims)`), and then the next dimension, and so on.
* `rev`: If `rev=true`, difference is taken backward starting for the last elements. By default, `rev=false`.

##### Examples

```julia
julia> diff(@nalift([10,NA,12,14,17]))
5-element DataCubes.AbstractArrayWrapper{Nullable{Int64},1,Array{Nullable{Int64},1}}:
 Nullable(10)
 Nullable{Int64}()
 Nullable{Int64}()
 Nullable(2)
 Nullable(3)

julia> diff(darr(a=[11 12 13;14 15 16], b=[10 9 8;7 6 5]), 1, 2)
2 x 3 DictArray

a  b  |a  b  |a  b
------+------+------
11 10 |-2 2  |-2 2
3  -3 |3  -3 |3  -3

julia> diff(darr(a=[11 12 13;14 15 16], b=[10 9 8;7 6 5]), 1, 2, rev=true)
2 x 3 DictArray

a  b  |a  b  |a  b
------+------+-----
-3 3  |-3 3  |-3 3
2  -2 |2  -2 |16 5

julia> diff(larr(a=[11 12 13;14 15 16], b=[10 9 8;7 6 5]), 1, 2, rev=true)
2 x 3 LabeledArray

  |1     |2     |3
--+------+------+-----
  |a  b  |a  b  |a  b
--+------+------+-----
1 |-3 3  |-3 3  |-3 3
2 |2  -2 |2  -2 |16 5
```

"""
function diff end

Base.diff{T}(arr::AbstractArrayWrapper{Nullable{T}}, dims::Integer...;rev=false) = begin
  function update!(tgt, src)
    acc = Nullable(zero(T))
    for i in eachindex(src)
      tgt[i] = naop_minus(src[i], acc)
      acc = src[i]
    end
  end
  result = similar(arr)
  map_array_preserve_shape!(update!, result.a, arr.a, dims...;rev=rev)
  result
end
Base.diff{T,N,A}(arr::AbstractArrayWrapper{Nullable{T},N,FloatNAArray{T,N,A}}, dims::Integer...;rev=false) = begin
  function update!(tgt, src)
    acc = zero(T)
    for i in eachindex(src)
      tgt[i] = src[i] - acc
      acc = src[i]
    end
  end
  result = similar(arr)
  map_array_preserve_shape!(update!, result.a.data, arr.a.data, dims...;rev=rev)
  result
end

"""

`msum(arr, dims... [; rev=false, window=0])` for `arr` of type `AbstractArrayWrapper`/`LabeledArray`/`DictArray`.

Calculate moving sum of `arr` using the last `window` elements, or cumulative sum if `window=0`.

##### Arguments

* `arr`: `AbstractArrayWrapper`/`LabeledArray`/`DictArray`, the input array. When applied to `DictArray`, `msum` is applied to each field. When applied to `LabeledArray`, `msum` is applied to the base.
* `dims`: by default `dims=(1,)`. That is, moving sum is performed in the first direction. If `dims=(n1, n2,...)`, for each slice spanned along the directions `n1`, `n2`, ..., moving sum is taken along the leading dimension in `dims` first (i.e. `sum(dims)`), and then the next dimension, and so on.
* `rev`: If `rev=true`, moving sum is calculated backward starting for the last elements. By default, `rev=false`.
* `window`: If `window>0`, only the last `window` elements, including the one in consideration, will be used to calculate moving sum. If `window=0`, `sum` calculates the cumulative sum. `NA` will be ignored.

##### Examples

```julia
julia> msum(@nalift([10,11,12,14,17]))
5-element DataCubes.AbstractArrayWrapper{Nullable{Int64},1,Array{Nullable{Int64},1}}:
 Nullable(10)
 Nullable(21)
 Nullable(33)
 Nullable(47)
 Nullable(64)

julia> msum(@nalift([10,NA,12,14,17]))
5-element DataCubes.AbstractArrayWrapper{Nullable{Int64},1,Array{Nullable{Int64},1}}:
 Nullable(10)
 Nullable(10)
 Nullable(22)
 Nullable(36)
 Nullable(53)

julia> msum(darr(a=[11 12 13;14 15 16], b=[10 9 8;7 6 5]), 1, 2)
2 x 3 DictArray

a  b  |a  b  |a  b
------+------+------
11 10 |37 26 |65 40
25 17 |52 32 |81 45

julia> msum(larr(a=[11 12 13;14 15 16], b=[10 9 8;7 6 5]), 2, 1, rev=true)
2 x 3 LabeledArray

  |1     |2     |3
--+------+------+------
  |a  b  |a  b  |a  b
--+------+------+------
1 |81 45 |56 28 |29 13
2 |70 35 |44 19 |16 5
```

"""
function msum end

msum{T}(arr::AbstractArrayWrapper{Nullable{T}}, dims::Integer...;rev=false,window=0) = begin
  if window == 0
    return cumsum(arr, dims...;rev=rev)
  end

  function update!(tgt, src)
    ringbuf = fill(zero(T), window)
    acc = zero(T)
    ringbuf_index = 1
    for i in eachindex(src)
      if src[i].isnull
        acc -= ringbuf[ringbuf_index]
        ringbuf[ringbuf_index] = zero(T)
      else
        acc += src[i].value - ringbuf[ringbuf_index]
        ringbuf[ringbuf_index] = src[i].value
      end
      if ringbuf_index == window
        ringbuf_index = 1
      else
        ringbuf_index += 1
      end
      tgt[i] = Nullable(acc)
    end
  end
  result = similar(arr)
  map_array_preserve_shape!(update!, result.a, arr.a, dims...;rev=rev)
  result
end
msum{T,N,A}(arr::AbstractArrayWrapper{Nullable{T},N,FloatNAArray{T,N,A}}, dims::Integer...;rev=false,window=0) = begin
  if window == 0
    return cumsum(arr, dims...;rev=rev)
  end

  function update!(tgt, src)
    ringbuf = fill(zero(T), window)
    acc = zero(T)
    ringbuf_index = 1
    for i in eachindex(src)
      if isnan(src[i])
        acc -= ringbuf[ringbuf_index]
        ringbuf[ringbuf_index] = zero(T)
      else
        acc += src[i] - ringbuf[ringbuf_index]
        ringbuf[ringbuf_index] = src[i]
      end
      if ringbuf_index == window
        ringbuf_index = 1
      else
        ringbuf_index += 1
      end
      tgt[i] = acc
    end
  end
  result = similar(arr)
  map_array_preserve_shape!(update!, result.a.data, arr.a.data, dims...;rev=rev)
  result
end

"""

`mprod(arr, dims... [; rev=false, window=0])` for `arr` of type `AbstractArrayWrapper`/`LabeledArray`/`DictArray`.

Calculate moving product of `arr` using the last `window` elements, or cumulative product if `window=0`.

##### Arguments

* `arr`: `AbstractArrayWrapper`/`LabeledArray`/`DictArray`, the input array. When applied to `DictArray`, `mprod` is applied to each field. When applied to `LabeledArray`, `mprod` is applied to the base.
* `dims`: by default `dims=(1,)`. That is, moving product is performed in the first direction. If `dims=(n1, n2,...)`, for each slice spanned along the directions `n1`, `n2`, ..., moving product is taken along the leading dimension in `dims` first (i.e. `prod(dims)`), and then the next dimension, and so on.
* `rev`: If `rev=true`, moving product is calculated backward starting for the last elements. By default, `rev=false`.
* `window`: If `window>0`, only the last `window` elements, including the one in consideration, will be used to calculate moving product. If `window=0`, `prod` calculates the cumulative product. `NA` will be ignored.

##### Examples

```julia
julia> mprod(@nalift([10,11,12,14,17]))
5-element DataCubes.AbstractArrayWrapper{Nullable{Int64},1,Array{Nullable{Int64},1}}:
 Nullable(10)
 Nullable(110)
 Nullable(1320)
 Nullable(18480)
 Nullable(314160)

julia> mprod(@nalift([10,NA,12,14,17]))
5-element DataCubes.AbstractArrayWrapper{Nullable{Int64},1,Array{Nullable{Int64},1}}:
 Nullable(10)
 Nullable(10)
 Nullable(120)
 Nullable(1680)
 Nullable(28560)

julia> mprod(darr(a=[11 12 13;14 15 16], b=[10 9 8;7 6 5]), 1, 2)
2 x 3 DictArray

a   b  |a     b    |a       b
-------+-----------+---------------
11  10 |1848  630  |360360  30240
154 70 |27720 3780 |5765760 151200

julia> mprod(larr(a=[11 12 13;14 15 16], b=[10 9 8;7 6 5]), 2, 1, rev=true)
2 x 3 LabeledArray

  |1              |2          |3
--+---------------+-----------+-------
  |a       b      |a     b    |a   b
--+---------------+-----------+-------
1 |5765760 151200 |37440 2160 |208 40
2 |524160  15120  |3120  240  |16  5
```

"""
function mprod end

mprod{T}(arr::AbstractArrayWrapper{Nullable{T}}, dims::Integer...;rev=false, window=0) = begin
  if window == 0
    return cumprod(arr, dims...;rev=rev)
  end
  function update!(tgt, src)
    ringbuf = fill(one(T), window)
    acc = one(T)
    ringbuf_index = 1
    for i in eachindex(src)
      if src[i].isnull
        acc /= ringbuf[ringbuf_index]
        ringbuf[ringbuf_index] = one(T)
      else
        acc *= src[i].value / ringbuf[ringbuf_index]
        ringbuf[ringbuf_index] = src[i].value
      end
      if ringbuf_index == window
        ringbuf_index = 1
      else
        ringbuf_index += 1
      end
      tgt[i] = Nullable(acc)
    end
  end
  result = similar(arr)
  map_array_preserve_shape!(update!, result.a, arr.a, dims...;rev=rev)
  result
end
mprod{T,N,A}(arr::AbstractArrayWrapper{Nullable{T},N,FloatNAArray{T,N,A}}, dims::Integer...;rev=false, window=0) = begin
  if window == 0
    return cumprod(arr, dims...;rev=rev)
  end
  function update!(tgt, src)
    ringbuf = fill(one(T), window)
    acc = one(T)
    ringbuf_index = 1
    for i in eachindex(src)
      if isnan(src[i])
        acc /= ringbuf[ringbuf_index]
        ringbuf[ringbuf_index] = one(T)
      else
        acc *= src[i] / ringbuf[ringbuf_index]
        ringbuf[ringbuf_index] = src[i]
      end
      if ringbuf_index == window
        ringbuf_index = 1
      else
        ringbuf_index += 1
      end
      tgt[i] = acc
    end
  end
  result = similar(arr)
  map_array_preserve_shape!(update!, result.a.data, arr.a.data, dims...;rev=rev)
  result
end

"""

`mmean(arr, dims... [; rev=false, window=0])` for `arr` of type `AbstractArrayWrapper`/`LabeledArray`/`DictArray`.

Calculate moving mean of `arr` using the last `window` elements, or cumulative mean if `window=0`.

##### Arguments

* `arr`: `AbstractArrayWrapper`/`LabeledArray`/`DictArray`, the input array. When applied to `DictArray`, `mmean` is applied to each field. When applied to `LabeledArray`, `mmean` is applied to the base.
* `dims`: by default `dims=(1,)`. That is, moving mean is performed in the first direction. If `dims=(n1, n2,...)`, for each slice spanned along the directions `n1`, `n2`, ..., moving mean is taken along the leading dimension in `dims` first (i.e. `mean(dims)`), and then the next dimension, and so on.
* `rev`: If `rev=true`, moving mean is calculated backward starting for the last elements. By default, `rev=false`.
* `window`: If `window>0`, only the last `window` elements, including the one in consideration, will be used to calculate moving mean. If `window=0`, `mean` calculates the cumulative mean. `NA` will be ignored.

##### Examples

```julia
julia> mmean(@nalift([10,11,12,14,17]))
5-element DataCubes.AbstractArrayWrapper{Nullable{Float64},1,Array{Nullable{Float64},1}}:
 Nullable(10.0)
 Nullable(10.5)
 Nullable(11.0)
 Nullable(11.75)
 Nullable(12.8)

julia> mmean(@nalift([10,NA,12,14,17]))
5-element DataCubes.AbstractArrayWrapper{Nullable{Float64},1,Array{Nullable{Float64},1}}:
 Nullable(10.0)
 Nullable(10.0)
 Nullable(11.0)
 Nullable(12.0)
 Nullable(13.25)

julia> mmean(darr(a=[11 12 13;14 15 16], b=[10 9 8;7 6 5]), 1, 2)
2 x 3 DictArray

a    b    |a                  b                 |a    b
----------+-------------------------------------+---------
11.0 10.0 |12.333333333333334 8.666666666666666 |13.0 8.0
12.5 8.5  |13.0               8.0               |13.5 7.5

julia> mmean(larr(a=[11 12 13;14 15 16], b=[10 9 8;7 6 5]), 2, 1, rev=true)
2 x 3 LabeledArray

  |1        |2                                    |3
--+---------+-------------------------------------+---------
  |a    b   |a                  b                 |a    b
--+---------+-------------------------------------+---------
1 |13.5 7.5 |14.0               7.0               |14.5 6.5
2 |14.0 7.0 |14.666666666666666 6.333333333333333 |16.0 5.0
```

"""
function mmean end

mmean{T}(arr::AbstractArrayWrapper{Nullable{T}}, dims::Integer...;rev=false, window=0) = begin
  DIVTYPE = typeof(one(T) / 1)
  result = similar(arr, Nullable{DIVTYPE})
  if window == 0
    map_array_preserve_shape!((tgt,src)->cummean_update!(tgt,src), result.a, arr.a, dims...;rev=rev)
  else
    map_array_preserve_shape!((tgt,src)->mmean_update!(tgt,src,window), result.a, arr.a, dims...;rev=rev)
  end
  result
end
mmean{T,N,A}(arr::AbstractArrayWrapper{Nullable{T},N,FloatNAArray{T,N,A}}, dims::Integer...;rev=false, window=0) = begin
  DIVTYPE = typeof(one(T) / 1)
  result = similar(arr, Nullable{DIVTYPE})
  if window == 0
    map_array_preserve_shape!((tgt,src)->cummean_update!(tgt,src), result.a.data, arr.a.data, dims...;rev=rev)
  else
    map_array_preserve_shape!((tgt,src)->mmean_update!(tgt,src,window), result.a.data, arr.a.data, dims...;rev=rev)
  end
  result
end

function mmean_update!{T,U}(tgt::AbstractArray{Nullable{T}}, src::AbstractArray{Nullable{U}}, window::Integer)
  ringbuf = fill(Nullable{U}(), window)
  acc = zero(U)
  num_non_nullable = 0
  ringbuf_index = 1
  for i in eachindex(src)
    srci = src[i]
    if !srci.isnull
      num_non_nullable += 1
      acc += srci.value
    end
    if !ringbuf[ringbuf_index].isnull
      num_non_nullable -= 1
      acc -= ringbuf[ringbuf_index].value
    end
    ringbuf[ringbuf_index] = srci
    if ringbuf_index == window
      ringbuf_index = 1
    else
      ringbuf_index += 1
    end
    tgt[i] = (num_non_nullable == 0 ? Nullable{T}() : Nullable(acc / num_non_nullable))::Nullable{T}
  end
end
function mmean_update!{T<:AbstractFloat,U<:AbstractFloat}(tgt::AbstractArray{T}, src::AbstractArray{U}, window::Integer)
  ringbuf = fill(convert(U,NaN), window)
  acc = zero(U)
  num_non_nullable = 0
  ringbuf_index = 1
  for i in eachindex(src)
    srci = src[i]
    if !isnan(srci)
      num_non_nullable += 1
      acc += srci
    end
    if !isnan(ringbuf[ringbuf_index])
      num_non_nullable -= 1
      acc -= ringbuf[ringbuf_index]
    end
    ringbuf[ringbuf_index] = srci
    if ringbuf_index == window
      ringbuf_index = 1
    else
      ringbuf_index += 1
    end
    tgt[i] = num_non_nullable == 0 ? convert(T,NaN) : acc / num_non_nullable
  end
end

function cummean_update!{T,U}(tgt::AbstractArray{Nullable{T}}, src::AbstractArray{Nullable{U}})
  acc = zero(U)
  num_non_nullable = 0
  for i in eachindex(src)
    srci = src[i]
    if !srci.isnull
      num_non_nullable += 1
      acc += srci.value
    end
    tgt[i] = (num_non_nullable == 0 ? Nullable{T}() : Nullable(acc / num_non_nullable))::Nullable{T}
  end
end
function cummean_update!{T<:AbstractFloat,U<:AbstractFloat}(tgt::AbstractArray{T}, src::AbstractArray{U})
  acc = zero(U)
  num_non_nullable = 0
  for i in eachindex(src)
    srci = src[i]
    if !isnan(srci)
      num_non_nullable += 1
      acc += srci
    end
    tgt[i] = num_non_nullable == 0 ? convert(T,NaN) : acc / num_non_nullable
  end
end

"""

`mminimum(arr, dims... [; rev=false, window=0])` for `arr` of type `AbstractArrayWrapper`/`LabeledArray`/`DictArray`.

Calculate moving minimum of `arr` using the last `window` elements, or cumulative minimum if `window=0`.

##### Arguments

* `arr`: `AbstractArrayWrapper`/`LabeledArray`/`DictArray`, the input array. When applied to `DictArray`, `mminimum` is applied to each field. When applied to `LabeledArray`, `mminimum` is applied to the base.
* `dims`: by default `dims=(1,)`. That is, moving minimum is performed in the first direction. If `dims=(n1, n2,...)`, for each slice spanned along the directions `n1`, `n2`, ..., moving minimum is taken along the leading dimension in `dims` first (i.e. `minimum(dims)`), and then the next dimension, and so on.
* `rev`: If `rev=true`, moving minimum is calculated backward starting for the last elements. By default, `rev=false`.
* `window`: If `window>0`, only the last `window` elements, including the one in consideration, will be used to calculate moving minimum. If `window=0`, `minimum` calculates the cumulative minimum. `NA` will be ignored.

##### Examples

```julia
julia> mminimum(@nalift([15,10,12,11,17]))
5-element DataCubes.AbstractArrayWrapper{Nullable{Int64},1,Array{Nullable{Int64},1}}:
 Nullable(15)
 Nullable(10)
 Nullable(10)
 Nullable(10)
 Nullable(10)

julia> mminimum(@nalift([15,NA,12,11,17]))
5-element DataCubes.AbstractArrayWrapper{Nullable{Int64},1,Array{Nullable{Int64},1}}:
 Nullable(15)
 Nullable(15)
 Nullable(12)
 Nullable(11)
 Nullable(11)

julia> mminimum(darr(a=[11 12 13;14 15 16], b=[10 9 8;7 6 5]), 1, 2)
2 x 3 DictArray

a  b  |a  b |a  b
------+-----+-----
11 10 |11 7 |11 6
11 7  |11 6 |11 5

julia> mminimum(larr(a=[11 12 13;14 15 16], b=[10 9 8;7 6 5]), 2, 1, rev=true)
2 x 3 LabeledArray

  |1    |2    |3
--+-----+-----+-----
  |a  b |a  b |a  b
--+-----+-----+-----
1 |11 5 |12 5 |13 5
2 |12 5 |13 5 |16 5
```

"""
function mminimum end

mminimum{T}(arr::AbstractArrayWrapper{Nullable{T}}, dims::Integer...;rev=false, window=0) = begin
  if window == 0
    return cummin(arr, dims...;rev=rev)
  end
  result = similar(arr, Nullable{T})
  map_array_preserve_shape!((tgt,src) -> moving_update!(minimum, window, tgt, src), result, arr, dims...;rev=rev)
  result
end

"""

`mmaximum(arr, dims... [; rev=false, window=0])` for `arr` of type `AbstractArrayWrapper`/`LabeledArray`/`DictArray`.

Calculate moving maximum of `arr` using the last `window` elements, or cumulative maximum if `window=0`.

##### Arguments

* `arr`: `AbstractArrayWrapper`/`LabeledArray`/`DictArray`, the input array. When applied to `DictArray`, `mmaximum` is applied to each field. When applied to `LabeledArray`, `mmaximum` is applied to the base.
* `dims`: by default `dims=(1,)`. That is, moving maximum is performed in the first direction. If `dims=(n1, n2,...)`, for each slice spanned along the directions `n1`, `n2`, ..., moving maximum is taken along the leading dimension in `dims` first (i.e. `minimum(dims)`), and then the next dimension, and so on.
* `rev`: If `rev=true`, moving maximum is calculated backward starting for the last elements. By default, `rev=false`.
* `window`: If `window>0`, only the last `window` elements, including the one in consideration, will be used to calculate moving maximum. If `window=0`, `maximum` calculates the cumulative maximum. `NA` will be ignored.

##### Examples

```julia
julia> mmaximum(@nalift([11,14,12,11,17]))
5-element DataCubes.AbstractArrayWrapper{Nullable{Int64},1,Array{Nullable{Int64},1}}:
 Nullable(11)
 Nullable(14)
 Nullable(14)
 Nullable(14)
 Nullable(17)

julia> mmaximum(@nalift([11,NA,12,11,17]))
5-element DataCubes.AbstractArrayWrapper{Nullable{Int64},1,Array{Nullable{Int64},1}}:
 Nullable(11)
 Nullable(11)
 Nullable(12)
 Nullable(12)
 Nullable(17)

julia> mmaximum(darr(a=[11 12 13;14 15 16], b=[10 9 8;7 6 5]), 1, 2)
2 x 3 DictArray

a  b  |a  b  |a  b
------+------+------
11 10 |14 10 |15 10
14 10 |15 10 |16 10

julia> mmaximum(larr(a=[11 12 13;14 15 16], b=[10 9 8;7 6 5]), 2, 1, rev=true)
2 x 3 LabeledArray

  |1     |2    |3
--+------+-----+-----
  |a  b  |a  b |a  b
--+------+-----+-----
1 |16 10 |16 9 |16 8
2 |16 9  |16 8 |16 5
```

"""
function mmaximum end

mmaximum{T}(arr::AbstractArrayWrapper{Nullable{T}}, dims::Integer...;rev=false, window=0) = begin
  if window == 0
    return cummax(arr, dims...;rev=rev)
  end
  result = similar(arr, Nullable{T})
  map_array_preserve_shape!((tgt,src) -> moving_update!(maximum, window, tgt, src), result, arr, dims...;rev=rev)
  result
end

"""

`mmedian(arr, dims... [; rev=false, window=0])` for `arr` of type `AbstractArrayWrapper`/`LabeledArray`/`DictArray`.

Calculate moving median of `arr` using the last `window` elements, or cumulative median if `window=0`.

##### Arguments

* `arr`: `AbstractArrayWrapper`/`LabeledArray`/`DictArray`, the input array. When applied to `DictArray`, `mmedian` is applied to each field. When applied to `LabeledArray`, `mmedian` is applied to the base.
* `dims`: by default `dims=(1,)`. That is, moving median is performed in the first direction. If `dims=(n1, n2,...)`, for each slice spanned along the directions `n1`, `n2`, ..., moving median is taken along the leading dimension in `dims` first (i.e. `minimum(dims)`), and then the next dimension, and so on.
* `rev`: If `rev=true`, moving median is calculated backward starting for the last elements. By default, `rev=false`.
* `window`: If `window>0`, only the last `window` elements, including the one in consideration, will be used to calculate moving median. If `window=0`, `median` calculates the cumulative median. `NA` will be ignored.

##### Examples

```julia
julia> mmedian(@nalift([11,14,12,11,17]))
5-element DataCubes.AbstractArrayWrapper{Nullable{Float64},1,Array{Nullable{Float64},1}}:
 Nullable(11.0)
 Nullable(12.5)
 Nullable(12.5)
 Nullable(12.5)
 Nullable(14.0)

julia> mmedian(@nalift([11,NA,12,11,17]))
5-element DataCubes.AbstractArrayWrapper{Nullable{Float64},1,Array{Nullable{Float64},1}}:
 Nullable(11.0)
 Nullable(11.0)
 Nullable(11.5)
 Nullable(11.5)
 Nullable(12.0)

julia> mmedian(darr(a=[11 12 13;14 15 16], b=[10 9 8;7 6 5]), 1, 2)
2 x 3 DictArray

a    b    |a    b   |a    b
----------+---------+---------
11.0 10.0 |12.5 8.5 |14.0 8.5
12.5 8.5  |14.0 8.5 |14.5 8.5

julia> mmedian(larr(a=[11 12 13;14 15 16], b=[10 9 8;7 6 5]), 2, 1, rev=true)
2 x 3 LabeledArray

  |1        |2        |3
--+---------+---------+---------
  |a    b   |a    b   |a    b
--+---------+---------+---------
1 |14.5 8.5 |14.5 8.0 |14.5 6.5
2 |14.5 8.0 |14.5 6.5 |16.0 5.0
```

"""
function mmedian end

mmedian{T}(arr::AbstractArrayWrapper{Nullable{T}}, dims::Integer...;rev=false, window=0) = begin
  if window == 0
    return cumquantile(arr, 0.5, dims...;rev=rev)
  end
  DIVTYPE = typeof(one(T) / 1)
  result = similar(arr, Nullable{DIVTYPE})
  map_array_preserve_shape!((tgt,src) -> moving_update!(median, window, tgt, src), result, arr, dims...;rev=rev)
  result
end

"""

`mmiddle(arr, dims... [; rev=false, window=0])` for `arr` of type `AbstractArrayWrapper`/`LabeledArray`/`DictArray`.

Calculate moving middle of `arr` using the last `window` elements, or cumulative middle if `window=0`.

##### Arguments

* `arr`: `AbstractArrayWrapper`/`LabeledArray`/`DictArray`, the input array. When applied to `DictArray`, `mmiddle` is applied to each field. When applied to `LabeledArray`, `mmiddle` is applied to the base.
* `dims`: by default `dims=(1,)`. That is, moving middle is performed in the first direction. If `dims=(n1, n2,...)`, for each slice spanned along the directions `n1`, `n2`, ..., moving middle is taken along the leading dimension in `dims` first (i.e. `minimum(dims)`), and then the next dimension, and so on.
* `rev`: If `rev=true`, moving middle is calculated backward starting for the last elements. By default, `rev=false`.
* `window`: If `window>0`, only the last `window` elements, including the one in consideration, will be used to calculate moving middle. If `window=0`, `middle` calculates the cumulative middle. `NA` will be ignored.

##### Examples

```julia
julia> mmiddle(@nalift([11,14,12,11,17]))
5-element DataCubes.AbstractArrayWrapper{Nullable{Float64},1,Array{Nullable{Float64},1}}:
 Nullable(11.0)
 Nullable(12.5)
 Nullable(12.5)
 Nullable(12.5)
 Nullable(14.0)

julia> mmiddle(@nalift([11,NA,12,11,17]))
5-element DataCubes.AbstractArrayWrapper{Nullable{Float64},1,Array{Nullable{Float64},1}}:
 Nullable(11.0)
 Nullable(11.0)
 Nullable(11.5)
 Nullable(11.5)
 Nullable(14.0)

julia> mmiddle(darr(a=[11 12 13;14 15 16], b=[10 9 8;7 6 5]), 1, 2)
2 x 3 DictArray

a    b    |a    b   |a    b
----------+---------+---------
11.0 10.0 |12.5 8.5 |13.0 8.0
12.5 8.5  |13.0 8.0 |13.5 7.5

julia> mmiddle(larr(a=[11 12 13;14 15 16], b=[10 9 8;7 6 5]), 2, 1, rev=true)
2 x 3 LabeledArray

  |1        |2        |3
--+---------+---------+---------
  |a    b   |a    b   |a    b
--+---------+---------+---------
1 |13.5 7.5 |14.0 7.0 |14.5 6.5
2 |14.0 7.0 |14.5 6.5 |16.0 5.0
```

"""
function mmiddle end

mmiddle{T}(arr::AbstractArrayWrapper{Nullable{T}}, dims::Integer...;rev=false, window=0) = begin
  if window == 0
    return cummiddle(arr, dims...;rev=rev)
  end
  DIVTYPE = typeof(one(T) / 1)
  result = similar(arr, Nullable{DIVTYPE})
  map_array_preserve_shape!((tgt,src) -> moving_update!(middle, window, tgt, src), result, arr, dims...;rev=rev)
  result
end

"""

`mquantile(arr, quantile, dims... [; rev=false, window=0])` for `arr` of type `AbstractArrayWrapper`/`LabeledArray`/`DictArray`.

Calculate moving quantile of `arr` using the last `window` elements, or cumulative quantile if `window=0`.

##### Arguments

* `arr`: `AbstractArrayWrapper`/`LabeledArray`/`DictArray`, the input array. When applied to `DictArray`, `mquantile` is applied to each field. When applied to `LabeledArray`, `mquantile` is applied to the base.
* `quantile`: a number between 0 and 1 for the quantile to calculate.
* `dims`: by default `dims=(1,)`. That is, moving quantile is performed in the first direction. If `dims=(n1, n2,...)`, for each slice spanned along the directions `n1`, `n2`, ..., moving quantile is taken along the leading dimension in `dims` first (i.e. `minimum(dims)`), and then the next dimension, and so on.
* `rev`: If `rev=true`, moving quantile is calculated backward starting for the last elements. By default, `rev=false`.
* `window`: If `window>0`, only the last `window` elements, including the one in consideration, will be used to calculate moving quantile. If `window=0`, `quantile` calculates the cumulative quantile. `NA` will be ignored.

##### Examples

```julia
julia> mquantile(@nalift([11,14,12,11,17]), 0.25)
5-element DataCubes.AbstractArrayWrapper{Nullable{Float64},1,Array{Nullable{Float64},1}}:
 Nullable(11.0)
 Nullable(11.75)
 Nullable(11.75)
 Nullable(11.75)
 Nullable(12.5)

julia> mquantile(@nalift([11,NA,12,11,17]), 0.25)
5-element DataCubes.AbstractArrayWrapper{Nullable{Float64},1,Array{Nullable{Float64},1}}:
 Nullable(11.0)
 Nullable(11.0)
 Nullable(11.25)
 Nullable(11.25)
 Nullable(11.5)

julia> mquantile(darr(a=[11 12 13;14 15 16], b=[10 9 8;7 6 5]), 0.25, 1, 2)
2 x 3 DictArray

a     b    |a     b    |a     b
-----------+-----------+-----------
11.0  10.0 |11.75 7.75 |12.5  7.75
11.75 7.75 |12.5  7.75 |13.25 7.75

julia> mquantile(larr(a=[11 12 13;14 15 16], b=[10 9 8;7 6 5]), 0.25, 2, 1, rev=true)
2 x 3 LabeledArray

  |1          |2          |3
--+-----------+-----------+-----------
  |a     b    |a     b    |a     b
--+-----------+-----------+-----------
1 |13.75 7.25 |13.75 6.5  |13.75 5.75
2 |13.75 6.5  |13.75 5.75 |16.0  5.0
```

"""
function mquantile end

mquantile{T}(arr::AbstractArrayWrapper{Nullable{T}}, q::Number, dims::Integer...;rev=false, window=0) = begin
  if window == 0
    return cumquantile(arr, q, dims...;rev=rev)
  end
  DIVTYPE = typeof(one(T) / 1)
  result = similar(arr, Nullable{DIVTYPE})
  # we don't mind changing the input src. native quantile function was not fast enough.
  map_array_preserve_shape!((tgt,src) -> moving_update!(x->mquantile_quickselect!(DIVTYPE, x, q), window, tgt, src), result, arr, dims...;rev=rev)
  result
end

mquantile_quickselect!{T,R}(::Type{T}, arr::AbstractArray, p::R) =
  mquantile_quickselect_inner!(T, arr, 1,
                               length(arr),
                               one(R) + p * (length(arr)-1))
mquantile_quickselect_inner!{T,R}(::Type{T}, arr::AbstractArray, left::Int, right::Int, expected_position::R) = begin
  while true
    if left == right
      return convert(T, arr[left])
    elseif left == right - 1
      return convert(T, (expected_position - left) * arr[right] + (right - expected_position) * arr[left])
    else
      pivotindex = div(right + left, 2)
      pivotindex = mquantile_partition(arr, left, right, pivotindex)
      if expected_position == pivotindex
        return convert(T, arr[pivotindex])
      elseif expected_position < pivotindex
        right = pivotindex
      else
        left = pivotindex
      end
    end
  end
  error("this should not be reacheable.")
end
mquantile_partition(arr::AbstractArray, left::Int, right::Int, pivotindex::Int) = begin
  pivotvalue = arr[pivotindex]
  temp = arr[right]; arr[right] = arr[pivotindex]; arr[pivotindex] = temp
  storeindex = left
  for i in left:right-1
    if arr[i] < pivotvalue
      temp = arr[storeindex]; arr[storeindex] = arr[i]; arr[i] = temp
      storeindex += 1
    end
  end
  temp = arr[right]; arr[right] = arr[storeindex]; arr[storeindex] = temp
  return storeindex
end
moving_update!{T,U}(f::Function, window::Integer, tgt::AbstractArrayWrapper{Nullable{T}}, src::AbstractArrayWrapper{Nullable{U}}) = moving_update!(f, window, tgt.a, src.a)
moving_update!{T,U}(f::Function, window::Integer, tgt::AbstractArrayWrapper{Nullable{T}}, src::AbstractArray{Nullable{U}}) = moving_update!(f, window, tgt.a, src)
moving_update!{T,U}(f::Function, window::Integer, tgt::AbstractArrayWrapper{Nullable{T}}, src::FloatNAArray{U}) = moving_update!(f, window, tgt.a, src)
moving_update!{T,U}(f::Function, window::Integer, tgt::AbstractArray{Nullable{T}}, src::AbstractArrayWrapper{Nullable{U}}) = moving_update!(f, window, tgt, src.a)
moving_update!{T,U}(f::Function, window::Integer, tgt::FloatNAArray{T}, src::AbstractArrayWrapper{Nullable{U}}) = moving_update!(f, window, tgt, src.a)

# generic case.
function moving_update!{T,U}(f::Function,
                             window::Integer,
                             tgt::AbstractArray{Nullable{T}},
                             src::AbstractArray{Nullable{U}})
  ringbuf = Array(Nullable{U}, window)
  projbuf = Array(U, window)
  ringbuf_index = 1
  eff_size = 1
  for i in eachindex(src)
    ringbuf[ringbuf_index] = src[i]
    index = 0
    for j in 1:eff_size
      if !ringbuf[j].isnull
        index += 1
        projbuf[index] = ringbuf[j].value
      end
    end
    tgt[i] = index==0 ? Nullable{T}() : Nullable(f(slice(projbuf, 1:index)))

    if ringbuf_index == window
      ringbuf_index = 1
    else
      ringbuf_index += 1
    end
    eff_size = min(window, eff_size+1)
  end
end

# tgt, src are float arrays.
function moving_update!{T,U}(f::Function,
                             window::Integer,
                             tgt::FloatNAArray{T},
                             src::FloatNAArray{U})
  ringbuf = Array(U, window)
  projbuf = Array(U, window)
  ringbuf_index = 1
  eff_size = 1
  srcdata = src.data
  na = convert(T, NaN)
  for i in eachindex(srcdata)
    ringbuf[ringbuf_index] = srcdata[i]
    index = 0
    for j in 1:eff_size
      if !isnan(ringbuf[j])
        index += 1
        projbuf[index] = ringbuf[j]
      end
    end
    tgt[i] = index==0 ? na : f(slice(projbuf, 1:index))

    if ringbuf_index == window
      ringbuf_index = 1
    else
      ringbuf_index += 1
    end
    eff_size = min(window, eff_size+1)
  end
end

# tgt is a float array.
function moving_update!{T,U}(f::Function,
                             window::Integer,
                             tgt::FloatNAArray{T},
                             src::AbstractArray{Nullable{U}})
  ringbuf = Array(Nullable{U}, window)
  projbuf = Array(U, window)
  ringbuf_index = 1
  eff_size = 1
  na = convert(T, NaN)
  for i in eachindex(src)
    ringbuf[ringbuf_index] = src[i]
    index = 0
    for j in 1:eff_size
      if !ringbuf[j].isnull
        index += 1
        projbuf[index] = ringbuf[j].value
      end
    end
    tgt[i] = index==0 ? na : f(slice(projbuf, 1:index))

    if ringbuf_index == window
      ringbuf_index = 1
    else
      ringbuf_index += 1
    end
    eff_size = min(window, eff_size+1)
  end
end

# src is a float array.
function moving_update!{T,U}(f::Function,
                             window::Integer,
                             tgt::AbstractArray{Nullable{T}},
                             src::FloatNAArray{U})
  ringbuf = Array(U, window)
  projbuf = Array(U, window)
  ringbuf_index = 1
  eff_size = 1
  srcdata = src.data
  for i in eachindex(srcdata)
    ringbuf[ringbuf_index] = srcdata[i]
    index = 0
    for j in 1:eff_size
      if !isnan(ringbuf[j])
        index += 1
        projbuf[index] = ringbuf[j]
      end
    end
    tgt[i] = index==0 ? Nullable{T}() : Nullable(f(slice(projbuf, 1:index)))

    if ringbuf_index == window
      ringbuf_index = 1
    else
      ringbuf_index += 1
    end
    eff_size = min(window, eff_size+1)
  end
end

cumquantile{T}(arr::AbstractArrayWrapper{Nullable{T}}, q::Number, dims::Integer...;rev=false) = begin
  DIVTYPE = typeof(one(T) / 1)
  result = similar(arr, Nullable{DIVTYPE})
  map_array_preserve_shape!((tgt,src) -> cumquantile_update!(tgt, src, q), result, arr, dims...;rev=rev)
  result
end

cumquantile_update!{T,U,R}(tgt::AbstractArray{Nullable{T}}, src::AbstractArray{Nullable{U}}, q::R) = begin
  lowqueue = Collections.PriorityQueue{U,U,Base.Order.Ordering}(Base.Order.Reverse) # queue for smaller elements. Easy to get the maximum element.
  highqueue = Collections.PriorityQueue{U,U,Base.Order.Ordering}(Base.Order.Forward) # queue for larger elements. Easy to get the minimum element.

  dummy_index = 0
  for i in eachindex(src)
    srci = src[i]
    if !srci.isnull
      dummy_index += 1
      if isempty(highqueue) || Collections.peek(highqueue)[2] < srci.value
        Collections.enqueue!(highqueue, dummy_index, srci.value)
      elseif isempty(lowqueue)
        Collections.enqueue!(lowqueue, dummy_index, srci.value)
      end
      # check the queues and shuffle elements if necessary.
      quantile_denom = length(lowqueue) + length(highqueue) - 1
      while q < (length(lowqueue)-1) / quantile_denom
        elem = Collections.peek(lowqueue)
        Collections.dequeue!(lowqueue)
        Collections.enqueue!(highqueue, elem...)
      end
      while q > length(lowqueue) / quantile_denom
        elem = Collections.peek(highqueue)
        Collections.dequeue!(highqueue)
        Collections.enqueue!(lowqueue, elem...)
      end
    end
    tgt[i] = if isempty(lowqueue)
      if isempty(highqueue)
        Nullable{T}()
      else
        Nullable(convert(T, Collections.peek(highqueue)[2]))
      end
    elseif isempty(highqueue)
      Nullable(convert(T, Collections.peek(lowqueue)[2]))
    else
      lowlen = length(lowqueue)
      highlen = length(highqueue)
      low_quantile = (lowlen - 1) / (lowlen + highlen - 1)
      high_quantile = lowlen / (lowlen + highlen - 1)
      lowpeek = Collections.peek(lowqueue)[2]
      highpeek = Collections.peek(highqueue)[2]
      Nullable((lowpeek + (q - low_quantile) * (highpeek - lowpeek) / (high_quantile - low_quantile))::T)
    end
  end
end

"""

`describe(arr)`

Generate a `LabeledArray` showing the overall statistics of the input.
If the input is a `Nullable` array, its summary statistics is calculated and the return value is of type `LDict`.
If the input is a `DictArray`, the summary is calculated for each field and the result is a `DictArray`.
If the input is a `LabeledArray`, `describe` returns the summary of its base.

##### Examples

```julia
julia> describe(@nalift([1,2,3,4,NA]))
DataCubes.LDict{Symbol,Any} with 10 entries:
  :min     => [Nullable(1)]
  :q1      => Nullable(1.75)
  :med     => Nullable(2.5)
  :q3      => Nullable(3.25)
  :max     => Nullable(4)
  :mean    => Nullable(2.5)
  :std     => Nullable(1.2909944487358056)
  :count   => Nullable(5)
  :nacount => Nullable(1)
  :naratio => Nullable(0.2)

julia> describe(@darr(a=[1,2,3,4,NA],b=[1,2,3,4,5]))
2 LabeledArray

  |min q1   med q3   max mean std                count nacount naratio
--+--------------------------------------------------------------------
a |1   1.75 2.5 3.25 4   2.5  1.2909944487358056 5     1       0.2
b |1   2.0  3.0 4.0  5   3.0  1.5811388300841898 5     0       0.0

julia> describe(@larr(a=[1,2,3,4,NA],b=[1,2,3,4,5],axis1[:m,:n,:p,:q,:r]))
2 LabeledArray

  |min q1   med q3   max mean std                count nacount naratio
--+--------------------------------------------------------------------
a |1   1.75 2.5 3.25 4   2.5  1.2909944487358056 5     1       0.2
b |1   2.0  3.0 4.0  5   3.0  1.5811388300841898 5     0       0.0
```

"""
function describe end

describe{T<:Nullable}(arr::AbstractArray{T}) = begin
  LDict(
    :min => try minimum(arr) catch Nullable{Float64}() end,
    :q1 => try quantile(arr, 0.25) catch Nullable{Float64}() end,
    :med => try quantile(arr, 0.5) catch Nullable{Float64}() end,
    :q3 => try quantile(arr, 0.75) catch Nullable{Float64}() end,
    :max => try maximum(arr) catch Nullable{Float64}() end,
    :mean => try mean(arr) catch Nullable{Float64}() end,
    :std => try std(arr) catch Nullable{Float64}() end,
    :count => Nullable(length(arr)),
    :nacount => Nullable(length(arr)-length(dropna(arr))),
    :naratio => Nullable((length(arr)-length(dropna(arr)))/length(arr)))
end

describe(arr::DictArray) = begin
  stat_func(v) = darr(
    :min => [try minimum(v) catch Nullable{Float64}() end],
    :q1 => try quantile(v, 0.25) catch Nullable{Float64}() end,
    :med => try quantile(v, 0.5) catch Nullable{Float64}() end,
    :q3 => try quantile(v, 0.75) catch Nullable{Float64}() end,
    :max => try maximum(v) catch Nullable{Float64}() end,
    :mean => try mean(v) catch Nullable{Float64}() end,
    :std => try std(v) catch Nullable{Float64}() end,
    :count => Nullable(length(v)),
    :nacount => Nullable(length(v)-length(dropna(v))),
    :naratio => Nullable((length(v)-length(dropna(v)))/length(v)))
  larr(cat(1, [stat_func(v) for (_,v) in peel(arr)]...), axis1=keys(arr))
end

describe(arr::LabeledArray) = describe(peel(arr))

"""

`shift(arr, offsets... [; isbound=false])`

Parallel shift the input array `arr` so that the element at `[1,...,1]` in `arr` shows up at `[1,...,1]+offsets` in the return array.

##### Arguments

* `arr`: `AbstractArrayWrapper`/`LabeledArray`/`DictArray`, the input array. When applied to `DictArray`, `shift` is applied to each field. When applied to `LabeledArray`, `shift` is applied to the base.
* `offsets`: integers to denote the amount of offset for each direction. It is assumed that there is no shift in the missing direcitons.
* `isbound`: default `false`. If `true`, the index is floored and capped between 1 and the maximum possible index along that direction. If `false`, any out of bound index due to shifting results in a nullable element.

##### Examples

```julia
julia> shift(nalift([1 2 3;4 5 6;7 8 9]), 1, 1)
3x3 DataCubes.AbstractArrayWrapper{Nullable{Int64},2,Array{Nullable{Int64},2}}:
 Nullable(5)        Nullable(6)        Nullable{Int64}()
 Nullable(8)        Nullable(9)        Nullable{Int64}()
 Nullable{Int64}()  Nullable{Int64}()  Nullable{Int64}()

julia> shift(nalift([1 2 3;4 5 6;7 8 9]), 1, 01)
3x3 DataCubes.AbstractArrayWrapper{Nullable{Int64},2,Array{Nullable{Int64},2}}:
 Nullable(5)        Nullable(6)        Nullable{Int64}()
 Nullable(8)        Nullable(9)        Nullable{Int64}()
 Nullable{Int64}()  Nullable{Int64}()  Nullable{Int64}()

julia> shift(nalift([1 2 3;4 5 6;7 8 9]), 1, 1)
3x3 DataCubes.AbstractArrayWrapper{Nullable{Int64},2,Array{Nullable{Int64},2}}:
 Nullable(5)        Nullable(6)        Nullable{Int64}()
 Nullable(8)        Nullable(9)        Nullable{Int64}()
 Nullable{Int64}()  Nullable{Int64}()  Nullable{Int64}()

julia> shift(nalift([1 2 3;4 5 6;7 8 9]), 1, -1)
3x3 DataCubes.AbstractArrayWrapper{Nullable{Int64},2,Array{Nullable{Int64},2}}:
 Nullable{Int64}()  Nullable(4)        Nullable(5)
 Nullable{Int64}()  Nullable(7)        Nullable(8)
 Nullable{Int64}()  Nullable{Int64}()  Nullable{Int64}()

julia> shift(darr(a=[1 2 3;4 5 6;7 8 9]), 1, -1)
3 x 3 DictArray

a |a |a
--+--+--
  |4 |5
  |7 |8
  |  |

julia> shift(larr(a=[1 2 3;4 5 6;7 8 9], axis2=[:X,:Y,:Z]), 1, -1)
3 x 3 LabeledArray

  |X |Y |Z
--+--+--+--
  |a |a |a
--+--+--+--
1 |  |4 |5
2 |  |7 |8
3 |  |  |
```

"""
function shift end

@generated shift{T,N}(arr::AbstractArrayWrapper{Nullable{T},N}, offsets::NTuple{N,Int};isbound=false) = quote
  result = similar(arr)
  sizearr = size(arr)
  if isbound
    @nloops $N i arr d->j_d=max(1,min(sizearr[d],i_d+offsets[d])) begin
      @nref($N,result,i) = @nref($N,arr,j)
    end
  else
    @nloops $N i arr d->j_d=i_d+offsets[d] begin
      @nref($N,result,i) = if @nall($N, d->checkbounds(Bool, sizearr[d], j_d))
        @nref($N,arr,j)
      else
        Nullable{T}()
      end
    end
  end
  result
end

@generated shift{T,N,A}(arr::AbstractArrayWrapper{Nullable{T},N,FloatNAArray{T,N,A}}, offsets::NTuple{N,Int};isbound=false) = quote
  arradata = arr.a.data
  result = similar(arradata)
  sizearr = size(arradata)
  na = convert(T, NaN)
  if isbound
    @nloops $N i arradata d->j_d=max(1,min(sizearr[d],i_d+offsets[d])) begin
      @nref($N,result,i) = @nref($N,arradata,j)
    end
  else
    @nloops $N i arr d->j_d=i_d+offsets[d] begin
      @nref($N,result,i) = if @nall($N, d->checkbounds(Bool, sizearr[d], j_d))
        @nref($N,arradata,j)
      else
        na
      end
    end
  end
  AbstractArrayWrapper(FloatNAArray(result))
end
shift{T}(arr::AbstractArrayWrapper{Nullable{T}}, offsets::Integer...;isbound=false) = begin
  lenoffsets = length(offsets)
  shift(arr, ntuple(d->lenoffsets<d ? 0:offsets[d], ndims(arr));isbound=isbound)
end
shift(arr::DictArray, offsets::Integer...;isbound=false) = DictArray(mapvalues(v->shift(v, offsets...;isbound=isbound), arr.data))
shift(arr::LabeledArray, offsets::Integer...;isbound=false) = LabeledArray(shift(arr.data, offsets...;isbound=isbound), arr.axes)


text:
using OrdinaryDiffEq, DiffEqSensitivity, Zygote

tspan = (0., 1.)
X = randn(3, 4)
p = randn(3, 4)
f(u,p,t) = u .* p
f(du,u,p,t) = (du .= u .* p)
prob_ube = ODEProblem{false}(f, X, tspan, p)
Zygote.gradient(p->sum(solve(prob_ube, Midpoint(), u0 = X, p = p)),p)

prob_ube = ODEProblem{true}(f, X, tspan, p)
Zygote.gradient(p->sum(solve(prob_ube, Midpoint(), u0 = X, p = p)),p)

function aug_dynamics!(dz, z, K, t)
    x = @view z[2:end]
    u = -K * x
    dz[1] = x' * x + u' * u
    dz[2:end] = x + u
end

policy_params = ones(2, 2)
z0 = zeros(3)
fwd_sol = solve(
    ODEProblem(aug_dynamics!, z0, (0.0, 1.0), policy_params),
    Tsit5(),
    u0 = z0,
    p = policy_params)

solve(
    ODEAdjointProblem(
        fwd_sol,
        InterpolatingAdjoint(),
        (out, x, p, t, i) -> (out .= 1),
        [1.0],
    ),Tsit5()
)

A = ones(2, 2)
B = ones(2, 2)
Q = ones(2, 2)
R = ones(2, 2)

function aug_dynamics!(dz, z, K, t)
    x = @view z[2:end]
    u = -K * x
    dz[1] = x' * Q * x + u' * R * u
    dz[2:end] = A * x + B * u # or just `x + u`
end

policy_params = ones(2, 2)
z0 = zeros(3)
fwd_sol = solve(
    ODEProblem(aug_dynamics!, z0, (0.0, 1.0), policy_params),
    u0 = z0,
    p = policy_params,
)

solve(
    ODEAdjointProblem(
        fwd_sol,
        InterpolatingAdjoint(),
        (out, x, p, t, i) -> (out .= 1),
        [1.0],
    ),
)

ext: jl


text:
## requirements in this file

using LinearAlgebra

## monocular camera on car tools

include(joinpath(@__DIR__, "CarSlamUtilsCommon.jl"))

##

# find joystick command segment
# get all inputs up to this next pose
function whichJoystickValues(fec, fromT, upToTime::ZonedDateTime)
  # upToTime = msgTime # getTimestamp(getVariable(fec.slam.dfg,addHist[end]))
  cmdTimes = (x->fec.synchronizer.cmdVal[x][2]).(1:length(fec.synchronizer.cmdVal)) .|> nanosecond2datetime
  addHist = getAddHistory(fec.slam.dfg)
  filter!(x->occursin(r"x\d",string(x)), addHist)
  mask = fromT .<= cmdTimes .< upToTime

  idxes = (1:length(fec.synchronizer.cmdVal))[mask]
  # get the actual times from the buffer
  joyValsTuple = fec.synchronizer.cmdVal[idxes]
  # @show typeof(joyValsTuple)
  # @show keys(joyValsTuple[1][3])
  # @show typeof(joyValsTuple[1][3])
  # @show joyValsTuple[1][3].axis
  # (x->x[3][:axis][[2;4]]).(joyValsTuple)
  # (x->x[3]).(joyValsTuple)

  # Tuple{Int,Int,Dict{Symbol, Any}}
  joyValsTuple
end

# for FluxModelsPose2Pose2
# ::AbstractArray{<:Real}
@everywhere function interpToOutx4(indata::AbstractMatrix{<:Real}; outlen=25)
  #
  len = size(indata,1)
  if 0 == len
    return zeros(outlen,4)
  end

  # assume evenly spaced
  idx = 1:len
  # tsLcl = range(indata[1,1],indata[end,1],length=outlen)
  iThrtTemp = DataInterpolations.LinearInterpolation(indata[:,1], idx)
  iSterTemp = DataInterpolations.LinearInterpolation(indata[:,2], idx)
  iVelxTemp = DataInterpolations.LinearInterpolation(indata[:,3], idx)
  iVelyTemp = DataInterpolations.LinearInterpolation(indata[:,4], idx)
  newArr = Matrix{Float64}(undef, outlen, 4)
  lookup = range(1,len, length=outlen)
  for i in 1:outlen
    lid = lookup[i]
    newArr[i, 1] = iThrtTemp(lid)
    newArr[i, 2] = iSterTemp(lid)
    newArr[i, 3] = iVelxTemp(lid)
    newArr[i, 4] = iVelyTemp(lid)
  end
  # currently have no velocity values
  return newArr
end

@everywhere function interpToOutx2(indata::AbstractMatrix{<:Real}; outlen=25)
  #
  @show typeof(indata)
  @show size(indata)
  len = size(indata,1)
  if 0 == len
    return zeros(outlen,2)
  end

  # assume evenly spaced
  idx = 1:len
  # tsLcl = range(indata[1,1],indata[end,1],length=outlen)
  iThrtTemp = DataInterpolations.LinearInterpolation(indata[:,1], idx)
  iSterTemp = DataInterpolations.LinearInterpolation(indata[:,2], idx)
  newArr = Matrix{Float64}(undef, outlen, 2)
  lookup = range(1,len, length=outlen)
  for i in 1:outlen
    lid = lookup[i]
    newArr[i, 1] = iThrtTemp(lid)
    newArr[i, 2] = iSterTemp(lid)
  end
  # currently have no velocity values
  return newArr
end

@everywhere function interpToOutx2(indata::Vector{<:AbstractVector{<:Real}}; outlen=25)
  if length(indata) == 1
    @warn "interpToOutx2 indata is too short, returning zeros"
    return zeros(25,2)
  end
  arr = zeros(length(indata),2)
  for i in 1:length(indata), j in 1:2
    arr[i, j] = indata[i][j]
  end
  interpToOutx2(arr, outlen=outlen)
end

@everywhere function interpToOutx2(indata::Vector{Any}; outlen=25)
  if length(indata) == 0
    return zeros(25,2)
  else
    @show indata
    error("interpToOutx2(indata::Vector{Any} does not know how to work with (see above)")
  end
end

# @everywhere function interpTo25x4Old(lclJD)
#   @warn "Unsure of interpTo25x4Old is doing the right thing (this is the old version)"
#   if 0 == size(lclJD,1)
#     return [zeros(4) for i in 1:25]
#   end
#
#   tsLcl = range(lclJD[1,1],lclJD[end,1],length=25)
#   intrTrTemp = DataInterpolations.LinearInterpolation(lclJD[:,2],lclJD[:,1])
#   intrStTemp = DataInterpolations.LinearInterpolation(lclJD[:,3],lclJD[:,1])
#   newVec = Vector{Vector{Float64}}()
#   for tsL in tsLcl
#     newVal = zeros(4)
#     newVal[1] = intrTrTemp(tsL)
#     newVal[2] = intrStTemp(tsL)
#     push!(newVec, newVal)
#   end
#   # currently have no velocity values
#   return newVec
# end

# Is this a general requirement?
function updateSLAMMono!(fec::FrontEndContainer,
                         WEIRDOFFSET::Dict,
                         syncList::Vector{Symbol}=fec.synchronizer.syncList;
                         useFluxModels::Bool=false,
                          allModels=nothing,
                          naiveFrac::Float64=0.6)
  #
  syncz = fec.synchronizer
  idxL,idxO = findSyncLatestIdx(syncz, syncList=[:leftFwdCam;:camOdo], weirdOffset=WEIRDOFFSET)
  if idxL == 0 || idxO == 0
    return nothing
  end
  # already have a counter
  fec.slam.frameCounter = syncz.camOdo[idxO][1]
  img = syncz.leftFwdCam[idxO][3]
  # @assert img isa Array "img is not an array, $(typeof(syncz.leftFwdCam[idxO]))"

  # calculate delta odo since last pose
  zTi = fec.slam.helpers.lastPoseOdomBuffer
  zT = syncz.camOdo[idxO][3]
  iT = zTi\zT
  # x-fwd, y-right, z-up
  psi = convert(Euler, iT.R).Y  # something not right here?
  DX = [iT.t[1];iT.t[2]; psi]

  # Tigger a new pose or skip out
  if !(0 < abs.(DX[3]) - parsed_args["pose_trigger_rotate"] ||
       0 < norm(DX[1:2]) - parsed_args["pose_trigger_distance"])
    # havent travelled far enough yet
    return nothing
  end
  @show round.(DX, digits=3)
  # psi = convert(Euler, zT.R).Y - convert(Euler, zTi.R).Y |> TU.wrapRad
  # @show convert(Euler, iT.R).Y  # something not right here?

  # reset interpose odometry counter for outside global odo integration process
  fec.slam.helpers.lastPoseOdomBuffer = syncz.camOdo[idxO][3] |> deepcopy

  # add new pose
  prevpose = Symbol("x$(fec.slam.poseCount)")
  fec.slam.poseCount += 1
  newpose = Symbol("x$(fec.slam.poseCount)")
  # take timestamp from image msg
  imgTime = nanosecond2datetime(syncz.leftFwdCam[idxL][2])
  addVariable!(fec.slam.dfg, newpose, Pose2, timestamp=imgTime)

  # previous pose time
  prevPoseT = getTimestamp(getVariable(fec.slam.dfg,prevpose))

  # find joystick values required
  joyVals = whichJoystickValues(fec, prevPoseT, imgTime)
  ## lets store the values while we're here
  addData!( fec.slam.dfg, :default_folder_store, prevpose, :JOYSTICK_CMD_VALS, Vector{UInt8}(JSON2.write( joyVals )), mimeType="application/json/octet-stream" )
  # addDataEntry!( fec.slam.dfg, prevpose, fec.datastore, :JOYSTICK_CMD_VALS, "application/json", Vector{UInt8}(JSON2.write( joyVals ))  )
  io = IOBuffer()
  pngSm = Stream(format"PNG", io)
  FileIO.save(pngSm, img)  # think FileIO is required for this
  addData!( fec.slam.dfg, :default_folder_store, prevpose, :KEYFRAME_IMG, take!(io), mimeType="image/png" )

  baselineOdo = MvNormal(DX, diagm([0.05; 0.05; 0.3].^2))
  pp = if useFluxModels
    # convert to NN required format
    # joyVelVal = catJoyVelData(joyVals, velVal)
    # interpolate joystick values to correct length
      # 2:throttle, 4:steering
    joyAs25x2 = ( x->[ x[3][:axis][2]; x[3][:axis][4] ] ).(joyVals)
      @show typeof(joyAs25x2), size(joyAs25x2)
    joyVals25x2 = interpToOutx2( joyAs25x2 )
    # @show size(joyVals25x2)
    joyVals25 = hcat(joyVals25x2, zeros(25,2))
    # build the factor
    FluxModelsPose2Pose2(allModels, joyVals25, baselineOdo, naiveFrac)
  else
    # TODO add covariance later
    Pose2Pose2(baselineOdo)
  end
  # take timestamp from image msg
  addFactor!(fec.slam.dfg, [prevpose;newpose], pp,
             timestamp=nanosecond2datetime(syncz.camOdo[idxO][2]) )

  # variables that can be initialized / solved
  put!(fec.slam.solveSettings.solvables, [newpose; ls(fec.slam.dfg,newpose)[1]])

  # add any potential landmarks and factors
  tagsL = detector(syncz.leftFwdCam[idxL][3] |> collect)
  # tagsR = detector(syncHdlr[:right][idxR][2])

  @show tagsL |> length

  for tag in tagsL
    @show lmid = Symbol("l$(tag.id)")
    !exists(slam.dfg, lmid) && addVariable!(slam.dfg, lmid, Pose2)
    # pose = homographytopose(tag.H, fx, fy, cx, cy, taglength = 160.)
    pose, err1 = tagOrthogonalIteration(tag, fx, fy, cx, cy, taglength = 0.25)
    cTt = LinearMap(pose[1:3, 1:3])∘Translation((pose[1:3,4])...)
    bRc = Rotations.Quat(1/sqrt(2),0,0,-1/sqrt(2))*Rotations.Quat(1/sqrt(2),-1/sqrt(2),0,0)
    bTt = LinearMap(bRc) ∘ cTt
    # wTb = LinearMap(zT.R.R) ∘ Translation(zT.t...)
    # wTt = wTb ∘ bTt

    ld = LinearAlgebra.cross(bTt.linear*[0;0;1], [0;0;1])
    theta = atan(ld[2],ld[1])
    Dtag = [bTt.translation[1:2,];theta]
    # println("$(tag.id), $(ld[3]), $(round.(Dtag, digits=3))")
    p2l2 = Pose2Pose2(MvNormal(Dtag,diagm([0.1;0.1;0.1].^2)))
    addFactor!(slam.dfg, [newpose;lmid], p2l2, timestamp=imgTime)

    # FIXME should be in FEC
    if parsed_args["vis3d"]
      setobject!(vis[:tags][Symbol("tag_$(tag.id)")][:triad], Triad(0.2))
      tagobj = HyperRectangle(Vec(0.,0,0), Vec(0.3,0.3,0.03))
      setobject!(vis[:tags][Symbol("tag_$(tag.id)")][:plate], tagobj)
      settransform!(vis[:tags][Symbol("tag_$(tag.id)")], bTt) # bTt
    end

    # setobject!(vis[:poses][Symbol(newpose)], Triad(0.5))
    # settransform!(vis[:poses][Symbol(newpose)], wTb)

  end

  # check whether a solve should be trigger
  checkSolveStrideTrigger!(fec.slam)

  nothing
end

function drawLatestImage(fec::FrontEndContainer; syncList=[:leftFwdCam;])
  # get synced images
  idxL = findSyncLatestIdx(fec.synchronizer, weirdOffset=WEIRDOFFSET, syncList=syncList )

  # @info "drawLatestImagePair, $idxL, $idxR, $idxO"
  if idxL == 0
    return nothing
  end

  # Do apriltag detection
  tagsL = detector(fec.synchronizer.leftFwdCam[idxL][3] |> collect)

  # @show poses = (T->homographytopose(T.H, fx, fy, cx, cy, taglength = 160.)).(tagsL)
  imgLt = AprilTags.drawTags(fec.synchronizer.leftFwdCam[idxL][3], K, tagsL) # showImage

  # imgL = syncImgs[:left][idxL][2] .|> Gray
  # imgR = syncImgs[:right][idxR][2] .|> Gray
  # featuresL = Keypoints(fastcorners(imgL, 12, 0.25))
  # featuresR = Features(fastcorners(imgR, 12, 0.25))
  # offsetx = CartesianIndex(0, 5)
  # offsety = CartesianIndex(5, 0)
  # map(m -> draw!(imgLt, LineSegment(m - offsetx, m + offsetx)), featuresL)
  # map(m -> draw!(imgLt, LineSegment(m - offsety, m + offsety)), featuresL)

  # draw both
  if parsed_args["imshow"]
    imshow(canvases[1,1], imgLt)
    Gtk.showall(gui["window"])
  end
  nothing
end

##  MessageHandler

function leftImgHdlr(msgdata, fec::FrontEndContainer)
  # @show "leftImgHdlr", msgdata[2].header.seq
  # @show "leftImgHdlr"
  mtms = getROSPyMsgTimestamp(msgdata[end])
  nsT = Int64(datetime2unix(mtms[1]))*1000_000_000 + mtms[2]

  leftdata = take!(IOBuffer(msgdata[2].data))

  # more direct image from ROS and zed
  nvi = normedview(leftdata)
  # im3 = last(fec.synchronizer.leftFwdCam)[2]
  arr = reshape(nvi, 3, 672,376) # im3
  pv = PermutedDimsArray(arr, (1,3,2))
  pvn = similar(pv)
  pvn[3,:,:] .= pv[1,:,:]
  pvn[2,:,:] .= pv[2,:,:]
  pvn[1,:,:] .= pv[3,:,:]
  ldImg = colorview(RGB, pvn)

  # for compressed images
  # ldImg = ImageMagick.load_(leftdata)
  push!(fec.synchronizer.leftFwdCam, (msgdata[2].header.seq, nsT, ldImg) )
  # img = last(SyncImages[:left])
  nothing
end

function odomHdlr(msgdata, fec::FrontEndContainer, WO, allModels=nothing)
  # @show "odomHdlr", msgdata[2].header.seq
  # @show "odomHdlr"
  mtms = getROSPyMsgTimestamp(msgdata[end])
  nsT = Int64(datetime2unix(mtms[1]))*1000_000_000 + mtms[2]

  pos = msgdata[2].pose.pose.position
  quat = msgdata[2].pose.pose.orientation

  zPi = [pos.x;pos.y;pos.z]
  zQi = TU.Quaternion(quat.w,[quat.x;quat.y;quat.z])

  zTi = SE3(zPi, zQi)
  push!(fec.synchronizer.camOdo, (msgdata[2].header.seq, nsT, zTi) )

  # show leftFwdCam image
  drawLatestImage(fec, syncList=[:leftFwdCam;])

  # update the factor graph object as required
  if allModels == nothing
    updateSLAMMono!(fec, WO, useFluxModels=false )
  else
    updateSLAMMono!(  fec, WO, useFluxModels=true, allModels=allModels, naiveFrac=parsed_args["naive_frac"]  )
  end

  nothing
end

function joystickHdlr(msgdata, fec::FrontEndContainer)
  # @show "joystickHdlr", msgdata[2].header.seq
  mtms = getROSPyMsgTimestamp(msgdata[end])
  nsT = Int64(datetime2unix(mtms[1]))*1000_000_000 + mtms[2]

  data = Dict(:axis => msgdata[2].axes,
              :buttons => msgdata[2].buttons)
  # @show data
  push!(fec.synchronizer.cmdVal, (msgdata[2].header.seq, nsT, data) )

  #
end

#


text:
"Meta data to attach to ω::Ω"
const Tags = NamedTuple

function mergetag end

# "tag value `val` with tag `tag` and merge tags"
hastag(::Type{Tags{K, V}}, tag::Symbol) where {K, V} = tag in K
hastag(::Type{Tags{K, V}}, tag::Type{Val{S}}) where {K, V, S} = S in K

"Trait type to denote Tagged type with tag"
struct HasTag{T} end

"Trait type to denote Tagged type with absence of tag"
struct NotHasTag{T} end

"`tag(x, tags)` tag value `x` with tags `tags`"
function tag end

"""
Trait function -- `traithastag(t, Val{:sometag})` returns `HasTag{:sometag}`
if `t` has that tag or `NotHasTag{:sometag}` otherwise
"""
traithastag(t::Type{T}, ::Type{Val{S}}) where {T <: Tags, S} = hastag(T, S) ? HasTag{S}() : NotHasTag{S}()
traithastag(t::T, s::Type{Val{S}}) where {T <: Tags, S} = traithastag(T, s)


text:
const TOKEN_PREFIX = '/'
"""
    JSONPointer(token)

Follows https://tools.ietf.org/html/rfc6901 standard

### Differences are
- Index numbers starts from '1' instead of '0'
- User can declare type with '::T' notation at the end

### Example
"""
struct JSONPointer{T}
    token::Tuple
end
function JSONPointer(token::AbstractString)
    if !startswith(token, TOKEN_PREFIX)
        throw(ArgumentError("JSONPointer must starts with '/' prefix"))
    end

    T = Any
    jk = convert(Array{Any, 1}, split(chop(token; head=1, tail=0), TOKEN_PREFIX))
    if occursin("::", jk[end])
        x = split(jk[end], "::")
        jk[end] = x[1]
        T = (x[2] == "Vector" ? "Vector{Any}" : x[2]) |> Meta.parse |> eval
    end
    @inbounds for i in 1:length(jk)
        if occursin(r"^\d+$", jk[i]) # index of a array
            jk[i] = parse(Int, string(jk[i]))
            if iszero(jk[i])
                throw(AssertionError("Julia uses 1-based indexing"))
            end
        end
    end

    JSONPointer{T}(tuple(jk...))
end

"""
    null_value(p::JSONPointer{T}) where T

provide appropriate value for 'T'
'Real' return 'zero(T)' and 'AbstractString' returns '""'

If user wants different null value for 'T' override 'null_value(p::JSONPointer{T})' method

"""
null_value(p::JSONPointer{T}) where T = missing
function null_value(p::JSONPointer{T}) where T <: Array
    eltype(T) <: Real ? eltype(T)[] :
    eltype(T) <: AbstractString ? eltype(T)[] :
    Any[]
end

Base.Dict(p::JSONPointer) = create_by_pointer(Dict, p)
DataStructures.OrderedDict(p::JSONPointer) = create_by_pointer(OrderedDict, p)

function create_by_pointer(::Type{T}, p::JSONPointer) where T <: AbstractDict
    val = missing

    @inbounds for i in length(p):-1:1
        k = p.token[i]
        if isa(k, Integer)
            tmp = Array{Any, 1}(missing, k)
            tmp[k] = if i == length(p)
                        null_value(p)
                    else
                        val
                    end
            val = tmp
        elseif isa(k, AbstractString)
            val = if i == length(p)
                    T{String, Any}(k => null_value(p))
                else
                    T{String, Any}(k => val)
                end
        end
    end
    return val
end

function create_by_pointer(::Type{T}, arr::Array) where T <: AbstractDict
    template = T(arr[1])
    if length(arr) > 1
        @inbounds for p in arr
            p2 = JSONPointer{Any}(p.token)
            template[p2] = null_value(p)
        end
    end
    return template
end

for T in (Dict, OrderedDict)
    @eval Base.haskey(dict::$T{K,V}, p::JSONPointer) where {K, V} = haskey_by_pointer(dict, p)
    @eval Base.getindex(dict::$T{K,V}, p::JSONPointer) where {K, V} = getindex_by_pointer(dict, p)
    @eval Base.setindex!(dict::$T{K,V}, v, p::JSONPointer) where {K <: AbstractString, V} = setindex_by_pointer!(dict, v, p)
    @eval Base.setindex!(dict::$T{K,V}, v, p::JSONPointer) where {K <: Integer, V} = setindex_by_pointer!(dict, v, p)
end
function haskey_by_pointer(collection, p::JSONPointer)::Bool
    b = true
    val = collection
    @inbounds for (i, k) in enumerate(p.token)
        val = begin
            if isa(val, Array)
                if !isa(k, Integer)
                    missing
                else
                    length(val) >= k ? val[k] : missing
                end
            else
                if isa(k, Integer)
                    missing
                else
                    haskey(val, k) ? val[k] : missing
                end
            end
        end

        if ismissing(val)
            b = false
            break
        end
    end
    return b
end
function getindex_by_pointer(collection, p::JSONPointer, i = 1)
    val = getindex(collection, p.token[i])
    if i < length(p)
        val = getindex_by_pointer(val, p, i+1)
    end
    return val
end

function setindex_by_pointer!(collection::T, v, p::JSONPointer{U}) where {T <: AbstractDict, U}
    v = ismissing(v) ? null_value(p) : v
    if !isa(v, U) &&
        try
            v = convert(eltype(p), v)
        catch e
            msg = isa(v, Array) ? "Vector" : "Any"
            error("$v is not valid value for $p use '::$msg' if you don't need static type")
            throw(e)
        end
    end
    prev = collection
    DT = OrderedDict{String, Any}

    @inbounds for (i, k) in enumerate(p.token)
        if isa(prev, AbstractDict)
            DT = typeof(prev)
        end
        if isa(prev, Array)
            if !isa(k, Integer)
                throw(MethodError(setindex!, k))
            end
            grow_array!(prev, k)
        else
            if isa(k, Integer)
                throw(MethodError(setindex!, k))
            end
            if !haskey(prev, k)
                setindex!(prev, missing, k)
            end
        end

        if i < length(p)
            tmp = getindex(prev, k)
            if ismissing(tmp)
                next_key = p.token[i+1]
                if isa(next_key, Integer)
                    new_data = Array{Any,1}(missing, next_key)
                else
                    new_data = DT(next_key => missing)
                end
                setindex!(prev, new_data, k)
            end
            prev = getindex(prev, k)
        end
    end
    setindex!(prev, v, p.token[end])
end

function grow_array!(arr, target_size)
    x = target_size - length(arr)
    if x > 0
        T = eltype(arr)
        new_arr = Array{T, 1}(undef, x)
        new_arr .= T <: Real ? zero(T) : missing
        append!(arr, new_arr)
    end
    return arr
end

Base.length(x::JSONPointer) = length(x.token)
Base.eltype(x::JSONPointer{T}) where T = T

macro j_str(token)
    JSONPointer(token)
end

function Base.show(io::IO, x::JSONPointer{T}) where T
    print(io,
    "JSONPointer{", T, "}(\"/", join(x.token, "/"), "\")")
end


text:
"""
    npdfnc(
        npd :: NASAPrecipitationDataset,
        geo :: GeoRegion,
        dt  :: TimeType
    ) -> String

Returns of the path of the file for the NASA Precipitation dataset specified by `npd` for a GeoRegion specified by `geo` at a date specified by `dt`.

Arguments
=========
- `npd` : a `NASAPrecipitationDataset` specifying the dataset details and date download range
- `geo` : a `GeoRegion` (see [GeoRegions.jl](https://github.com/JuliaClimate/GeoRegions.jl)) that sets the geographic bounds of the data array in lon-lat
- `dt`  : A specified date. The NCDataset retrieved may will contain data for the date, although it may also contain data for other dates depending on the `NASAPrecipitationDataset` specified by `npd`
"""
function npdfnc(
    npd :: Union{IMERGHalfHourly{ST,DT},TRMM3Hourly{ST,DT}},
	geo :: GeoRegion,
    dt  :: TimeType
) where {ST<:AbstractString, DT<:TimeType}

    fol = joinpath(npd.sroot,geo.regID,yrmo2dir(dt))
    fnc = npd.npdID * "-" * geo.regID * "-" * ymd2str(dt) * ".nc"
    return joinpath(fol,fnc)

end

function npdfnc(
    npd :: Union{IMERGDaily{ST,DT},TRMMDaily{ST,DT}},
	geo :: GeoRegion,
    dt  :: TimeType
) where {ST<:AbstractString, DT<:TimeType}

    fol = joinpath(npd.sroot,geo.regID,yr2str(dt))
    fnc = npd.npdID * "-" * geo.regID * "-" * yrmo2str(dt) * ".nc"
    return joinpath(fol,fnc)

end

function npdfnc(
    npd :: Union{IMERGMonthly{ST,DT},TRMMMonthly{ST,DT}},
	geo :: GeoRegion,
    dt  :: TimeType
) where {ST<:AbstractString, DT<:TimeType}

    fol = joinpath(npd.sroot,geo.regID)
    fnc = npd.npdID * "-" * geo.regID * "-" * yr2str(dt) * ".nc"
    return joinpath(fol,fnc)

end

max_stars_repo_name: natgeo-wong/NASAPrecipitation

max_issues_repo_name: natgeo-wong/NASAPrecipitation

max_forks_repo_name: natgeo-wong/NASAPrecipitation