Support overlapping tiles

.travis.yml

@@ -3,8 +3,8 @@ language: julia
 os:
   - linux
   - osx
+  - windows
 julia:
-  - 0.7
   - 1.0
   - 1
   - nightly

Project.toml

@@ -8,10 +8,11 @@ OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
 
 [compat]
 OffsetArrays = "0.8, 0.9, 0.10, 0.11, 1"
-julia = "0.7, 1"
+julia = "1"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 
 [targets]
-test = ["Test"]
+test = ["Test", "Documenter"]

docs/src/note.txt

@@ -0,0 +1 @@
+Documenter needs this directory to exist. Otherwise it refuses to run doctests.

src/TiledIteration.jl

@@ -11,55 +11,13 @@ else
     _inc(state, iter) = inc(state, iter.indices)
 end
 
-export TileIterator, EdgeIterator, padded_tilesize, TileBuffer
+export TileIterator, EdgeIterator, padded_tilesize, TileBuffer, TileIterator, RelaxStride, RelaxLastTileSize
 
-### TileIterator ###
+include("tileiterator.jl")
 
 const L1cachesize = 2^15
 const cachelinesize = 64
 
-struct TileIterator{N,I,UR}
-    inds::I
-    sz::Dims{N}
-    R::CartesianIndices{N,UR}
-end
-rangetype(::CartesianIndices{N,T}) where {N,T} = T
-function TileIterator(inds::Indices{N}, sz::Dims{N}) where N
-    ls = map(length, inds)
-    ns = map(ceildiv, ls, sz)
-    R = CartesianIndices(ns)
-    TileIterator{N,typeof(inds),rangetype(R)}(inds, sz, R)
-end
-
-Iterators.IteratorEltype(::Type{<:TileIterator}) = Iterators.HasEltype()
-
-ceildiv(l, s) = ceil(Int, l/s)
-
-Base.length(iter::TileIterator) = length(iter.R)
-Base.eltype(iter::TileIterator{N}) where {N} = NTuple{N,UnitRange{Int}}
-
-function Base.iterate(iter::TileIterator)
-    iterR = iterate(iter.R)
-    iterR === nothing && return nothing
-    I, state = iterR
-    return getindices(iter, I), state
-end
-function Base.iterate(iter::TileIterator, state)
-    iterR = iterate(iter.R, state)
-    iterR === nothing && return nothing
-    I, newstate = iterR
-    return getindices(iter, I), newstate
-end
-
-Base.show(io::IO, iter::TileIterator) = print(io, "TileIterator(", iter.inds, ", ", iter.sz, ')')
-
-@inline function getindices(iter::TileIterator, I::CartesianIndex)
-    map3(_getindices, iter.inds, iter.sz, I.I)
-end
-_getindices(ind, s, i) = first(ind)+(i-1)*s : min(last(ind),first(ind)+i*s-1)
-map3(f, ::Tuple{}, ::Tuple{}, ::Tuple{}) = ()
-@inline map3(f, a::Tuple, b::Tuple, c::Tuple) = (f(a[1], b[1], c[1]), map3(f, tail(a), tail(b), tail(c))...)
-
 ### EdgeIterator ###
 
 struct EdgeIterator{N,UR1,UR2}

src/tileiterator.jl

@@ -0,0 +1,224 @@
+################################################################################
+##### CoveredRange
+################################################################################
+
+# a range covered by subranges
+struct CoveredRange{R,S} <: AbstractVector{UnitRange{Int}}
+    offsets::R
+    stopping::S
+end
+
+struct FixedLength
+    length::Int
+end
+
+struct LengthAtMost{N <: Union{Int, Nothing}}
+    maxlength::Int
+    maxstop::N
+end
+
+function compute_stop(offset, stopping::FixedLength)
+    return offset+stopping.length
+end
+function compute_stop(offset, stopping::LengthAtMost)
+    return min(offset+stopping.maxlength, stopping.maxstop)
+end
+
+compute_range(offset, stopping)::UnitRange{Int} = (offset+1):compute_stop(offset, stopping)
+
+Base.size(o::CoveredRange,args...) = size(o.offsets, args...)
+Base.@propagate_inbounds function Base.getindex(o::CoveredRange, inds...)
+    offset = o.offsets[inds...]
+    return compute_range(offset, o.stopping)
+end
+
+################################################################################
+##### RoundedRange
+################################################################################
+struct RoundedRange{R} <: AbstractVector{Int}
+    range::R
+end
+
+Base.@propagate_inbounds function Base.getindex(r::RoundedRange, i)
+    rough = r.range[i]
+    return round(Int, rough)
+end
+
+function roundedrange(start; stop, length)
+    inner = LinRange(start, stop, length)
+    return RoundedRange(inner)
+end
+
+Base.size(o::RoundedRange, args...) = size(o.range, args...)
+
+################################################################################
+##### TileIterator
+################################################################################
+struct TileIterator{N,C} <: AbstractArray{NTuple{N, UnitRange{Int}}, N}
+    covers1d::C
+end
+
+function TileIterator(covers1d::NTuple{N, AbstractVector{UnitRange{Int}}}) where {N}
+    C = typeof(covers1d)
+    return TileIterator{N, C}(covers1d)
+end
+
+function TileIterator(axes::Indices{N}, tilesize::Dims{N}) where {N}
+    TileIterator(axes, RelaxLastTile(tilesize))::TileIterator{N}
+end
+
+"""
+    titr = TileIterator(axes::NTuple{N, AbstractUnitRange}, strategy)
+
+Decompose `axes` into an iterator `titr` of smaller axes according to `strategy`.
+
+The `strategy` argument controls the details of the tiling. For instance
+if the length of an axis is not divisible by the tile size, what should happen?
+One approach would be to relax the size requirement for the last tile.
+Another possibility to relax the `stride` so that all tiles are of the requested size,
+but tiles may be slightly overlapping.
+These two possibilities are implemented by [`RelaxLastTile`](@ref) and [`RelaxStride`](@ref).
+
+# Examples
+```jldoctest
+julia> using TiledIteration
+
+julia> collect(TileIterator((1:3, 0:5), RelaxLastTile((2, 3))))
+2×2 Array{Tuple{UnitRange{Int64},UnitRange{Int64}},2}:
+ (1:2, 0:2)  (1:2, 3:5)
+ (3:3, 0:2)  (3:3, 3:5)
+
+julia> collect(TileIterator((1:3, 0:5), (2, 3))) # defaults to RelaxLastTile
+2×2 Array{Tuple{UnitRange{Int64},UnitRange{Int64}},2}:
+ (1:2, 0:2)  (1:2, 3:5)
+ (3:3, 0:2)  (3:3, 3:5)
+
+julia> collect(TileIterator((1:3, 0:5), RelaxStride((2, 3))))
+2×2 Array{Tuple{UnitRange{Int64},UnitRange{Int64}},2}:
+ (1:2, 0:2)  (1:2, 3:5)
+ (2:3, 0:2)  (2:3, 3:5)
+```
+"""
+function TileIterator(axes, strategy)
+    covers1d = map(cover1d, axes, split(strategy))
+    return TileIterator(covers1d)
+end
+
+# strategies
+export RelaxStride
+"""
+    RelaxStride(tilesize)
+
+Tiling strategy, that guarantees each tile of size `tilesize`.
+If the needed tiles will slightly overlap, to cover everything.
+
+# Examples
+```jldoctest
+julia> using TiledIteration
+
+julia> collect(TileIterator((1:4,), RelaxStride((2,))))
+2-element Array{Tuple{UnitRange{Int64}},1}:
+ (1:2,)
+ (3:4,)
+
+julia> collect(TileIterator((1:4,), RelaxStride((3,))))
+2-element Array{Tuple{UnitRange{Int64}},1}:
+ (1:3,)
+ (2:4,)
+```
+
+See also [`TileIterator`](@ref).
+"""
+struct RelaxStride{N}
+    tilesize::Dims{N}
+end
+
+export RelaxLastTile
+
+"""
+    RelaxLastTile(tilesize)
+
+Tiling strategy, that premits the size of the last tiles along each dimension to be smaller
+then `tilesize` if needed. All other tiles are of size `tilesize`.
+
+# Examples
+```jldoctest
+julia> using TiledIteration
+
+julia> collect(TileIterator((1:4,), RelaxLastTile((2,))))
+2-element Array{Tuple{UnitRange{Int64}},1}:
+ (1:2,)
+ (3:4,)
+
+julia> collect(TileIterator((1:7,), RelaxLastTile((2,))))
+4-element Array{Tuple{UnitRange{Int64}},1}:
+ (1:2,)
+ (3:4,)
+ (5:6,)
+ (7:7,)
+```
+
+See also [`TileIterator`](@ref).
+"""
+struct RelaxLastTile{N}
+    tilesize::Dims{N}
+end
+
+"""
+    split(strategy)
+
+Split an N dimensional strategy into an NTuple of 1 dimensional strategies.
+"""
+function split end
+
+function split(strategy::RelaxStride)
+    map(strategy.tilesize) do s
+        RelaxStride((s,))
+    end
+end
+
+function split(strategy::RelaxLastTile)
+    map(strategy.tilesize) do s
+        RelaxLastTile((s,))
+    end
+end
+
+function cover1d(ax, strategy::RelaxStride{1})::CoveredRange
+    tilelen = stride = first(strategy.tilesize)
+    firstoffset = first(ax)-1
+    lastoffset = last(ax) - tilelen
+    stepcount = ceil(Int, (lastoffset - firstoffset) / stride) + 1
+    offsets = roundedrange(firstoffset, stop=lastoffset, length=stepcount)
+    @assert last(ax) - tilelen <= last(offsets) <= last(ax)
+
+    stopping = FixedLength(tilelen)
+    return CoveredRange(offsets, stopping)
+end
+
+function cover1d(ax, strategy::RelaxLastTile{1})::CoveredRange
+    tilelen = stride = first(strategy.tilesize)
+    maxstop = last(ax)
+    stopping = LengthAtMost(tilelen, maxstop)
+
+    lo = first(ax)
+    hi = last(ax)
+    stepcount = if tilelen <= stride
+        floor(Int, (hi - lo) / stride) + 1
+    else
+        ceil(Int, (hi + 1 - lo - tilelen) / stride) + 1
+    end
+    firstoffset = lo - 1
+    offsets = range(firstoffset, step=tilelen, length=stepcount)
+    return CoveredRange(offsets, stopping)
+end
+
+Base.@propagate_inbounds function Base.getindex(o::TileIterator, inds::Integer...)
+    cis = CartesianIndices(o)[inds...]
+    o[cis]
+end
+Base.@propagate_inbounds function Base.getindex(o::TileIterator, ci::CartesianIndex)
+    map(getindex, o.covers1d, Tuple(ci))
+end
+
+Base.size(o::TileIterator) = map(length, o.covers1d)
+Base.IndexStyle(o::TileIterator) = IndexCartesian()