adding slicedim() function

adding intersect() on Ranges
adding squeeze()

Author: JeffBezanson

j/abstractarray.j

@@ -464,6 +464,43 @@ function ref(A::AbstractArray, I::Indices...)
 end
 end
 
+# index A[:,:,...,i,:,:,...] where "i" is in dimension "d"
+# TODO: more optimized special cases
+function slicedim(A::AbstractArray, d::Int, i)
+    A[ntuple(ndims(A), n->(n==d ? i : (1:size(A,n))))...]
+end
+
+function slicedim(A::AbstractArray, d::Int, i::Int)
+    d_in = size(A)
+    leading = d_in[1:(d-1)]
+    d_out = append(leading, (1,), d_in[(d+1):end])
+
+    M = prod(leading)
+    N = numel(A)
+    stride = M * d_in[d]
+
+    B = similar(A, d_out)
+    index_offset = 1 + (i-1)*M
+
+    l = 1
+
+    if M==1
+        for j=0:stride:(N-stride)
+            B[l] = A[j + index_offset]
+            l += 1
+        end
+    else
+        for j=0:stride:(N-stride)
+            offs = j + index_offset
+            for k=0:(M-1)
+                B[l] = A[offs + k]
+                l += 1
+            end
+        end
+    end
+    return B
+end
+
 ## Indexing: assign ##
 
 # 1-d indexing is assumed defined on subtypes
@@ -1317,6 +1354,16 @@ function ind2sub(dims, ind::Int)
     return tuple(ind, sub...)
 end
 
+function squeeze(A::AbstractArray)
+    d = ()
+    for i = size(A)
+        if i != 1
+            d = tuple(d..., i)
+        end
+    end
+    reshape(A, d)
+end
+
 ## subarrays ##
 
 type SubArray{T,N,A<:AbstractArray,I<:(RangeIndex...,)} <: AbstractArray{T,N}

j/darray.j

@@ -89,29 +89,14 @@ function localize(s::SubDArray)
     l = localize(d)
     if isa(sdi,Int)
         if lo <= sdi <= hi
-            return l[ntuple(ndims(l), i->(i==d.distdim ? sdi-lo+1 :
-                                          1:size(l,i)))...]
+            return slicedim(l, d.distdim, sdi-lo+1)
         else
-            return Array(eltype(l), ntuple(ndims(l),
-                                           i->(i==d.distdim ? 0 :
-                                               size(l,i))))
+            return Array(eltype(l), ntuple(ndims(l), i->(i==d.distdim ? 0 :
+                                                         size(l,i))))
         end
     else
-        sta = start(sdi)
-        ste = step(sdi)
-        sto = sdi[end]
-        i0 = sta + ste*div((lo-sta), ste)
-        i1 = sta + ste*div((hi-sta), ste)
-        i0 = max(i0,sta)
-        i1 = min(i1,sto)
-        if i0 > sto || i1 < sta
-            return Array(eltype(l), ntuple(ndims(l),
-                                           i->(i==d.distdim ? 0 :
-                                               size(l,i))))
-        else
-            return l[ntuple(ndims(l), i->(i==d.distdim ? i0:ste:i1 :
-                                          1:size(l,i)))...]
-        end
+        r = intersect(lo:hi, sdi)
+        return l[ntuple(ndims(l), i->(i==d.distdim ? r : 1:size(l,i)))...]
     end
 end
 

j/range.j

@@ -24,11 +24,6 @@ similar(r::Range1, T::Type, dims::Dims) = Range1(convert(T, r.start), convert(T,
 
 typealias Ranges Union(Range,Range1)
 
-ref(r::Range, s::Range{Index}) = Range(r[s[1]],r.step*s.step,r[s[end]])
-ref(r::Range1, s::Range{Index}) = Range(r[s[1]],s.step,r[s[end]])
-ref(r::Range, s::Range1{Index}) = Range(r[s[1]],r.step,r[s[end]])
-ref(r::Range1, s::Range1{Index}) = Range1(r[s[1]],r[s[end]])
-
 step(r::Range)  = r.step
 step(r::Range1) = one(r.start)
 
@@ -63,6 +58,11 @@ next{T}(r::Range{T}, st) = (r.start+st*r.step, st+1)
 colon(start::Real, stop::Real, step::Real) = Range(start, step, stop)
 colon(start::Real, stop::Real) = Range1(start, stop)
 
+ref(r::Range, s::Range{Index}) = Range(r[s[1]],r.step*s.step,r[s[end]])
+ref(r::Range1, s::Range{Index}) = Range(r[s[1]],s.step,r[s[end]])
+ref(r::Range, s::Range1{Index}) = Range(r[s[1]],r.step,r[s[end]])
+ref(r::Range1, s::Range1{Index}) = Range1(r[s[1]],r[s[end]])
+
 function ref(r::Ranges, i::Int)
     if i < 1; error(BoundsError); end
     x = r.start + (i-1)*step(r)
@@ -72,6 +72,23 @@ function ref(r::Ranges, i::Int)
     return x
 end
 
+intersect(r::Range1, s::Range1) = max(r.start,s.start):min(r.stop,s.stop)
+
+intersect(r::Range, s::Range1) = intersect(s, r)
+
+function intersect(r::Range1, s::Range)
+    sta = start(s)
+    ste = step(s)
+    sto = s[end]
+    lo = r.start
+    hi = r.stop
+    i0 = max(lo, sta + ste*div((lo-sta)+ste-1, ste))
+    i1 = min(hi, sta + ste*div((hi-sta), ste))
+    i0 = max(i0, sta)
+    i1 = min(i1, sto)
+    i0:ste:i1
+end
+
 ## linear operations on 1-d ranges ##
 
 -(r::Ranges) = Range(-r.start, -step(r), -r.stop)