Skip to content

rework OffsetArrays constructor #148

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 9 commits into from
Sep 23, 2020
Merged

rework OffsetArrays constructor #148

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
67977d8
[WIP] rework OffsetArrays constructor
johnnychen94 Sep 18, 2020
07f9b3c
apply suggestion from comments
johnnychen94 Sep 20, 2020
9166dfc
add more test cases for OffsetArray constructors
johnnychen94 Sep 20, 2020
622c487
more simplification
johnnychen94 Sep 20, 2020
5a804c9
remove unrelated changes
johnnychen94 Sep 20, 2020
88061ad
Revert "more simplification"
johnnychen94 Sep 21, 2020
ffebce9
apply suggestion from comments and fix a dispatching error
johnnychen94 Sep 23, 2020
b88fb74
more test coverage
johnnychen94 Sep 23, 2020
6469a8c
more test to improve coverage
johnnychen94 Sep 23, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
209 changes: 70 additions & 139 deletions src/OffsetArrays.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,62 +10,15 @@ end
export OffsetArray, OffsetMatrix, OffsetVector

include("axes.jl")
include("utils.jl")

## OffsetArray
struct OffsetArray{T,N,AA<:AbstractArray} <: AbstractArray{T,N}
parent::AA
offsets::NTuple{N,Int}
function OffsetArray{T, N, AA}(parent::AA, offsets::NTuple{N, Int}) where {T, N, AA<:AbstractArray}
overflow_check.(axes(parent), offsets)
new{T, N, AA}(parent, offsets)
end
end
OffsetVector{T,AA<:AbstractArray} = OffsetArray{T,1,AA}
OffsetMatrix{T,AA<:AbstractArray} = OffsetArray{T,2,AA}

function overflow_check(r, offset::T) where T
# This gives some performance boost https://github.com/JuliaLang/julia/issues/33273
throw_upper_overflow_error() = throw(ArgumentError("Boundary overflow detected: offset $offset should be equal or less than $(typemax(T) - last(r))"))
throw_lower_overflow_error() = throw(ArgumentError("Boundary overflow detected: offset $offset should be equal or greater than $(typemin(T) - first(r))"))

if offset > 0 && last(r) > typemax(T) - offset
throw_upper_overflow_error()
elseif offset < 0 && first(r) < typemin(T) - offset
throw_lower_overflow_error()
end
end

## OffsetArray constructors

offset(axparent::AbstractUnitRange, ax::AbstractUnitRange) = first(ax) - first(axparent)
offset(axparent::AbstractUnitRange, ax::Integer) = 1 - first(axparent)

function OffsetArray(A::AbstractArray{T,N}, offsets::NTuple{N,Int}) where {T,N}
OffsetArray{T,N,typeof(A)}(A, offsets)
end
OffsetArray(A::AbstractArray{T,0}, offsets::Tuple{}) where T =
OffsetArray{T,0,typeof(A)}(A, ())

OffsetArray(A::AbstractArray{T,N}, offsets::Vararg{Int,N}) where {T,N} =
OffsetArray(A, offsets)
OffsetArray(A::AbstractArray{T,0}) where {T} = OffsetArray(A, ())

# Technically we know the length of CartesianIndices but we need to convert it first, so here we
# don't put it in OffsetAxisKnownLength.
const OffsetAxisKnownLength = Union{Integer, AbstractUnitRange, IdOffsetRange}
const OffsetAxis = Union{OffsetAxisKnownLength, CartesianIndices, Colon}
const ArrayInitializer = Union{UndefInitializer, Missing, Nothing}
OffsetArray{T,N}(init::ArrayInitializer, inds::Indices{N}) where {T,N} =
OffsetArray(Array{T,N}(init, map(indexlength, inds)), map(indexoffset, inds))
OffsetArray{T}(init::ArrayInitializer, inds::Indices{N}) where {T,N} = OffsetArray{T,N}(init, inds)
OffsetArray{T,N}(init::ArrayInitializer, inds::Vararg{AbstractUnitRange,N}) where {T,N} = OffsetArray{T,N}(init, inds)
OffsetArray{T}(init::ArrayInitializer, inds::Vararg{AbstractUnitRange,N}) where {T,N} = OffsetArray{T,N}(init, inds)

# OffsetVector constructors
OffsetVector(A::AbstractVector, offset) = OffsetArray(A, offset)
OffsetVector{T}(init::ArrayInitializer, inds::AbstractUnitRange) where {T} = OffsetArray{T}(init, inds)

# OffsetMatrix constructors
OffsetMatrix(A::AbstractMatrix, offset1, offset2) = OffsetArray(A, offset1, offset2)
OffsetMatrix(A::AbstractMatrix, I::CartesianIndices{2}) = OffsetArray(A, I)
OffsetMatrix{T}(init::ArrayInitializer, inds1::AbstractUnitRange, inds2::AbstractUnitRange) where {T} = OffsetArray{T}(init, inds1, inds2)

## OffsetArray
"""
OffsetArray(A, indices...)

Expand Down Expand Up @@ -116,84 +69,74 @@ ERROR: [...]
```

"""
function OffsetArray(A::AbstractArray{T,N}, inds::NTuple{N,AbstractUnitRange}) where {T,N}
axparent = axes(A)
lA = map(length, axparent)
lI = map(length, inds)
lA == lI || throw(DimensionMismatch("supplied axes do not agree with the size of the array (got size $lA for the array and $lI for the indices"))
OffsetArray(A, map(offset, axparent, inds))
struct OffsetArray{T,N,AA<:AbstractArray} <: AbstractArray{T,N}
parent::AA
offsets::NTuple{N,Int}
function OffsetArray{T, N, AA}(parent::AA, offsets::NTuple{N, Int}) where {T, N, AA<:AbstractArray}
@boundscheck overflow_check.(axes(parent), offsets)
new{T, N, AA}(parent, offsets)
end
end
OffsetArray(A::AbstractArray{T,N}, inds::Vararg{AbstractUnitRange,N}) where {T,N} =
OffsetArray(A, inds)
const OffsetVector{T,AA<:AbstractArray} = OffsetArray{T,1,AA}
const OffsetMatrix{T,AA<:AbstractArray} = OffsetArray{T,2,AA}

uncolonindices(A::AbstractArray{<:Any,N}, inds::NTuple{N,Any}) where {N} = uncolonindices(axes(A), inds)
uncolonindices(ax::Tuple, inds::Tuple) = (first(inds), uncolonindices(tail(ax), tail(inds))...)
uncolonindices(ax::Tuple, inds::Tuple{Colon, Vararg{Any}}) = (first(ax), uncolonindices(tail(ax), tail(inds))...)
uncolonindices(::Tuple{}, ::Tuple{}) = ()
function overflow_check(r, offset::T) where T
# This gives some performance boost https://github.com/JuliaLang/julia/issues/33273
throw_upper_overflow_error() = throw(ArgumentError("Boundary overflow detected: offset $offset should be equal or less than $(typemax(T) - last(r))"))
throw_lower_overflow_error() = throw(ArgumentError("Boundary overflow detected: offset $offset should be equal or greater than $(typemin(T) - first(r))"))

function OffsetArray(A::AbstractArray{T,N}, inds::NTuple{N,Union{AbstractUnitRange, CartesianIndices{1}, Colon}}) where {T,N}
OffsetArray(A, uncolonindices(A, inds))
end
OffsetArray(A::AbstractArray{T,N}, inds::Vararg{Union{AbstractUnitRange, CartesianIndices{1}, Colon},N}) where {T,N} =
OffsetArray(A, uncolonindices(A, inds))

# Specify offsets using CartesianIndices (issue #71)
# Support a mix of AbstractUnitRanges and CartesianIndices{1}
# Extract the range r from CartesianIndices((r,))
function stripCartesianIndices(inds::Tuple{CartesianIndices{1},Vararg{Any}})
I = first(inds)
Ir = convert(Tuple{AbstractUnitRange{Int}}, I) |> first
(Ir, stripCartesianIndices(tail(inds))...)
end
stripCartesianIndices(inds::Tuple)= (first(inds), stripCartesianIndices(tail(inds))...)
stripCartesianIndices(::Tuple{}) = ()

OffsetArray(A::AbstractArray{<:Any,N}, inds::NTuple{N,Union{CartesianIndices{1}, AbstractUnitRange}}) where {N} =
OffsetArray(A, stripCartesianIndices(inds))
OffsetArray(A::AbstractArray{<:Any,N}, inds::Vararg{Union{CartesianIndices{1}, AbstractUnitRange},N}) where {N} =
OffsetArray(A, inds)

# Support an arbitrary CartesianIndices alongside colons and ranges
# The total number of indices should equal ndims(arr)
# We split the CartesianIndices{N} into N CartesianIndices{1} indices to facilitate dispatch
splitCartesianIndices(c::CartesianIndices{0}) = ()
function splitCartesianIndices(c::CartesianIndices)
c1, ct = Base.IteratorsMD.split(c, Val(1))
(c1, splitCartesianIndices(ct)...)
end
function splitCartesianIndices(t::Tuple{CartesianIndices, Vararg{Any}})
(splitCartesianIndices(first(t))..., splitCartesianIndices(tail(t))...)
end
function splitCartesianIndices(t::Tuple)
(first(t), splitCartesianIndices(tail(t))...)
if offset > 0 && last(r) > typemax(T) - offset
throw_upper_overflow_error()
elseif offset < 0 && first(r) < typemin(T) - offset
throw_lower_overflow_error()
end
end
splitCartesianIndices(::Tuple{}) = ()
## OffsetArray constructors

function OffsetArray(A::AbstractArray, inds::Tuple{Vararg{Union{AbstractUnitRange, CartesianIndices, Colon}}})
OffsetArray(A, splitCartesianIndices(inds))
end
function OffsetArray(A::AbstractArray, inds::Vararg{Union{AbstractUnitRange, CartesianIndices, Colon}})
OffsetArray(A, inds)
for FT in (:OffsetArray, :OffsetVector, :OffsetMatrix)
# The only route out to inner constructor
@eval function $FT(A::AbstractArray{T, N}, offsets::NTuple{N, Integer}) where {T, N}
ndims(A) == N || throw(DimensionMismatch("The number of offsets $(N) should equal ndims(A) = $(ndims(A))"))
OffsetArray{T, ndims(A), typeof(A)}(A, offsets)
end
# nested OffsetArrays
@eval $FT(A::OffsetArray{T, N}, offsets::NTuple{N, Integer}) where {T,N} = $FT(parent(A), A.offsets .+ offsets)
# convert ranges to offsets
@eval function $FT(A::AbstractArray{T}, inds::NTuple{N,OffsetAxisKnownLength}) where {T,N}
axparent = axes(A)
lA = map(length, axparent)
lI = map(length, inds)
lA == lI || throw(DimensionMismatch("supplied axes do not agree with the size of the array (got size $lA for the array and $lI for the indices"))
$FT(A, map(_offset, axparent, inds))
end
# lower CartesianIndices and Colon
@eval function $FT(A::AbstractArray{T}, inds::NTuple{N, OffsetAxis}) where {T, N}
indsN = _uncolonindices(A, _expandCartesianIndices(inds))
$FT(A, indsN)
end
@eval $FT(A::AbstractArray{T}, inds::Vararg{OffsetAxis,N}) where {T, N} = $FT(A, inds)
end

# Add methods to initialize OffsetArrays using CartesianIndices (issue #71)
function OffsetArray{T,N}(init::ArrayInitializer, inds::Tuple{Vararg{Union{AbstractUnitRange, CartesianIndices}}}) where {T,N}
indsN = stripCartesianIndices(splitCartesianIndices(inds))
OffsetArray{T,N}(init, indsN)
# array initialization
function OffsetArray{T,N}(init::ArrayInitializer, inds::NTuple{N, OffsetAxisKnownLength}) where {T,N}
AA = Array{T,N}(init, map(_indexlength, inds))
OffsetArray{T, N, typeof(AA)}(AA, map(_indexoffset, inds))
end
function OffsetArray{T}(init::ArrayInitializer, inds::Tuple{Vararg{Union{AbstractUnitRange, CartesianIndices}}}) where {T}
indsN = stripCartesianIndices(splitCartesianIndices(inds))
OffsetArray{T}(init, indsN)
function OffsetArray{T, N}(init::ArrayInitializer, inds::NTuple{NT, Union{OffsetAxisKnownLength, CartesianIndices}}) where {T, N, NT}
# NT is probably not the actual dimension of the array; CartesianIndices might contain multiple dimensions
indsN = _expandCartesianIndices(inds)
length(indsN) == N || throw(DimensionMismatch("The number of offsets $(length(indsN)) should equal ndims(A) = $N"))
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For reference,

julia> Array{Float32,3}(undef, (2, 3))
ERROR: MethodError: no method matching Array{Float32, 3}(::UndefInitializer, ::Tuple{Int64, Int64})
Closest candidates are:
  Array{T, N}(::UndefInitializer, ::Tuple{Vararg{Int64, N}}) where {T, N} at boot.jl:445
  Array{T, N}(::UndefInitializer, ::Tuple{Vararg{Integer, N}}) where {T, N} at baseext.jl:25
  Array{T, 3}(::UndefInitializer, ::Int64, ::Int64, ::Int64) where T at boot.jl:437
  ...
Stacktrace:
 [1] top-level scope
   @ REPL[1]:1

Personally I like the DimensionMismatch better. If this were to get fixed in Julia I'd suggest deferring to that error, but since it's not good now and won't be fixable on earlier releases, this seem like the right choice.

OffsetArray{T, N}(init, indsN)
end
OffsetArray{T,N}(init::ArrayInitializer, inds::Vararg{Union{AbstractUnitRange, CartesianIndices}}) where {T,N} = OffsetArray{T,N}(init, inds)
OffsetArray{T}(init::ArrayInitializer, inds::Vararg{Union{AbstractUnitRange, CartesianIndices}}) where {T} = OffsetArray{T}(init, inds)
OffsetArray{T,N}(init::ArrayInitializer, inds::Union{OffsetAxisKnownLength, CartesianIndices}...) where {T,N} = OffsetArray{T,N}(init, inds)

# avoid a level of indirection when nesting OffsetArrays
function OffsetArray(A::OffsetArray, offsets::NTuple{N,Int}) where {N}
OffsetArray(parent(A), offsets .+ A.offsets)
OffsetArray{T}(init::ArrayInitializer, inds::NTuple{N, OffsetAxisKnownLength}) where {T,N} = OffsetArray{T,N}(init, inds)
function OffsetArray{T}(init::ArrayInitializer, inds::NTuple{N, Union{OffsetAxisKnownLength, CartesianIndices}}) where {T, N}
# N is probably not the actual dimension of the array; CartesianIndices might contain multiple dimensions
indsN = _expandCartesianIndices(inds)
OffsetArray{T, length(indsN)}(init, indsN)
end
OffsetArray(A::OffsetArray{T,0}, inds::Tuple{}) where {T} = OffsetArray(parent(A), ())
# OffsetArray(A::OffsetArray{T,N}, inds::Tuple{}) where {T,N} = error("this should never be called")
OffsetArray{T}(init::ArrayInitializer, inds::Union{OffsetAxisKnownLength, CartesianIndices}...) where {T} = OffsetArray{T}(init, inds)

Base.IndexStyle(::Type{OA}) where {OA<:OffsetArray} = IndexStyle(parenttype(OA))
parenttype(::Type{OffsetArray{T,N,AA}}) where {T,N,AA} = AA
Expand All @@ -211,27 +154,24 @@ Base.eachindex(::IndexLinear, A::OffsetVector) = axes(A, 1)
@inline Base.axes(A::OffsetArray, d) = d <= ndims(A) ? IdOffsetRange(axes(parent(A), d), A.offsets[d]) : IdOffsetRange(axes(parent(A), d))
@inline Base.axes1(A::OffsetArray{T,0}) where {T} = IdOffsetRange(axes(parent(A), 1)) # we only need to specialize this one

const OffsetAxisKnownLength = Union{Integer, UnitRange, Base.OneTo, IdentityUnitRange, IdOffsetRange}

Base.similar(A::OffsetArray, ::Type{T}, dims::Dims) where T =
similar(parent(A), T, dims)
function Base.similar(A::AbstractArray, ::Type{T}, inds::Tuple{OffsetAxisKnownLength,Vararg{OffsetAxisKnownLength}}) where T
B = similar(A, T, map(indexlength, inds))
return OffsetArray(B, map(offset, axes(B), inds))
B = similar(A, T, map(_indexlength, inds))
return OffsetArray(B, map(_offset, axes(B), inds))
end

# reshape accepts a single colon
const OffsetAxis = Union{OffsetAxisKnownLength, Colon}
Base.reshape(A::AbstractArray, inds::OffsetAxis...) = reshape(A, inds)
function Base.reshape(A::AbstractArray, inds::Tuple{OffsetAxis,Vararg{OffsetAxis}})
AR = reshape(A, map(indexlength, inds))
return OffsetArray(AR, map(offset, axes(AR), inds))
AR = reshape(A, map(_indexlength, inds))
return OffsetArray(AR, map(_offset, axes(AR), inds))
end

# Reshaping OffsetArrays can "pop" the original OffsetArray wrapper and return
# an OffsetArray(reshape(...)) instead of an OffsetArray(reshape(OffsetArray(...)))
Base.reshape(A::OffsetArray, inds::Tuple{OffsetAxis,Vararg{OffsetAxis}}) =
OffsetArray(reshape(parent(A), map(indexlength, inds)), map(indexoffset, inds))
OffsetArray(reshape(parent(A), map(_indexlength, inds)), map(_indexoffset, inds))
# And for non-offset axes, we can just return a reshape of the parent directly
Base.reshape(A::OffsetArray, inds::Tuple{Union{Integer,Base.OneTo},Vararg{Union{Integer,Base.OneTo}}}) = reshape(parent(A), inds)
Base.reshape(A::OffsetArray, inds::Dims) = reshape(parent(A), inds)
Expand All @@ -241,8 +181,8 @@ Base.reshape(A::OffsetArray, inds::Union{Int,Colon}...) = reshape(parent(A), ind
Base.reshape(A::OffsetArray, inds::Tuple{Vararg{Union{Int,Colon}}}) = reshape(parent(A), inds)

function Base.similar(::Type{T}, shape::Tuple{OffsetAxis,Vararg{OffsetAxis}}) where {T<:AbstractArray}
P = T(undef, map(indexlength, shape))
OffsetArray(P, map(offset, axes(P), shape))
P = T(undef, map(_indexlength, shape))
OffsetArray(P, map(_offset, axes(P), shape))
end

Base.fill(v, inds::NTuple{N, Union{Integer, AbstractUnitRange}}) where {N} =
Expand Down Expand Up @@ -346,15 +286,6 @@ Base.pop!(A::OffsetVector) = pop!(A.parent)
Base.append!(A::OffsetVector, items) = (append!(A.parent, items); A)
Base.empty!(A::OffsetVector) = (empty!(A.parent); A)

### Low-level utilities ###

indexoffset(r::AbstractRange) = first(r) - 1
indexoffset(i::Integer) = 0
indexoffset(i::Colon) = 0
indexlength(r::AbstractRange) = length(r)
indexlength(i::Integer) = i
indexlength(i::Colon) = Colon()

# These functions keep the summary compact
function Base.inds2string(inds::Tuple{Vararg{Union{IdOffsetRange, IdentityUnitRange{<:IdOffsetRange}}}})
Base.inds2string(map(UnitRange, inds))
Expand Down
21 changes: 21 additions & 0 deletions src/utils.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,21 @@
### Low-level utilities ###

_indexoffset(r::AbstractRange) = first(r) - 1
_indexoffset(i::Integer) = 0
_indexoffset(i::Colon) = 0
_indexlength(r::AbstractRange) = length(r)
_indexlength(i::Integer) = i
_indexlength(i::Colon) = Colon()

_offset(axparent::AbstractUnitRange, ax::AbstractUnitRange) = first(ax) - first(axparent)
_offset(axparent::AbstractUnitRange, ax::CartesianIndices) = _offset(axparent, first(ax.indices))
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not sure I understand the justification here, and the fact that it's flagged by codecov is interesting.

Copy link
Member Author

@johnnychen94 johnnychen94 Sep 23, 2020 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Given that _indexlength doesn't support CartesianIndices, this single method seems not possible to be called from anywhere.

I still leave it here for potential future usage.

_offset(axparent::AbstractUnitRange, ax::Integer) = 1 - first(axparent)

_uncolonindices(A::AbstractArray{<:Any,N}, inds::NTuple{N,Any}) where {N} = _uncolonindices(axes(A), inds)
_uncolonindices(ax::Tuple, inds::Tuple) = (first(inds), _uncolonindices(tail(ax), tail(inds))...)
_uncolonindices(ax::Tuple, inds::Tuple{Colon, Vararg{Any}}) = (first(ax), _uncolonindices(tail(ax), tail(inds))...)
_uncolonindices(::Tuple{}, ::Tuple{}) = ()

_expandCartesianIndices(inds::Tuple{<:CartesianIndices, Vararg{Any}}) = (convert(Tuple{Vararg{AbstractUnitRange{Int}}}, inds[1])..., _expandCartesianIndices(Base.tail(inds))...)
_expandCartesianIndices(inds::Tuple{Any,Vararg{Any}}) = (inds[1], _expandCartesianIndices(Base.tail(inds))...)
_expandCartesianIndices(::Tuple{}) = ()
Loading