Merge pull request #423 from gridap/kernel_refactoring_fv_comments

Make Reindex and PosNegReindex resilient to empty array of values

src/Arrays/Arrays.jl

@@ -89,6 +89,7 @@ export Reindex
 #export reindex
 
 export PosNegReindex
+export PosNegPartition
 #export posneg_reindex
 
 export FilterMap

src/Arrays/CompressedArrays.jl

@@ -53,13 +53,13 @@ function IndexStyle(a::Type{CompressedArray{T,N,A,P}}) where {T,N,A,P}
   IndexStyle(P)
 end
 
-function lazy_map(::typeof(evaluate),g::CompressedArray...)
-  if _have_same_ptrs(g)
-    _lazy_map_compressed(g...)
-  else
-    LazyArray(g...)
-  end
-end
+#function lazy_map(::typeof(evaluate),g::CompressedArray...)
+#  if _have_same_ptrs(g)
+#    _lazy_map_compressed(g...)
+#  else
+#    LazyArray(g...)
+#  end
+#end
 
 function lazy_map(::typeof(evaluate),::Type{T},g::CompressedArray...) where T
   if _have_same_ptrs(g)
@@ -69,16 +69,16 @@ function lazy_map(::typeof(evaluate),::Type{T},g::CompressedArray...) where T
   end
 end
 
-function lazy_map(::typeof(evaluate),g::Union{CompressedArray,Fill}...)
-  g_compressed = _find_compressed_ones(g)
-  if _have_same_ptrs(g_compressed)
-    g1 = first(g_compressed)
-    g_all_compressed = map(gi->_compress(gi,g1),g)
-    _lazy_map_compressed(g_all_compressed...)
-  else
-    LazyArray(g...)
-  end
-end
+#function lazy_map(::typeof(evaluate),g::Union{CompressedArray,Fill}...)
+#  g_compressed = _find_compressed_ones(g)
+#  if _have_same_ptrs(g_compressed)
+#    g1 = first(g_compressed)
+#    g_all_compressed = map(gi->_compress(gi,g1),g)
+#    _lazy_map_compressed(g_all_compressed...)
+#  else
+#    LazyArray(g...)
+#  end
+#end
 
 function lazy_map(::typeof(evaluate),::Type{T},g::Union{CompressedArray,Fill}...) where T
   g_compressed = _find_compressed_ones(g)

src/Arrays/Interface.jl

@@ -135,6 +135,7 @@ bi = testitem(b)
 
 """
 function testitem(a::AbstractArray{T}) where T
+  #@check isconcretetype(T) "This array is type-instable"
   if length(a) >0
     first(a)
   else

src/Arrays/LazyArrays.jl

@@ -53,14 +53,16 @@ println(c)
 ```
 """
 @inline function lazy_map(k,f::AbstractArray...)
-  s = _common_size(f...)
-  lazy_map(evaluate,Fill(k, s), f...)
+  fi = map(testitem,f)
+  T = return_type(k, fi...)
+  lazy_map(k,T,f...)
 end
 
-@inline lazy_map(::typeof(evaluate),k::AbstractArray,f::AbstractArray...) = LazyArray(k,f...)
+#@inline lazy_map(::typeof(evaluate),k::AbstractArray,f::AbstractArray...) = LazyArray(k,f...)
 
+# This is the function to be overload to specialize on the Map f
 """
-    lazy_map(::Type{T},f,a::AbstractArray...) where T
+    lazy_map(f,::Type{T},a::AbstractArray...) where T
 
 Like [`lazy_map(f,a::AbstractArray...)`](@ref), but the user provides the element type
 of the resulting array in order to circumvent type inference.
@@ -70,6 +72,7 @@ of the resulting array in order to circumvent type inference.
   lazy_map(evaluate,T,Fill(k, s), f...)
 end
 
+# This is the function to be overload to specialize on the array types
 @inline lazy_map(::typeof(evaluate),T::Type,k::AbstractArray,f::AbstractArray...) = LazyArray(T,k,f...)
 
 """
@@ -94,12 +97,12 @@ struct LazyArray{G,T,N,F} <: AbstractArray{T,N}
   end
 end
 
-function LazyArray(g::AbstractArray{S}, f::AbstractArray...) where S
-  isconcretetype(S) ? gi = testitem(g) : @notimplemented
-  fi = map(testitem,f)
-  T = return_type(gi, fi...)
-  LazyArray(T, g, f...)
-end
+#function LazyArray(g::AbstractArray{S}, f::AbstractArray...) where S
+#  isconcretetype(S) ? gi = testitem(g) : @notimplemented
+#  fi = map(testitem,f)
+#  T = return_type(gi, fi...)
+#  LazyArray(T, g, f...)
+#end
 
 IndexStyle(::Type{<:LazyArray}) = IndexCartesian()
 
@@ -138,7 +141,8 @@ function _array_cache!(hash::Dict,a::LazyArray)
   cf = map(fi->array_cache(hash,fi),a.f)
   cgi = return_cache(gi, fi...)
   index = -1
-  item = evaluate!(cgi,gi,testargs(gi,fi...)...)
+  #item = evaluate!(cgi,gi,testargs(gi,fi...)...)
+  item = return_value(gi,fi...)
   (cg, cgi, cf), IndexItemPair(index, item)
 end
 
@@ -185,17 +189,21 @@ Base.size(a::LazyArray) = size(a.g)
 
 # Particular implementations for Fill
 
-function lazy_map(::typeof(evaluate),f::Fill, a::Fill...)
+#function lazy_map(::typeof(evaluate),f::Fill, a::Fill...)
+#  ai = map(ai->ai.value,a)
+#  r = evaluate(f.value, ai...)
+#  s = _common_size(f, a...)
+#  Fill(r, s)
+#end
+
+function lazy_map(::typeof(evaluate),::Type{T}, f::Fill, a::Fill...) where T
+  #lazy_map(evaluate, f, a...)
   ai = map(ai->ai.value,a)
   r = evaluate(f.value, ai...)
   s = _common_size(f, a...)
   Fill(r, s)
 end
 
-function lazy_map(::typeof(evaluate),::Type{T}, f::Fill, a::Fill...) where T
-  lazy_map(evaluate, f, a...)
-end
-
 function _common_size(a::AbstractArray...)
   a1, = a
   @check all(map(ai->length(a1) == length(ai),a)) "Array sizes $(map(size,a)) are not compatible."

src/Arrays/Maps.jl

@@ -51,7 +51,9 @@ evaluate!(cache,f,x...) = @abstractmethod
 Returns the type of the result of calling mapping `f` with
 arguments of the types of the objects `x`.
 """
-return_type(f,x...) = typeof(evaluate(f,testargs(f,x...)...))
+return_type(f,x...) = typeof(return_value(f,x...))
+
+return_value(f,x...) = evaluate(f,testargs(f,x...)...)
 
 """
     testargs(f,x...)

src/Arrays/PosNegReindex.jl

@@ -71,9 +71,19 @@ end
 
 function testargs(k::PosNegReindex,i::Integer)
   @check length(k.values_pos) !=0 || length(k.values_neg) != 0 "This map has empty domain"
+  @check eltype(k.values_pos) == eltype(k.values_neg) "This map is type-instable"
   length(k.values_pos) !=0 ? (one(i),) : (-one(i))
 end
 
+function return_value(k::PosNegReindex,i::Integer)
+  if length(k.values_pos)==0 && length(k.values_neg)==0
+    @check eltype(k.values_pos) == eltype(k.values_neg) "This map is type-instable"
+    testitem(k.values_pos)
+  else
+    evaluate(k,testargs(k,i)...)
+  end
+end
+
 function return_cache(k::PosNegReindex,i::Integer)
   c_p = array_cache(k.values_pos)
   c_n = array_cache(k.values_neg)
@@ -85,6 +95,187 @@ end
   i>0 ? getindex!(c_p,k.values_pos,i) : getindex!(c_n,k.values_neg,-i)
 end
 
+#function lazy_map(::typeof(evaluate),a::LazyArray{<:Fill{<:PosNegReindex}}...)
+#  i_to_iposneg = a[1].f[1]
+#  if all(map( ai-> is_exhaustive(a[1].f[1]),a)) && all( map( ai-> i_to_iposneg==a[1].f[1],a) )
+#    bpos = map(ai->ai.g.value.values_pos,a)
+#    bneg = map(ai->ai.g.value.values_neg,a)
+#    cpos = lazy_map(evaluate,bpos...)
+#    cneg = lazy_map(evaluate,bneg...)
+#    lazy_map(PosNegReindex(cpos,cneg),i_to_iposneg)
+#  else
+#    LazyArray(a...)
+#  end
+#end
+
+function lazy_map(::typeof(evaluate),::Type{T},a::LazyArray{<:Fill{<:PosNegReindex}}...) where T
+  i_to_iposneg = a[1].f[1]
+  if all(map( ai-> is_exhaustive(a[1].f[1]),a)) && all( map( ai-> i_to_iposneg==a[1].f[1],a) )
+    bpos = map(ai->ai.g.value.values_pos,a)
+    bneg = map(ai->ai.g.value.values_neg,a)
+    cpos = lazy_map(evaluate,T,bpos...)
+    cneg = lazy_map(evaluate,T,bneg...)
+    lazy_map(PosNegReindex(cpos,cneg),T,i_to_iposneg)
+  else
+    LazyArray(T,a...)
+  end
+end
+
+#function lazy_map(::typeof(evaluate),b::Fill,a::LazyArray{<:Fill{<:PosNegReindex}}...)
+#  i_to_iposneg = a[1].f[1]
+#  if all(map( ai-> is_exhaustive(a[1].f[1]),a)) && all( map( ai-> i_to_iposneg==a[1].f[1],a) )
+#    k = b.value
+#    bpos = map(ai->ai.g.value.values_pos,a)
+#    bneg = map(ai->ai.g.value.values_neg,a)
+#    cpos = lazy_map(k,bpos...)
+#    cneg = lazy_map(k,bneg...)
+#    lazy_map(PosNegReindex(cpos,cneg),i_to_iposneg)
+#  else
+#    LazyArray(b,a...)
+#  end
+#end
+
+function lazy_map(::typeof(evaluate),::Type{T},b::Fill,a::LazyArray{<:Fill{<:PosNegReindex}}...) where T
+  i_to_iposneg = a[1].f[1]
+  if all(map( ai-> is_exhaustive(a[1].f[1]),a)) && all( map( ai-> i_to_iposneg==a[1].f[1],a) )
+    k = b.value
+    bpos = map(ai->ai.g.value.values_pos,a)
+    bneg = map(ai->ai.g.value.values_neg,a)
+    cpos = lazy_map(k,T,bpos...)
+    cneg = lazy_map(k,T,bneg...)
+    lazy_map(PosNegReindex(cpos,cneg),T,i_to_iposneg)
+  else
+    LazyArray(T,b,a...)
+  end
+end
+
+# Helper functions to work with arrays representing a binary partition
+
+function pos_and_neg_length(i_to_iposneg)
+  Npos = maximum(i_to_iposneg)
+  Nneg = -minimum(i_to_iposneg)
+  Npos, Nneg
+end
+
+function is_exhaustive(i_to_iposneg)
+  Npos, Nneg = pos_and_neg_length(i_to_iposneg)
+  if length(i_to_iposneg) != Npos+Nneg
+    return false
+  end
+  ipos_to_touched = fill(false,Npos)
+  ineg_to_touched = fill(false,Nneg)
+  for iposneg in i_to_iposneg
+    iposneg>0 ? ipos_to_touched[iposneg] = true : ineg_to_touched[-iposneg] = true
+  end
+  all(ipos_to_touched) && all(ineg_to_touched)
+end
+
+function pos_and_neg_indices(i_to_iposneg)
+  @check is_exhaustive(i_to_iposneg)
+  Npos, Nneg = pos_and_neg_length(i_to_iposneg)
+  ipos_to_i = zeros(Int,Npos)
+  ineg_to_i = zeros(Int,Nneg)
+  @inbounds for (i,iposneg) in enumerate(i_to_iposneg)
+    iposneg>0 ? ipos_to_i[iposneg] = i : ineg_to_i[-iposneg] = i
+  end
+  ipos_to_i, ineg_to_i
+end
+
+function aligned_with_pos(i_to_iposneg,j_to_i,npos)
+  j_to_iposneg = lazy_map(Reindex(i_to_iposneg),j_to_i)
+  j_to_iposneg == 1:npos
+end
+
+function aligned_with_neg(i_to_iposneg,j_to_i,nneg)
+  j_to_iposneg = lazy_map(Reindex(i_to_iposneg),j_to_i)
+  j_to_iposneg == -(1:nneg)
+end
+
+function all_pos(i_to_iposneg,j_to_i)
+  j_to_iposneg = lazy_map(Reindex(i_to_iposneg),j_to_i)
+  all( lazy_map(iposneg -> iposneg>0, j_to_iposneg) )
+end
+
+function all_neg(i_to_iposneg,j_to_i)
+  j_to_iposneg = lazy_map(Reindex(i_to_iposneg),j_to_i)
+  all( lazy_map(iposneg -> iposneg<0, j_to_iposneg) )
+end
+
+# This is important to do optimizations associated with ExtendedFESpace
+"""
+struct representing a binary partition of a range of indices
+
+Using this allows one to do a number of important optimizations when working with `PosNegReindex`
+"""
+struct PosNegPartition{T,V,Vp,Vn} <: AbstractVector{T}
+  i_to_iposneg::V
+  ipos_to_i::Vp
+  ineg_to_i::Vn
+  function PosNegPartition(
+    i_to_iposneg::AbstractVector,
+    ipos_to_i::AbstractVector,
+    ineg_to_i::AbstractVector)
+
+    @check eltype(ipos_to_i) == eltype(ineg_to_i)
+    @check all((ipos_to_i, ineg_to_i) .== pos_and_neg_indices(i_to_iposneg))
+    T = eltype(i_to_iposneg)
+    @check T <: Integer
+    V = typeof(i_to_iposneg)
+    Vp = typeof(ipos_to_i)
+    Vn = typeof(ineg_to_i)
+    new{T,V,Vp,Vn}(i_to_iposneg,ipos_to_i,ineg_to_i)
+  end
+end
+
+function PosNegPartition(ipos_to_i::AbstractArray,Ni::Integer)
+  i_to_iposneg = zeros(typeof(Ni),Ni)
+  i_to_iposneg[ipos_to_i] .= 1:length(ipos_to_i)
+  Nneg = count(k->k==0,i_to_iposneg)
+  ineg_to_i = similar(ipos_to_i,eltype(ipos_to_i),(Nneg,))
+  ineg = 1
+  for (i,iposneg) in enumerate(i_to_iposneg)
+    if iposneg==0
+      i_to_iposneg[i] = -ineg
+      ineg_to_i[ineg] = i
+      ineg += 1
+    end
+  end
+  PosNegPartition(i_to_iposneg,ipos_to_i,ineg_to_i)
+end
+
+Base.IndexStyle(::Type{<:PosNegPartition}) = IndexLinear()
+@propagate_inbounds Base.getindex(a::PosNegPartition,i::Integer) = a.i_to_iposneg[i]
+Base.size(a::PosNegPartition) = size(a.i_to_iposneg)
+get_array(a::PosNegPartition) = a.i_to_iposneg
+
+pos_and_neg_length(a::PosNegPartition) = (length(a.ipos_to_i),length(a.ineg_to_i))
+is_exhaustive(a::PosNegPartition) = true
+pos_and_neg_indices(a::PosNegPartition) = (a.ipos_to_i,a.ineg_to_i)
+function aligned_with_pos(a::PosNegPartition,j_to_i,npos)
+  @check length(a.ipos_to_i) == npos
+  a.ipos_to_i === j_to_i || a.ipos_to_i == j_to_i
+end
+function aligned_with_neg(a::PosNegPartition,j_to_i,nneg)
+  @check length(a.ineg_to_i) == nneg
+  a.ineg_to_i === j_to_i || a.ineg_to_i == j_to_i
+end
+
+#function lazy_map(::typeof(evaluate),a::LazyArray{<:Fill{<:PosNegReindex}},b::AbstractArray...)
+#  bpos = map(bi->lazy_map(Reindex(),bi),b)
+#
+#end
+#
+#i_to_iposneg
+#
+#Ni = length(i_to_iposneg)
+#Nipos = count(iposneg->iposneg>0,i_to_iposneg)
+#Nineg = Ni - Nipos
+#ipos_to_i = zeros(Int,Nipos)
+#ineg_to_i = zeros(Int,Nipos)
+#for (i,iposneg) in enumerate(i_to_iposneg)
+#  iposneg>0 ? ipos_to_i[iposneg]=i : ineg_to_i[-iposneg]=i
+#end
+
 
 #@inline return_type(k::PosNegReindex,x...) = typeof(first(k.values))
 #