Merge pull request #4 from FelixBenning/1-lazy_product

Lazy Product

Author: FelixBenning

Project.toml

@@ -9,6 +9,7 @@ KernelFunctions = "ec8451be-7e33-11e9-00cf-bbf324bd1392"
 LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MappedArrays = "dbb5928d-eab1-5f90-85c2-b9b0edb7c900"
+ProductArrays = "39a08dd0-b4b7-4086-afbb-eb1796240846"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 
 [compat]

src/multiOutput.jl

@@ -1,39 +1,7 @@
 
 using MappedArrays: mappedarray, ReadonlyMappedArray
 
-function ensure_all_linear_indexed(vecs::T) where {T<:Tuple}
-    linear_indexed = ntuple(
-        n -> Base.IndexStyle(fieldtype(T,n)) === IndexLinear(),
-        Base._counttuple(T)
-    )
-    all(linear_indexed) || throw(ArgumentError(
-        "$(vecs[findfirst(x->!x, linear_indexed)]) cannot be linearly accessed. All inputs need to implement `Base.getindex(::T, ::Int)`"
-    ))
-end
-
-struct ProductArray
-    arrays
-end
-
-"""
-    lazy_product(vectors...)
-
-The output is a lazy form of
-```julia
-collect(Iterators.product(vectors...))
-```
-i.e. it is an AbstractArray in contrast to `Iterators.product(vectors...)`. So
-is accessible with `getindex` and gets default Array implementations for free.
-In particular it can be passed to `Base.PermutedDimsArray`` for lazy permutation
-and `vec()` to obtain a lazy `Base.ReshapedArray`.
-"""
-function lazy_product(vectors...)
-    ensure_all_linear_indexed(vectors)
-    indices = CartesianIndices(ntuple(n -> axes(vec(vectors[n]), 1), length(vectors)))
-    return mappedarray(indices) do idx
-        return ntuple(n -> vec(vectors[n])[idx[n]], length(vectors))
-    end
-end
+using ProductArrays: productArray
 
 """
     lazy_flatten(vectors...)
@@ -49,19 +17,19 @@ and `vec()` to obtain a lazy `Base.ReshapedArray`.
 """
 function lazy_flatten(vectors...)
     ensure_all_linear_indexed(vectors)
-    lengths = cumsum(length.(vectors)) 
+    lengths = cumsum(length.(vectors))
     return mappedarray(1:lengths[end]) do idx
         # this is not efficient for iteration - maybe go with LazyArrays.jl -> Vcat instead.
         v_idx = searchsortedfirst(lengths, idx)
-        return vectors[v_idx][idx - get(lengths, v_idx-1, 0)]
+        return vectors[v_idx][idx-get(lengths, v_idx - 1, 0)]
     end
 end
 
 
 function multi_out_byFeatures(positions, out_dims)
-    return vec(PermutedDimsArray(lazy_product(positions, 1:out_dims), (2, 1)))
+    return vec(PermutedDimsArray(productArray(positions, 1:out_dims), (2, 1)))
 end
 function multi_out_byOutput(positions, out_dims)
-    return vec(lazy_product(positions, 1:out_dims))
+    return vec(productArray(positions, 1:out_dims))
 end
 

src/partial.jl

@@ -55,8 +55,8 @@ end
 const DiffPt{T} = Tuple{T,Partial}
 
 gradient(dim::Integer) = mappedarray(partial, Base.OneTo(dim))
-hessian(dim::Integer) = mappedarray(partial, lazy_product(Base.OneTo(dim), Base.OneTo(dim)))
-fullderivative(order::Integer,dim::Integer) = mappedarray(partial, lazy_product(ntuple(_->Base.OneTo(dim), order)...))
+hessian(dim::Integer) = mappedarray(partial, productArray(Base.OneTo(dim), Base.OneTo(dim)))
+fullderivative(order::Integer,dim::Integer) = mappedarray(partial, productArray(ntuple(_->Base.OneTo(dim), order)...))
 
 # idea: lazy mappings can be undone (extract original range -> towards a specialization speedup of broadcasting over multiple derivatives using backwardsdiff)
 const MappedPartialVec{T} = ReadonlyMappedArray{Partial{1,Int},1,T,typeof(partial)}

test/multiOutput.jl

@@ -6,7 +6,6 @@ using Test
 List of Testfiles without extension. `\$(test).jl"` should be a file for every test in AVAILABLE_TESTS
 """
 const AVAILABLE_TESTS = [
-    "multiOutput",
     "diffKernel",
 ]