diff --git a/README.md b/README.md
index da63166..b534374 100644
--- a/README.md
+++ b/README.md
@@ -14,13 +14,13 @@ Collaboration is welcome! This is still a work-in-progress. See [the roadmap](ht
 
 ```julia
 julia> Pkg.clone("https://github.com/JuliaArrays/AxisArrays.jl")
-       using AxisArrays, SIUnits
-       import SIUnits.ShortUnits: s, ms, µs
+       using AxisArrays, Unitful
+       import Unitful: s, ms, µs
 
 julia> fs = 40000 # Generate a 40kHz noisy signal, with spike-like stuff added for testing
        y = randn(60*fs+1)*3
-       for spk = (sin(0.8:0.2:8.6) .* [0:0.01:.1; .15:.1:.95; 1:-.05:.05]   .* 50,
-                  sin(0.8:0.4:8.6) .* [0:0.02:.1; .15:.1:1; 1:-.2:.1] .* 50)
+       for spk = (sin.(0.8:0.2:8.6) .* [0:0.01:.1; .15:.1:.95; 1:-.05:.05]   .* 50,
+                  sin.(0.8:0.4:8.6) .* [0:0.02:.1; .15:.1:1; 1:-.2:.1] .* 50)
            i = rand(round(Int,.001fs):1fs)
            while i+length(spk)-1 < length(y)
                y[i:i+length(spk)-1] += spk
@@ -54,16 +54,15 @@ indices in *any* order, just so long as we annotate them with the axis name:
 
 ```jl
 julia> A[Axis{:time}(4)]
-2-dimensional AxisArray{Float64,2,...} with axes:
-    :time, 7.5e-5 s:2.5e-5 s:7.5e-5 s
-    :chan, [:c1,:c2]
-And data, a 1x2 SubArray{Float64,2,Array{Float64,2},Tuple{UnitRange{Int64},Colon},2}:
+2-dimensional AxisArray{Float64,1,...} with axes:
+    :chan, Symbol[:c1,:c2]
+And data, a 2-element Array{Float64,1}:
  -1.4144  -2.82879
 
 julia> A[Axis{:chan}(:c2), Axis{:time}(1:5)]
 1-dimensional AxisArray{Float64,1,...} with axes:
     :time, 0.0 s:2.5e-5 s:0.0001 s
-And data, a 5-element SubArray{Float64,1,Array{Float64,2},Tuple{UnitRange{Int64},Int64},2}:
+A[Axis{:chan}(:c2), Axis{:time}(1:5)]:
  -6.12181
   0.304668
  15.7366
@@ -80,7 +79,7 @@ still has the correct time information for those datapoints!
 julia> A[40µs .. 220µs, :c1]
 1-dimensional AxisArray{Float64,1,...} with axes:
     :time, 5.0e-5 s:2.5e-5 s:0.0002 s
-And data, a 7-element SubArray{Float64,1,Array{Float64,2},Tuple{UnitRange{Int64},Int64},2}:
+And data, a 7-element Array{Float64,1}:
   7.86831
  -1.4144
  -2.02881
@@ -90,7 +89,7 @@ And data, a 7-element SubArray{Float64,1,Array{Float64,2},Tuple{UnitRange{Int64}
  -1.97716
 
 julia> axes(ans, 1)
-AxisArrays.Axis{:time,SIUnits.SIRange{FloatRange{Float64},Float64,0,0,1,0,0,0,0,0,0}}(5.0e-5 s:2.5e-5 s:0.0002 s)
+AxisArrays.Axis{:time,StepRangeLen{Quantity{Float64, Dimensions:{𝐓}, Units:{s}},Base.TwicePrecision{Quantity{Float64, Dimensions:{𝐓}, Units:{s}}},Base.TwicePrecision{Quantity{Float64, Dimensions:{𝐓}, Units:{s}}}}}(5.0e-5 s:2.5e-5 s:0.0002 s)
 ```
 
 Sometimes, though, what we're really interested in is a window of time about a
@@ -125,7 +124,7 @@ julia> idxs = find(diff(A[:,:c1] .< -15) .> 0)
 julia> spks = A[atindex(-200µs .. 800µs, idxs), :c1]
 2-dimensional AxisArray{Float64,2,...} with axes:
     :time_sub, -0.000175 s:2.5e-5 s:0.000775 s
-    :time_rep, SIUnits.SIQuantity{Float64,0,0,1,0,0,0,0,0,0}[0.178725 s,0.806825 s,0.88305 s,1.47485 s,1.50465 s,1.53805 s,1.541025 s,2.16365 s,2.368425 s,2.739 s  …  57.797925 s,57.924075 s,58.06075 s,58.215125 s,58.6403 s,58.96215 s,58.990225 s,59.001325 s,59.48395 s,59.611525 s]
+    :time_rep, Quantity{Float64, Dimensions:{𝐓}, Units:{s}}[0.178725 s,0.806825 s,0.88305 s,1.47485 s,1.50465 s,1.53805 s,1.541025 s,2.16365 s,2.368425 s,2.739 s  …  57.797925 s,57.924075 s,58.06075 s,58.215125 s,58.6403 s,58.96215 s,58.990225 s,59.001325 s,59.48395 s,59.611525 s]
 And data, a 39x242 Array{Float64,2}:
  -1.53038     4.72882     5.8706    …  -0.231564      0.624714   3.44076
  -2.24961     2.12414     5.69936       7.00179       2.30993    5.20432
diff --git a/REQUIRE b/REQUIRE
index fbc022f..dbc6f97 100644
--- a/REQUIRE
+++ b/REQUIRE
@@ -2,3 +2,4 @@ julia 0.5
 IntervalSets
 Iterators
 RangeArrays
+Compat 0.19
diff --git a/src/AxisArrays.jl b/src/AxisArrays.jl
index 01a7548..f754091 100644
--- a/src/AxisArrays.jl
+++ b/src/AxisArrays.jl
@@ -4,6 +4,7 @@ module AxisArrays
 
 using Base: tail
 using RangeArrays, Iterators, IntervalSets
+using Compat
 
 export AxisArray, Axis, axisnames, axisvalues, axisdim, axes, atindex
 
diff --git a/src/combine.jl b/src/combine.jl
index 03d3a07..85f43e4 100644
--- a/src/combine.jl
+++ b/src/combine.jl
@@ -13,7 +13,7 @@ sizes{T<:AxisArray}(As::T...) = tuple(zip(map(size, As)...)...)
 matchingdims{N,T<:AxisArray}(As::NTuple{N,T}) = all(equalvalued, sizes(As...))
 matchingdimsexcept{N,T<:AxisArray}(As::NTuple{N,T}, n::Int) = all(equalvalued, sizes(As[[1:n-1; n+1:end]]...))
 
-function Base.cat{T<:AxisArray}(n::Integer, As::T...)
+function Base.cat{T}(n::Integer, As::AxisArray{T}...)
     if n <= ndims(As[1])
         matchingdimsexcept(As, n) || error("All non-concatenated axes must be identically-valued")
         newaxis = Axis{axisnames(As[1])[n]}(vcat(map(A -> A.axes[n].val, As)...))
diff --git a/src/core.jl b/src/core.jl
index 1b190d9..f6e670f 100644
--- a/src/core.jl
+++ b/src/core.jl
@@ -8,7 +8,7 @@ else
     using Base: @pure
 end
 
-typealias Symbols Tuple{Symbol,Vararg{Symbol}}
+const Symbols = Tuple{Symbol,Vararg{Symbol}}
 
 @doc """
 Type-stable axis-specific indexing and identification with a
@@ -158,17 +158,17 @@ A[ClosedInterval(0.,.3), [:a, :c]]   # select an interval and two columns
 immutable AxisArray{T,N,D,Ax} <: AbstractArray{T,N}
     data::D  # D <:AbstractArray, enforced in constructor to avoid dispatch bugs (https://github.com/JuliaLang/julia/issues/6383)
     axes::Ax # Ax<:NTuple{N, Axis}, but with specialized Axis{...} types
-    AxisArray(data::AbstractArray, axs) = new{T,N,D,Ax}(data, axs)
+    (::Type{AxisArray{T,N,D,Ax}}){T,N,D,Ax}(data::AbstractArray{T,N}, axs::Tuple{Vararg{Axis,N}}) = new{T,N,D,Ax}(data, axs)
 end
 #
 _defaultdimname(i) = i == 1 ? (:row) : i == 2 ? (:col) : i == 3 ? (:page) : Symbol(:dim_, i)
 
 default_axes(A::AbstractArray) = _default_axes(A, indices(A), ())
-_default_axes{T,N}(A::AbstractArray{T,N}, inds, axs::NTuple{N}) = axs
-@inline _default_axes{T,N,M}(A::AbstractArray{T,N}, inds, axs::NTuple{M}) =
+_default_axes{T,N}(A::AbstractArray{T,N}, inds, axs::NTuple{N,Axis}) = axs
+@inline _default_axes{T,N,M}(A::AbstractArray{T,N}, inds, axs::NTuple{M,Axis}) =
     _default_axes(A, inds, (axs..., _nextaxistype(A, axs)(inds[M+1])))
 # Why doesn't @pure work here?
-@generated function _nextaxistype{T,M}(A::AbstractArray{T}, axs::NTuple{M})
+@generated function _nextaxistype{T,M}(A::AbstractArray{T}, axs::NTuple{M,Axis})
     name = _defaultdimname(M+1)
     :(Axis{$(Expr(:quote, name))})
 end
@@ -245,7 +245,6 @@ Base.size{Ax<:Axis}(A::AxisArray, ::Type{Ax}) = size(A.data, axisdim(A, Ax))
 Base.indices(A::AxisArray) = indices(A.data)
 Base.indices(A::AxisArray, Ax::Axis) = indices(A.data, axisdim(A, Ax))
 Base.indices{Ax<:Axis}(A::AxisArray, ::Type{Ax}) = indices(A.data, axisdim(A, Ax))
-Base.linearindexing(A::AxisArray) = Base.linearindexing(A.data)
 Base.convert{T,N}(::Type{Array{T,N}}, A::AxisArray{T,N}) = convert(Array{T,N}, A.data)
 # Similar is tricky. If we're just changing the element type, it can stay as an
 # AxisArray. But if we're changing dimensions, there's no way it can know how
@@ -308,7 +307,7 @@ function permutation(to::Symbols, from::Symbols)
     li = linearindices(from)
     d = Dict(from[i]=>i for i in li)
     covered = similar(dims->falses(length(li)), li)
-    ind = Array(Int, max(n, nf))
+    ind = Array{Int}(max(n, nf))
     for (i,toi) in enumerate(to)
         j = get(d, toi, 0)
         ind[i] = j
@@ -421,7 +420,7 @@ axes(A::AbstractArray) = default_axes(A)
 axes(A::AbstractArray, dim::Int) = default_axes(A)[dim]
 
 ### Axis traits ###
-abstract AxisTrait
+@compat abstract type AxisTrait end
 immutable Dimensional <: AxisTrait end
 immutable Categorical <: AxisTrait end
 immutable Unsupported <: AxisTrait end
diff --git a/src/indexing.jl b/src/indexing.jl
index 8251e7a..ca2ca3f 100644
--- a/src/indexing.jl
+++ b/src/indexing.jl
@@ -1,9 +1,9 @@
-typealias Idx Union{Real,Colon,AbstractArray{Int}}
+const Idx = Union{Real,Colon,AbstractArray{Int}}
 
-using Base: ViewIndex, linearindexing, unsafe_getindex, unsafe_setindex!
+using Base: ViewIndex, unsafe_getindex, unsafe_setindex!
 
-# Defer linearindexing to the wrapped array
-Base.linearindexing{T,N,D}(::AxisArray{T,N,D}) = linearindexing(D)
+# Defer IndexStyle to the wrapped array
+@compat Base.IndexStyle{T,N,D,Ax}(::Type{AxisArray{T,N,D,Ax}}) = IndexStyle(D)
 
 # Simple scalar indexing where we just set or return scalars
 @inline Base.getindex(A::AxisArray, idxs::Int...) = A.data[idxs...]
@@ -26,7 +26,8 @@ Base.setindex!(A::AxisArray, v, idx::Base.IteratorsMD.CartesianIndex) = (A.data[
     end
     names = axisnames(A)
     newaxes = Expr[]
-    for d=1:lastnonscalar-droplastaxis
+    drange = 1:lastnonscalar-droplastaxis
+    for d=drange
         if I[d] <: AxisArray
             # Indexing with an AxisArray joins the axis names
             idxnames = axisnames(I[d])
@@ -34,8 +35,15 @@ Base.setindex!(A::AxisArray, v, idx::Base.IteratorsMD.CartesianIndex) = (A.data[
                 push!(newaxes, :($(Axis{Symbol(names[d], "_", idxnames[i])})(I[$d].axes[$i].val)))
             end
         elseif I[d] <: Real
-        elseif I[d] <: Union{AbstractVector,Colon}
+        elseif I[d] <: AbstractVector
             push!(newaxes, :($(Axis{names[d]})(A.axes[$d].val[Base.to_index(I[$d])])))
+        elseif I[d] <: Colon
+            if d < length(I) || d <= ndims(A)
+                push!(newaxes, :($(Axis{names[d]})(A.axes[$d].val)))
+            else
+                dimname = _defaultdimname(d)
+                push!(newaxes, :($(Axis{dimname})(Base.OneTo(Base.trailingsize(A, $d)))))
+            end
         elseif I[d] <: AbstractArray
             for i=1:ndims(I[d])
                 # When we index with non-vector arrays, we *add* dimensions.
@@ -117,6 +125,12 @@ end
     return :($meta; to_index(A, $(idxs...)))
 end
 
+function Base.reshape{N}(A::AxisArray, ::Type{Val{N}})
+    # axN, _ = Base.IteratorsMD.split(axes(A), Val{N})
+    # AxisArray(reshape(A.data, Val{N}), reaxis(A, Base.fill_to_length(axN, :, Val{N})...))
+    AxisArray(reshape(A.data, Val{N}), reaxis(A, ntuple(d->Colon(), Val{N})...))
+end
+
 ### Indexing along values of the axes ###
 
 # Default axes indexing throws an error
diff --git a/src/intervals.jl b/src/intervals.jl
index 8595d02..68b55e3 100644
--- a/src/intervals.jl
+++ b/src/intervals.jl
@@ -17,7 +17,7 @@
 # downside is that Intervals are not as useful as they could be; they really
 # could be considered as <: Number themselves. We do this in general for any
 # supported Scalar:
-typealias Scalar Union{Number, Dates.AbstractTime}
+const Scalar = Union{Number, Dates.AbstractTime}
 Base.promote_rule{T<:Scalar}(::Type{ClosedInterval{T}}, ::Type{T}) = ClosedInterval{T}
 Base.promote_rule{T,S<:Scalar}(::Type{ClosedInterval{T}}, ::Type{S}) = ClosedInterval{promote_type(T,S)}
 Base.promote_rule{T,S}(::Type{ClosedInterval{T}}, ::Type{ClosedInterval{S}}) = ClosedInterval{promote_type(T,S)}
@@ -62,7 +62,7 @@ immutable RepeatedInterval{T,S,A} <: AbstractVector{T}
 end
 RepeatedInterval{S,A<:AbstractVector}(window::ClosedInterval{S}, offsets::A) = RepeatedInterval{promote_type(ClosedInterval{S}, eltype(A)), S, A}(window, offsets)
 Base.size(r::RepeatedInterval) = size(r.offsets)
-Base.linearindexing{R<:RepeatedInterval}(::Type{R}) = Base.LinearFast()
+@compat Base.IndexStyle(::Type{<:RepeatedInterval}) = IndexLinear()
 Base.getindex(r::RepeatedInterval, i::Int) = r.window + r.offsets[i]
 +(window::ClosedInterval, offsets::AbstractVector) = RepeatedInterval(window, offsets)
 +(offsets::AbstractVector, window::ClosedInterval) = RepeatedInterval(window, offsets)
@@ -84,5 +84,5 @@ immutable RepeatedIntervalAtIndexes{T,A<:AbstractVector{Int}} <: AbstractVector{
 end
 atindex(window::ClosedInterval, indexes::AbstractVector) = RepeatedIntervalAtIndexes(window, indexes)
 Base.size(r::RepeatedIntervalAtIndexes) = size(r.indexes)
-Base.linearindexing{R<:RepeatedIntervalAtIndexes}(::Type{R}) = Base.LinearFast()
+@compat Base.IndexStyle(::Type{<:RepeatedIntervalAtIndexes}) = IndexLinear()
 Base.getindex(r::RepeatedIntervalAtIndexes, i::Int) = IntervalAtIndex(r.window, r.indexes[i])
diff --git a/src/search.jl b/src/search.jl
index 22d11d6..e2cca18 100644
--- a/src/search.jl
+++ b/src/search.jl
@@ -36,21 +36,30 @@ function Base.searchsorted(a::Range, I::ClosedInterval)
     searchsortedfirst(a, I.left):searchsortedlast(a, I.right)
 end
 
-if VERSION > v"0.5.0-dev+4557"
-    # When running with "--check-bounds=yes" (like on Travis), the bounds-check isn't elided
-    @inline function Base.unsafe_getindex{T}(v::Range{T}, i::Integer)
-        convert(T, first(v) + (i-1)*step(v))
-    end
+# When running with "--check-bounds=yes" (like on Travis), the bounds-check isn't elided
+@inline function Base.unsafe_getindex{T}(v::Range{T}, i::Integer)
+    convert(T, first(v) + (i-1)*step(v))
+end
+@inline function Base.unsafe_getindex{T<:Integer}(r::StepRange, s::Range{T})
+    st = oftype(r.start, r.start + (first(s)-1)*step(r))
+    range(st, step(r)*step(s), length(s))
+end
+if VERSION < v"0.6.0-dev.2390"
+    include_string("""
     @inline function Base.unsafe_getindex{T}(r::FloatRange{T}, i::Integer)
         convert(T, (r.start + (i-1)*r.step)/r.divisor)
     end
-    @inline function Base.unsafe_getindex{T<:Integer}(r::StepRange, s::Range{T})
-        st = oftype(r.start, r.start + (first(s)-1)*step(r))
-        range(st, step(r)*step(s), length(s))
-    end
     @inline function Base.unsafe_getindex(r::FloatRange, s::OrdinalRange)
         FloatRange(r.start + (first(s)-1)*r.step, step(s)*r.step, length(s), r.divisor)
     end
+    """)
+else
+    include_string("""
+    @inline function Base.unsafe_getindex(r::StepRangeLen, s::OrdinalRange)
+        vfirst = unsafe_getindex(r, first(s))
+        StepRangeLen(vfirst, r.step*step(s), length(s))
+    end
+    """)
 end
 
 function unsafe_searchsortedlast{T<:Number}(a::Range{T}, x::Number)
diff --git a/test/REQUIRE b/test/REQUIRE
index bcad44a..0082847 100644
--- a/test/REQUIRE
+++ b/test/REQUIRE
@@ -1,2 +1,2 @@
 OffsetArrays
-SIUnits
+Unitful
diff --git a/test/indexing.jl b/test/indexing.jl
index bb357ad..2f9d753 100644
--- a/test/indexing.jl
+++ b/test/indexing.jl
@@ -37,8 +37,13 @@ D[1,1,1,1,1] = 10
 
 # Linear indexing across multiple dimensions drops tracking of those dims
 @test A[:].axes[1].val == 1:length(A)
-@test A[1:2,:].axes[1].val == A.axes[1].val[1:2]
-@test A[1:2,:].axes[2].val == 1:Base.trailingsize(A,2)
+B = A[1:2,:]
+@test B.axes[1].val == A.axes[1].val[1:2]
+@test B.axes[2].val == 1:Base.trailingsize(A,2)
+B2 = reshape(A, Val{2})
+B = B2[1:2,:]
+@test B.axes[1].val == A.axes[1].val[1:2]
+@test B.axes[2].val == 1:Base.trailingsize(A,2)
 
 B = AxisArray(reshape(1:15, 5,3), .1:.1:0.5, [:a, :b, :c])
 
diff --git a/test/readme.jl b/test/readme.jl
index 67fe7ca..053468e 100644
--- a/test/readme.jl
+++ b/test/readme.jl
@@ -1,12 +1,12 @@
 # Intended to ensure the README stays working (this is a copy)
 
-using AxisArrays, SIUnits
-import SIUnits.ShortUnits: s, ms, µs
+using AxisArrays, Unitful
+import Unitful: s, ms, µs
 
 fs = 40000
 y = randn(60*fs+1)*3
-for spk = (sin(0.8:0.2:8.6) .* [0:0.01:.1; .15:.1:.95; 1:-.05:.05]   .* 50,
-           sin(0.8:0.4:8.6) .* [0:0.02:.1; .15:.1:1; 1:-.2:.1] .* 50)
+for spk = (sin.(0.8:0.2:8.6) .* [0:0.01:.1; .15:.1:.95; 1:-.05:.05]   .* 50,
+           sin.(0.8:0.4:8.6) .* [0:0.02:.1; .15:.1:1; 1:-.2:.1] .* 50)
     i = rand(round(Int,.001fs):1fs)
     while i+length(spk)-1 < length(y)
         y[i:i+length(spk)-1] += spk
diff --git a/test/runtests.jl b/test/runtests.jl
index 1167722..776a138 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,7 +1,9 @@
 using AxisArrays
 using Base.Test
 
-@test isempty(detect_ambiguities(AxisArrays, Base, Core))
+if VERSION < v"0.6.0-dev"
+    @test isempty(detect_ambiguities(AxisArrays, Base, Core))
+end
 
 include("core.jl")
 include("intervals.jl")